|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00013 MODULE mc_types 00014 USE atomic_kind_types, ONLY: atomic_kind_type,& 00015 get_atomic_kind 00016 USE cell_types, ONLY: cell_type,& 00017 get_cell 00018 USE cp_files, ONLY: close_file,& 00019 open_file 00020 USE cp_para_types, ONLY: cp_para_env_type 00021 USE cp_subsys_types, ONLY: cp_subsys_get,& 00022 cp_subsys_type 00023 USE cp_units, ONLY: cp_unit_to_cp2k 00024 USE f77_blas 00025 USE fist_environment_types, ONLY: fist_env_get,& 00026 fist_environment_type 00027 USE fist_nonbond_env_types, ONLY: fist_nonbond_env_get,& 00028 fist_nonbond_env_type 00029 USE force_env_types, ONLY: force_env_get,& 00030 force_env_p_type,& 00031 force_env_type,& 00032 use_fist_force,& 00033 use_qs_force 00034 USE global_types, ONLY: global_environment_type 00035 USE input_constants, ONLY: do_fist,& 00036 do_qs 00037 USE input_section_types, ONLY: section_vals_get_subs_vals,& 00038 section_vals_type,& 00039 section_vals_val_get,& 00040 section_vals_val_set 00041 USE kinds, ONLY: default_path_length,& 00042 default_string_length,& 00043 dp,& 00044 dp_size,& 00045 int_size 00046 USE mathconstants, ONLY: pi 00047 USE mol_kind_new_list_types, ONLY: mol_kind_new_list_p_type 00048 USE molecule_kind_types, ONLY: atom_type,& 00049 get_molecule_kind,& 00050 molecule_kind_type 00051 USE pair_potential_types, ONLY: pair_potential_pp_type 00052 USE particle_list_types, ONLY: particle_list_p_type 00053 USE physcon, ONLY: boltzmann,& 00054 joule 00055 USE string_utilities, ONLY: uppercase,& 00056 xstring 00057 USE termination, ONLY: stop_memory,& 00058 stop_program 00059 USE timings, ONLY: timeset,& 00060 timestop 00061 #include "cp_common_uses.h" 00062 00063 IMPLICIT NONE 00064 00065 PRIVATE 00066 00067 ! *** Global parameters *** 00068 00069 PUBLIC :: mc_simpar_type,& 00070 mc_simulation_parameters_p_type,& 00071 mc_averages_type,mc_averages_p_type,& 00072 mc_moves_type, mc_moves_p_type,accattempt,& 00073 get_mc_par,set_mc_par,read_mc_section,& 00074 find_mc_rcut,mc_molecule_info_type,& 00075 mc_determine_molecule_info,mc_molecule_info_destroy,& 00076 mc_sim_par_create,mc_sim_par_destroy,& 00077 get_mc_molecule_info,& 00078 mc_input_file_type,mc_input_file_create,& 00079 get_mc_input_file,mc_input_file_destroy,& 00080 mc_input_parameters_check,& 00081 mc_ekin_type 00082 00083 ! ***************************************************************************** 00084 TYPE mc_simpar_type 00085 ! contains all the parameters needed for running Monte Carlo simulations 00086 PRIVATE 00087 INTEGER,DIMENSION(:),POINTER :: avbmc_atom 00088 INTEGER :: nstep 00089 INTEGER :: iupvolume 00090 INTEGER :: iuptrans 00091 INTEGER :: nvirial 00092 INTEGER :: nbox 00093 INTEGER :: nmoves 00094 INTEGER :: nswapmoves 00095 INTEGER :: rm 00096 INTEGER :: cl 00097 INTEGER :: diff 00098 INTEGER :: nstart 00099 INTEGER :: source 00100 INTEGER :: group 00101 INTEGER :: iprint 00102 LOGICAL :: ldiscrete 00103 LOGICAL :: lbias 00104 LOGICAL :: ionode 00105 LOGICAL :: lrestart 00106 LOGICAL :: lstop 00107 LOGICAL :: lhmc 00108 CHARACTER ( LEN = 20 ) :: ensemble 00109 CHARACTER ( LEN = default_path_length ) :: restart_file_name 00110 CHARACTER ( LEN = default_path_length ) :: molecules_file 00111 CHARACTER ( LEN = default_path_length ) :: moves_file 00112 CHARACTER ( LEN = default_path_length ) :: coords_file 00113 CHARACTER ( LEN = default_path_length ) :: energy_file 00114 CHARACTER ( LEN = default_path_length ) :: displacement_file 00115 CHARACTER ( LEN = default_path_length ) :: cell_file 00116 CHARACTER ( LEN = default_path_length ) :: dat_file 00117 CHARACTER ( LEN = default_path_length ) :: data_file 00118 CHARACTER ( LEN = default_path_length ) :: box2_file 00119 CHARACTER ( LEN = 200 ) :: fft_lib 00120 CHARACTER ( LEN = 50 ) :: PROGRAM 00121 REAL ( dp ),DIMENSION(:),POINTER :: pmtraion_mol,pmtrans_mol,pmrot_mol,pmswap_mol, 00122 pbias,pmavbmc_mol 00123 REAL ( dp ) :: discrete_step 00124 REAL ( dp ) :: rmvolume 00125 REAL ( dp ),DIMENSION(:),POINTER :: rmbond,rmangle,rmdihedral,rmrot,rmtrans 00126 REAL ( dp ),DIMENSION(:),POINTER :: eta 00127 REAL ( dp ) :: temperature 00128 REAL ( dp ) :: pressure 00129 REAL ( dp ) :: pmavbmc 00130 REAL ( dp ) :: pmswap 00131 REAL ( dp ) :: pmvolume,pmvol_box 00132 REAL ( dp ) :: pmhmc,pmhmc_box 00133 REAL ( dp ) :: pmtraion 00134 REAL ( dp ) :: pmtrans 00135 REAL ( dp ) :: BETA 00136 REAL ( dp ) :: rcut 00137 REAL ( dp ),DIMENSION(:),POINTER :: avbmc_rmin,avbmc_rmax 00138 REAL ( dp ),DIMENSION(:),POINTER :: virial_temps 00139 TYPE(mc_input_file_type),POINTER :: 00140 mc_input_file,mc_bias_file 00141 TYPE(section_vals_type),POINTER :: input_file 00142 TYPE(mc_molecule_info_type),POINTER :: mc_molecule_info 00143 REAL ( dp ) :: exp_min_val, exp_max_val, min_val, max_val 00144 INTEGER :: rand2skip 00145 END TYPE mc_simpar_type 00146 00147 ! ***************************************************************************** 00148 TYPE mc_ekin_type 00149 ! contains the kinetic energy of an MD sequence from hybrid MC 00150 REAL ( dp ) :: initial_ekin,final_ekin 00151 END TYPE mc_ekin_type 00152 ! ***************************************************************************** 00153 TYPE mc_input_file_type 00154 ! contains all the text of the input file 00155 PRIVATE 00156 INTEGER :: run_type_row,run_type_column,coord_row_start,coord_row_end, 00157 cell_row,cell_column,global_row_end,in_use,nunits_empty, 00158 motion_row_end,motion_row_start 00159 INTEGER,DIMENSION(:),POINTER :: mol_set_nmol_row,mol_set_nmol_column 00160 CHARACTER ( LEN = default_string_length ),DIMENSION(:),POINTER :: text 00161 CHARACTER ( LEN = default_string_length ),DIMENSION(:),POINTER :: atom_names_empty 00162 REAL(dp),DIMENSION(:,:),POINTER :: coordinates_empty 00163 END TYPE mc_input_file_type 00164 00165 ! ***************************************************************************** 00166 TYPE mc_molecule_info_type 00167 ! contains information on molecules...I had to do this because I use multiple force 00168 ! environments, and they know nothing about each other 00169 PRIVATE 00170 INTEGER :: nboxes 00171 CHARACTER ( LEN = default_string_length ),DIMENSION(:),POINTER :: names 00172 CHARACTER ( LEN = default_string_length ), 00173 DIMENSION(:,:),POINTER :: atom_names 00174 REAL ( dp ),DIMENSION(:,:),POINTER :: conf_prob,mass 00175 INTEGER, DIMENSION(:,:),POINTER :: nchains 00176 INTEGER :: nmol_types,nchain_total 00177 INTEGER,DIMENSION(:),POINTER :: nunits,mol_type,nunits_tot,in_box 00178 END TYPE mc_molecule_info_type 00179 00180 ! ***************************************************************************** 00181 TYPE mc_simulation_parameters_p_type 00182 ! a pointer for multiple boxes 00183 TYPE (mc_simpar_type),POINTER :: mc_par 00184 END TYPE mc_simulation_parameters_p_type 00185 00186 ! ***************************************************************************** 00187 TYPE mc_averages_type 00188 ! contains some data that is averaged throughout the course of a run 00189 REAL(KIND = dp) :: ave_energy 00190 REAL(KIND = dp) :: ave_energy_squared 00191 REAL(KIND = dp) :: ave_volume 00192 REAL(KIND = dp) :: molecules 00193 END TYPE mc_averages_type 00194 00195 ! ***************************************************************************** 00196 TYPE mc_averages_p_type 00197 TYPE (mc_averages_type),POINTER :: averages 00198 END TYPE mc_averages_p_type 00199 00200 ! ***************************************************************************** 00201 TYPE mc_moves_type 00202 ! contains data for how many moves of a give type we've accepted/rejected 00203 TYPE ( accattempt ), POINTER :: bias_bond 00204 TYPE ( accattempt ), POINTER :: bias_angle 00205 TYPE ( accattempt ), POINTER :: bias_dihedral 00206 TYPE ( accattempt ), POINTER :: bias_trans 00207 TYPE ( accattempt ), POINTER :: bias_rot 00208 TYPE ( accattempt ), POINTER :: bond 00209 TYPE ( accattempt ), POINTER :: angle 00210 TYPE ( accattempt ), POINTER :: dihedral 00211 TYPE ( accattempt ), POINTER :: trans 00212 TYPE ( accattempt ), POINTER :: rot 00213 TYPE ( accattempt ), POINTER :: swap 00214 TYPE ( accattempt ), POINTER :: volume 00215 TYPE ( accattempt ), POINTER :: hmc 00216 TYPE ( accattempt ), POINTER :: avbmc_inin 00217 TYPE ( accattempt ), POINTER :: avbmc_outin 00218 TYPE ( accattempt ), POINTER :: avbmc_inout 00219 TYPE ( accattempt ), POINTER :: avbmc_outout 00220 TYPE ( accattempt ), POINTER :: Quickstep 00221 REAL(KIND = dp) :: trans_dis,qtrans_dis 00222 INTEGER :: empty,grown,empty_conf,empty_avbmc 00223 END TYPE mc_moves_type 00224 00225 ! ***************************************************************************** 00226 TYPE accattempt 00227 INTEGER :: successes 00228 INTEGER :: qsuccesses 00229 INTEGER :: attempts 00230 END TYPE accattempt 00231 00232 ! ***************************************************************************** 00233 TYPE mc_moves_p_type 00234 TYPE(mc_moves_type), POINTER :: moves 00235 END TYPE mc_moves_p_type 00236 00237 ! *** Global parameters *** 00238 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mc_types' 00239 00240 CONTAINS 00241 00242 ! ***************************************************************************** 00249 SUBROUTINE get_mc_input_file ( mc_input_file,run_type_row,run_type_column,& 00250 coord_row_start,coord_row_end,cell_row,cell_column,global_row_end,& 00251 mol_set_nmol_row,mol_set_nmol_column,in_use,text,atom_names_empty,& 00252 nunits_empty,coordinates_empty,motion_row_end,motion_row_start) 00253 00254 TYPE(mc_input_file_type), POINTER :: mc_input_file 00255 INTEGER, INTENT(OUT), OPTIONAL :: run_type_row, run_type_column, 00256 coord_row_start, coord_row_end, cell_row, cell_column, global_row_end 00257 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: mol_set_nmol_row, 00258 mol_set_nmol_column 00259 INTEGER, INTENT(OUT), OPTIONAL :: in_use 00260 CHARACTER(LEN=default_string_length), 00261 DIMENSION(:), OPTIONAL, POINTER :: text, atom_names_empty 00262 INTEGER, INTENT(OUT), OPTIONAL :: nunits_empty 00263 REAL(dp), DIMENSION(:, :), OPTIONAL, 00264 POINTER :: coordinates_empty 00265 INTEGER, INTENT(OUT), OPTIONAL :: motion_row_end, 00266 motion_row_start 00267 00268 IF(PRESENT(in_use)) in_use=mc_input_file%in_use 00269 IF(PRESENT(run_type_row)) run_type_row=mc_input_file%run_type_row 00270 IF(PRESENT(run_type_column)) run_type_column=mc_input_file%run_type_column 00271 IF(PRESENT(coord_row_start)) coord_row_start=mc_input_file%coord_row_start 00272 IF(PRESENT(coord_row_end)) coord_row_end=mc_input_file%coord_row_end 00273 IF(PRESENT(cell_row)) cell_row=mc_input_file%cell_row 00274 IF(PRESENT(cell_column)) cell_column=mc_input_file%cell_column 00275 IF(PRESENT(global_row_end)) global_row_end=mc_input_file%global_row_end 00276 IF(PRESENT(nunits_empty)) nunits_empty=mc_input_file%nunits_empty 00277 IF(PRESENT(motion_row_start)) motion_row_start=mc_input_file%motion_row_start 00278 IF(PRESENT(motion_row_end)) motion_row_end=mc_input_file%motion_row_end 00279 00280 IF(PRESENT(mol_set_nmol_row)) mol_set_nmol_row=> mc_input_file%mol_set_nmol_row 00281 IF(PRESENT(mol_set_nmol_column)) mol_set_nmol_column=> mc_input_file%mol_set_nmol_column 00282 IF(PRESENT(text)) text=> mc_input_file%text 00283 IF(PRESENT(atom_names_empty)) atom_names_empty=> mc_input_file%atom_names_empty 00284 IF(PRESENT(coordinates_empty)) coordinates_empty=> mc_input_file%coordinates_empty 00285 00286 END SUBROUTINE GET_MC_INPUT_FILE 00287 ! ***************************************************************************** 00288 SUBROUTINE get_mc_par ( mc_par, nstep, nvirial, iuptrans, iupvolume, & 00289 nmoves,nswapmoves,rm,cl,diff,nstart,& 00290 source,group,lbias,ionode,lrestart,lstop,rmvolume,rmbond,rmangle,& 00291 rmrot,rmtrans,temperature,pressure,BETA,pmswap,pmvolume,pmtraion,pmtrans,& 00292 ensemble,PROGRAM,restart_file_name,molecules_file,moves_file,coords_file,& 00293 energy_file,displacement_file,cell_file,dat_file,data_file,box2_file,& 00294 fft_lib,iprint,rcut,ldiscrete,discrete_step,& 00295 pmavbmc,pbias,avbmc_atom,avbmc_rmin,avbmc_rmax,rmdihedral,& 00296 input_file,mc_molecule_info,pmswap_mol,pmavbmc_mol,pmtrans_mol,pmrot_mol,& 00297 pmtraion_mol,mc_input_file,mc_bias_file,pmvol_box,virial_temps,& 00298 exp_min_val, exp_max_val, min_val, max_val, eta, pmhmc, pmhmc_box,lhmc,rand2skip) 00299 00300 00301 ! ***************************************************************************** 00307 TYPE(mc_simpar_type), POINTER :: mc_par 00308 INTEGER, INTENT(OUT), OPTIONAL :: nstep, nvirial, iuptrans, 00309 iupvolume, nmoves, 00310 nswapmoves, rm, cl, diff, 00311 nstart, source, group 00312 LOGICAL, INTENT(OUT), OPTIONAL :: lbias, ionode, lrestart, lstop 00313 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: rmvolume 00314 REAL(KIND=dp), DIMENSION(:), OPTIONAL, 00315 POINTER :: rmbond, rmangle, rmrot, 00316 rmtrans 00317 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: temperature, pressure, BETA, 00318 pmswap, pmvolume, pmtraion, 00319 pmtrans 00320 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: ensemble, PROGRAM, 00321 restart_file_name, molecules_file, moves_file, coords_file, 00322 energy_file, displacement_file, cell_file, dat_file, data_file, 00323 box2_file, fft_lib 00324 INTEGER, INTENT(OUT), OPTIONAL :: iprint 00325 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: rcut 00326 LOGICAL, INTENT(OUT), OPTIONAL :: ldiscrete 00327 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: discrete_step, pmavbmc 00328 REAL(KIND=dp), DIMENSION(:), OPTIONAL, 00329 POINTER :: pbias 00330 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: avbmc_atom 00331 REAL(KIND=dp), DIMENSION(:), OPTIONAL, 00332 POINTER :: avbmc_rmin, avbmc_rmax, 00333 rmdihedral 00334 TYPE(section_vals_type), OPTIONAL, 00335 POINTER :: input_file 00336 TYPE(mc_molecule_info_type), OPTIONAL, 00337 POINTER :: mc_molecule_info 00338 REAL(dp), DIMENSION(:), OPTIONAL, 00339 POINTER :: pmswap_mol 00340 REAL(KIND=dp), DIMENSION(:), OPTIONAL, 00341 POINTER :: pmavbmc_mol, pmtrans_mol, 00342 pmrot_mol, pmtraion_mol 00343 TYPE(mc_input_file_type), OPTIONAL, 00344 POINTER :: mc_input_file, mc_bias_file 00345 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: pmvol_box 00346 REAL(KIND=dp), DIMENSION(:), OPTIONAL, 00347 POINTER :: virial_temps 00348 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: exp_min_val, exp_max_val, 00349 min_val, max_val 00350 REAL(KIND=dp), DIMENSION(:), OPTIONAL, 00351 POINTER :: eta 00352 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: pmhmc, pmhmc_box 00353 LOGICAL, INTENT(OUT), OPTIONAL :: lhmc 00354 INTEGER, INTENT(OUT), OPTIONAL :: rand2skip 00355 00356 CHARACTER(LEN=*), PARAMETER :: routineN = 'get_mc_par', 00357 routineP = moduleN//':'//routineN 00358 00359 IF ( PRESENT ( nstep ) ) nstep = mc_par % nstep 00360 IF ( PRESENT ( nvirial ) ) nvirial = mc_par % nvirial 00361 IF ( PRESENT ( iuptrans ) ) iuptrans = mc_par % iuptrans 00362 IF ( PRESENT ( iupvolume ) ) iupvolume = mc_par % iupvolume 00363 IF ( PRESENT ( nmoves ) ) nmoves = mc_par % nmoves 00364 IF ( PRESENT ( nswapmoves ) ) nswapmoves = mc_par % nswapmoves 00365 IF ( PRESENT ( rm ) ) rm = mc_par % rm 00366 IF ( PRESENT ( cl ) ) cl = mc_par % cl 00367 IF ( PRESENT ( diff ) ) diff = mc_par % diff 00368 IF ( PRESENT ( nstart ) ) nstart = mc_par % nstart 00369 IF ( PRESENT ( source ) ) source = mc_par % source 00370 IF ( PRESENT ( group ) ) group = mc_par % group 00371 IF ( PRESENT ( iprint ) ) iprint = mc_par % iprint 00372 00373 IF ( PRESENT ( lbias ) ) lbias = mc_par % lbias 00374 IF ( PRESENT ( ionode ) ) ionode = mc_par % ionode 00375 IF ( PRESENT ( lrestart ) ) lrestart = mc_par % lrestart 00376 IF ( PRESENT ( lstop ) ) lstop = mc_par % lstop 00377 IF ( PRESENT ( ldiscrete ) ) ldiscrete = mc_par % ldiscrete 00378 00379 IF ( PRESENT ( virial_temps ) ) virial_temps => mc_par % virial_temps 00380 IF ( PRESENT ( avbmc_atom ) ) avbmc_atom => mc_par % avbmc_atom 00381 IF ( PRESENT ( avbmc_rmin ) ) avbmc_rmin => mc_par % avbmc_rmin 00382 IF ( PRESENT ( avbmc_rmax ) ) avbmc_rmax => mc_par % avbmc_rmax 00383 IF ( PRESENT ( rcut ) ) rcut = mc_par % rcut 00384 IF ( PRESENT ( discrete_step ) ) discrete_step = mc_par % discrete_step 00385 IF ( PRESENT ( rmvolume ) ) rmvolume = mc_par % rmvolume 00386 IF ( PRESENT ( temperature ) ) temperature = mc_par % temperature 00387 IF ( PRESENT ( pressure ) ) pressure = mc_par % pressure 00388 IF ( PRESENT ( BETA ) ) BETA = mc_par % BETA 00389 IF ( PRESENT ( exp_max_val ) ) exp_max_val = mc_par % exp_max_val 00390 IF ( PRESENT ( exp_min_val ) ) exp_min_val = mc_par % exp_min_val 00391 IF ( PRESENT ( max_val ) ) max_val = mc_par % max_val 00392 IF ( PRESENT ( min_val ) ) min_val = mc_par % min_val 00393 IF ( PRESENT ( pmswap ) ) pmswap = mc_par % pmswap 00394 IF ( PRESENT ( pmvolume) ) pmvolume = mc_par % pmvolume 00395 IF ( PRESENT ( lhmc) ) lhmc = mc_par % lhmc 00396 IF ( PRESENT ( pmhmc) ) pmhmc = mc_par % pmhmc 00397 IF ( PRESENT ( pmhmc_box) ) pmhmc_box = mc_par % pmhmc_box 00398 IF ( PRESENT ( pmvol_box) ) pmvol_box = mc_par % pmvol_box 00399 IF ( PRESENT ( pmtraion ) ) pmtraion = mc_par % pmtraion 00400 IF ( PRESENT ( pmtrans ) ) pmtrans = mc_par % pmtrans 00401 IF ( PRESENT ( pmavbmc) ) pmavbmc = mc_par % pmavbmc 00402 IF ( PRESENT ( pmrot_mol) ) pmrot_mol => mc_par % pmrot_mol 00403 IF ( PRESENT ( pmtrans_mol) ) pmtrans_mol => mc_par % pmtrans_mol 00404 IF ( PRESENT ( pmtraion_mol) ) pmtraion_mol => mc_par % pmtraion_mol 00405 IF ( PRESENT ( pmavbmc_mol) ) pmavbmc_mol => mc_par % pmavbmc_mol 00406 IF ( PRESENT ( pbias) ) pbias => mc_par % pbias 00407 00408 IF ( PRESENT ( rmbond ) ) rmbond => mc_par % rmbond 00409 IF ( PRESENT ( rmangle ) ) rmangle => mc_par % rmangle 00410 IF ( PRESENT ( rmdihedral ) ) rmdihedral => mc_par % rmdihedral 00411 IF ( PRESENT ( rmrot ) ) rmrot => mc_par % rmrot 00412 IF ( PRESENT ( rmtrans ) ) rmtrans => mc_par % rmtrans 00413 00414 IF ( PRESENT ( eta ) ) eta => mc_par % eta 00415 00416 IF ( PRESENT ( ensemble ) ) ensemble = mc_par % ensemble 00417 IF ( PRESENT ( PROGRAM ) ) PROGRAM = mc_par % program 00418 IF ( PRESENT ( restart_file_name ) ) restart_file_name = & 00419 mc_par % restart_file_name 00420 IF ( PRESENT ( moves_file ) ) moves_file = mc_par % moves_file 00421 IF ( PRESENT ( coords_file ) ) coords_file = mc_par % coords_file 00422 IF ( PRESENT ( molecules_file ) ) molecules_file = mc_par % molecules_file 00423 IF ( PRESENT ( energy_file ) ) energy_file = mc_par % energy_file 00424 IF ( PRESENT ( displacement_file ) ) displacement_file = & 00425 mc_par % displacement_file 00426 IF ( PRESENT ( cell_file ) ) cell_file = mc_par % cell_file 00427 IF ( PRESENT ( dat_file ) ) dat_file = mc_par % dat_file 00428 IF ( PRESENT ( data_file ) ) data_file = mc_par % data_file 00429 IF ( PRESENT ( box2_file ) ) box2_file = mc_par % box2_file 00430 IF ( PRESENT ( fft_lib ) ) fft_lib = mc_par % fft_lib 00431 00432 IF ( PRESENT ( input_file ) ) input_file => mc_par % input_file 00433 IF ( PRESENT ( mc_molecule_info ) ) mc_molecule_info => mc_par % mc_molecule_info 00434 IF ( PRESENT ( mc_input_file ) ) mc_input_file => mc_par % mc_input_file 00435 IF ( PRESENT ( mc_bias_file ) ) mc_bias_file => mc_par % mc_bias_file 00436 00437 IF ( PRESENT ( pmswap_mol ) ) pmswap_mol => mc_par % pmswap_mol 00438 IF ( PRESENT ( rand2skip ) ) rand2skip = mc_par % rand2skip 00439 END SUBROUTINE get_mc_par 00440 00441 ! ***************************************************************************** 00442 SUBROUTINE get_mc_molecule_info ( mc_molecule_info, nmol_types,nchain_total,nboxes,& 00443 names,conf_prob,nchains,nunits,mol_type,nunits_tot,in_box,atom_names,& 00444 mass) 00445 00446 ! ***************************************************************************** 00452 TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info 00453 INTEGER, INTENT(OUT), OPTIONAL :: nmol_types, nchain_total, 00454 nboxes 00455 CHARACTER(LEN=default_string_length), 00456 DIMENSION(:), OPTIONAL, POINTER :: names 00457 REAL(dp), DIMENSION(:, :), OPTIONAL, 00458 POINTER :: conf_prob 00459 INTEGER, DIMENSION(:, :), OPTIONAL, 00460 POINTER :: nchains 00461 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: nunits, mol_type, nunits_tot, 00462 in_box 00463 CHARACTER(LEN=default_string_length), 00464 DIMENSION(:, :), OPTIONAL, POINTER :: atom_names 00465 REAL(dp), DIMENSION(:, :), OPTIONAL, 00466 POINTER :: mass 00467 00468 CHARACTER(LEN=*), PARAMETER :: routineN = 'get_mc_molecule_info', 00469 routineP = moduleN//':'//routineN 00470 00471 IF ( PRESENT ( nboxes ) ) nboxes = mc_molecule_info % nboxes 00472 IF ( PRESENT ( nchain_total ) ) nchain_total = mc_molecule_info % nchain_total 00473 IF ( PRESENT ( nmol_types ) ) nmol_types = mc_molecule_info % nmol_types 00474 00475 IF ( PRESENT ( names ) ) names => mc_molecule_info % names 00476 IF ( PRESENT ( atom_names ) ) atom_names => mc_molecule_info % atom_names 00477 00478 IF ( PRESENT ( conf_prob ) ) conf_prob => mc_molecule_info % conf_prob 00479 IF ( PRESENT ( mass ) ) mass => mc_molecule_info % mass 00480 00481 IF ( PRESENT ( nchains ) ) nchains => mc_molecule_info % nchains 00482 IF ( PRESENT ( nunits ) ) nunits => mc_molecule_info % nunits 00483 IF ( PRESENT ( mol_type ) ) mol_type => mc_molecule_info % mol_type 00484 IF ( PRESENT ( nunits_tot ) ) nunits_tot => mc_molecule_info % nunits_tot 00485 IF ( PRESENT ( in_box ) ) in_box => mc_molecule_info % in_box 00486 00487 END SUBROUTINE get_mc_molecule_info 00488 00489 ! ***************************************************************************** 00490 SUBROUTINE set_mc_par ( mc_par, rm,cl,& 00491 diff,nstart,rmvolume,rmbond,rmangle,rmdihedral,rmrot,rmtrans,PROGRAM,& 00492 nmoves,nswapmoves,lstop,temperature,pressure,iuptrans,iupvolume,& 00493 pmswap,pmvolume,pmtraion,pmtrans,BETA,rcut,iprint,lbias,nstep,& 00494 lrestart,ldiscrete,discrete_step,pmavbmc,mc_molecule_info,& 00495 pmavbmc_mol,pmtrans_mol,pmrot_mol,pmtraion_mol,pmswap_mol,& 00496 avbmc_rmin,avbmc_rmax,avbmc_atom,pbias,ensemble,pmvol_box,eta,& 00497 mc_input_file,mc_bias_file,exp_max_val,exp_min_val,min_val,max_val,& 00498 pmhmc,pmhmc_box,lhmc,ionode,source,group,rand2skip) 00499 00500 00501 ! ***************************************************************************** 00507 TYPE(mc_simpar_type), POINTER :: mc_par 00508 INTEGER, INTENT(IN), OPTIONAL :: rm, cl, diff, nstart 00509 REAL(KIND=dp), INTENT(IN), OPTIONAL :: rmvolume 00510 REAL(KIND=dp), DIMENSION(:), OPTIONAL, 00511 POINTER :: rmbond, rmangle, rmdihedral, 00512 rmrot, rmtrans 00513 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: PROGRAM 00514 INTEGER, INTENT(IN), OPTIONAL :: nmoves, nswapmoves 00515 LOGICAL, INTENT(IN), OPTIONAL :: lstop 00516 REAL(KIND=dp), INTENT(IN), OPTIONAL :: temperature, pressure 00517 INTEGER, INTENT(IN), OPTIONAL :: iuptrans, iupvolume 00518 REAL(KIND=dp), INTENT(IN), OPTIONAL :: pmswap, pmvolume, pmtraion, 00519 pmtrans, BETA, rcut 00520 INTEGER, INTENT(IN), OPTIONAL :: iprint 00521 LOGICAL, INTENT(IN), OPTIONAL :: lbias 00522 INTEGER, INTENT(IN), OPTIONAL :: nstep 00523 LOGICAL, INTENT(IN), OPTIONAL :: lrestart, ldiscrete 00524 REAL(KIND=dp), INTENT(IN), OPTIONAL :: discrete_step, pmavbmc 00525 TYPE(mc_molecule_info_type), OPTIONAL, 00526 POINTER :: mc_molecule_info 00527 REAL(dp), DIMENSION(:), OPTIONAL, 00528 POINTER :: pmavbmc_mol, pmtrans_mol, 00529 pmrot_mol, pmtraion_mol, 00530 pmswap_mol, avbmc_rmin, 00531 avbmc_rmax 00532 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: avbmc_atom 00533 REAL(dp), DIMENSION(:), OPTIONAL, 00534 POINTER :: pbias 00535 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: ensemble 00536 REAL(KIND=dp), INTENT(IN), OPTIONAL :: pmvol_box 00537 REAL(KIND=dp), DIMENSION(:), OPTIONAL, 00538 POINTER :: eta 00539 TYPE(mc_input_file_type), OPTIONAL, 00540 POINTER :: mc_input_file, mc_bias_file 00541 REAL(KIND=dp), INTENT(IN), OPTIONAL :: exp_max_val, exp_min_val, 00542 min_val, max_val, pmhmc, 00543 pmhmc_box 00544 LOGICAL, INTENT(IN), OPTIONAL :: lhmc, ionode 00545 INTEGER, INTENT(IN), OPTIONAL :: source, group, rand2skip 00546 00547 CHARACTER(LEN=*), PARAMETER :: routineN = 'set_mc_par', 00548 routineP = moduleN//':'//routineN 00549 00550 INTEGER :: iparm 00551 00552 ! These are the only values that change during the course of the simulation 00553 ! or are computed outside of this module 00554 00555 IF ( PRESENT ( nstep ) ) mc_par % nstep = nstep 00556 IF ( PRESENT ( rm ) ) mc_par % rm = rm 00557 IF ( PRESENT ( cl ) ) mc_par % cl = cl 00558 IF ( PRESENT ( diff ) ) mc_par % diff = diff 00559 IF ( PRESENT ( nstart ) ) mc_par % nstart = nstart 00560 IF ( PRESENT ( nmoves ) ) mc_par % nmoves = nmoves 00561 IF ( PRESENT ( nswapmoves ) ) mc_par % nswapmoves = nswapmoves 00562 IF ( PRESENT ( iprint ) ) mc_par % iprint = iprint 00563 IF ( PRESENT ( iuptrans ) ) mc_par % iuptrans = iuptrans 00564 IF ( PRESENT ( iupvolume ) ) mc_par % iupvolume = iupvolume 00565 00566 IF ( PRESENT ( ldiscrete ) ) mc_par % ldiscrete = ldiscrete 00567 IF ( PRESENT ( lstop ) ) mc_par % lstop = lstop 00568 IF ( PRESENT ( lbias ) ) mc_par % lbias = lbias 00569 IF ( PRESENT ( lrestart ) ) mc_par % lrestart = lrestart 00570 00571 IF ( PRESENT ( exp_max_val ) ) mc_par % exp_max_val= exp_max_val 00572 IF ( PRESENT ( exp_min_val ) ) mc_par % exp_min_val= exp_min_val 00573 IF ( PRESENT ( max_val ) ) mc_par % max_val = max_val 00574 IF ( PRESENT ( min_val ) ) mc_par % min_val = min_val 00575 IF ( PRESENT ( BETA ) ) mc_par % BETA = BETA 00576 IF ( PRESENT ( temperature ) ) mc_par % temperature = temperature 00577 IF ( PRESENT ( rcut ) ) mc_par % rcut = rcut 00578 IF ( PRESENT ( pressure ) ) mc_par % pressure = pressure 00579 IF ( PRESENT ( pmvolume ) ) mc_par % pmvolume = pmvolume 00580 IF ( PRESENT ( lhmc ) ) mc_par % lhmc = lhmc 00581 IF ( PRESENT ( pmhmc ) ) mc_par % pmhmc = pmhmc 00582 IF ( PRESENT ( pmhmc_box ) ) mc_par % pmhmc_box = pmhmc_box 00583 IF ( PRESENT ( pmvol_box ) ) mc_par % pmvol_box = pmvol_box 00584 IF ( PRESENT ( pmswap ) ) mc_par % pmswap = pmswap 00585 IF ( PRESENT ( pmtrans ) ) mc_par % pmtrans = pmtrans 00586 IF ( PRESENT ( pmtraion ) ) mc_par % pmtraion = pmtraion 00587 IF ( PRESENT ( pmavbmc ) ) mc_par % pmavbmc = pmavbmc 00588 00589 IF ( PRESENT ( discrete_step ) ) mc_par % discrete_step = discrete_step 00590 IF ( PRESENT ( rmvolume ) ) mc_par % rmvolume = rmvolume 00591 00592 IF ( PRESENT ( mc_input_file ) ) mc_par % mc_input_file => mc_input_file 00593 IF ( PRESENT ( mc_bias_file ) ) mc_par % mc_bias_file => mc_bias_file 00594 00595 ! cannot just change the pointers, because then we have memory leaks 00596 ! and the program crashes if we try to release it (in certain cases) 00597 IF ( PRESENT ( eta ) ) THEN 00598 DO iparm=1,SIZE(eta) 00599 mc_par % eta(iparm) = eta(iparm) 00600 ENDDO 00601 ENDIF 00602 IF ( PRESENT ( rmbond ) ) THEN 00603 DO iparm=1,SIZE(rmbond) 00604 mc_par % rmbond(iparm) = rmbond(iparm) 00605 ENDDO 00606 ENDIF 00607 IF ( PRESENT ( rmangle ) ) THEN 00608 DO iparm=1,SIZE(rmangle) 00609 mc_par % rmangle(iparm) = rmangle(iparm) 00610 ENDDO 00611 ENDIF 00612 IF ( PRESENT ( rmdihedral ) ) THEN 00613 DO iparm=1,SIZE(rmdihedral) 00614 mc_par % rmdihedral(iparm) = rmdihedral(iparm) 00615 ENDDO 00616 ENDIF 00617 IF ( PRESENT ( rmrot ) ) THEN 00618 DO iparm=1,SIZE(rmrot) 00619 mc_par % rmrot(iparm) = rmrot(iparm) 00620 ENDDO 00621 ENDIF 00622 IF ( PRESENT ( rmtrans ) ) THEN 00623 DO iparm=1,SIZE(rmtrans) 00624 mc_par % rmtrans(iparm) = rmtrans(iparm) 00625 ENDDO 00626 ENDIF 00627 00628 IF ( PRESENT ( pmavbmc_mol ) ) THEN 00629 DO iparm=1,SIZE(pmavbmc_mol) 00630 mc_par % pmavbmc_mol(iparm) = pmavbmc_mol(iparm) 00631 ENDDO 00632 ENDIF 00633 IF ( PRESENT ( pmtrans_mol ) ) THEN 00634 DO iparm=1,SIZE(pmtrans_mol) 00635 mc_par % pmtrans_mol(iparm) = pmtrans_mol(iparm) 00636 ENDDO 00637 ENDIF 00638 IF ( PRESENT ( pmrot_mol ) ) THEN 00639 DO iparm=1,SIZE(pmrot_mol) 00640 mc_par % pmrot_mol(iparm) = pmrot_mol(iparm) 00641 ENDDO 00642 ENDIF 00643 IF ( PRESENT ( pmtraion_mol ) ) THEN 00644 DO iparm=1,SIZE(pmtraion_mol) 00645 mc_par % pmtraion_mol(iparm) = pmtraion_mol(iparm) 00646 ENDDO 00647 ENDIF 00648 IF ( PRESENT ( pmswap_mol ) ) THEN 00649 DO iparm=1,SIZE(pmswap_mol) 00650 mc_par % pmswap_mol(iparm) = pmswap_mol(iparm) 00651 ENDDO 00652 ENDIF 00653 00654 IF ( PRESENT ( avbmc_atom ) ) THEN 00655 DO iparm=1,SIZE(avbmc_atom) 00656 mc_par % avbmc_atom(iparm) = avbmc_atom(iparm) 00657 ENDDO 00658 ENDIF 00659 IF ( PRESENT ( avbmc_rmin ) ) THEN 00660 DO iparm=1,SIZE(avbmc_rmin) 00661 mc_par % avbmc_rmin(iparm) = avbmc_rmin(iparm) 00662 ENDDO 00663 ENDIF 00664 IF ( PRESENT ( avbmc_rmax ) ) THEN 00665 DO iparm=1,SIZE(avbmc_rmax) 00666 mc_par % avbmc_rmax(iparm) = avbmc_rmax(iparm) 00667 ENDDO 00668 ENDIF 00669 IF ( PRESENT ( pbias ) ) THEN 00670 DO iparm=1,SIZE(pbias) 00671 mc_par % pbias(iparm) = pbias(iparm) 00672 ENDDO 00673 ENDIF 00674 00675 IF ( PRESENT ( PROGRAM ) ) mc_par % program = PROGRAM 00676 IF ( PRESENT ( ensemble ) ) mc_par % ensemble = ensemble 00677 00678 IF(PRESENT(mc_molecule_info)) mc_par%mc_molecule_info => mc_molecule_info 00679 IF(PRESENT(source)) mc_par%source=source 00680 IF(PRESENT(group)) mc_par%group=group 00681 IF(PRESENT(ionode)) mc_par%ionode=ionode 00682 00683 IF(PRESENT(rand2skip)) mc_par%rand2skip=rand2skip 00684 END SUBROUTINE set_mc_par 00685 00686 ! ***************************************************************************** 00692 SUBROUTINE mc_sim_par_create ( mc_par ,nmol_types) 00693 00694 TYPE(mc_simpar_type), POINTER :: mc_par 00695 INTEGER, INTENT(IN) :: nmol_types 00696 00697 CHARACTER(len=*), PARAMETER :: routineN = 'mc_sim_par_create', 00698 routineP = moduleN//':'//routineN 00699 00700 INTEGER :: istat 00701 00702 NULLIFY(mc_par) 00703 00704 ALLOCATE(mc_par) 00705 ALLOCATE (mc_par%pmtraion_mol(1:nmol_types),STAT=istat) 00706 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00707 "mc_par%pmtraion_mol",nmol_types*dp_size) 00708 ALLOCATE (mc_par%pmtrans_mol(1:nmol_types),STAT=istat) 00709 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00710 "mc_par%pmtrans_mol",nmol_types*dp_size) 00711 ALLOCATE (mc_par%pmrot_mol(1:nmol_types),STAT=istat) 00712 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00713 "mc_par%pmrot_mol",nmol_types*dp_size) 00714 ALLOCATE (mc_par%pmswap_mol(1:nmol_types),STAT=istat) 00715 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00716 "mc_par%pmswap_mol",nmol_types*dp_size) 00717 00718 ALLOCATE (mc_par%eta(1:nmol_types),STAT=istat) 00719 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00720 "mc_par%eta",nmol_types*dp_size) 00721 00722 ALLOCATE (mc_par%rmbond(1:nmol_types),STAT=istat) 00723 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00724 "mc_par%rmbond",nmol_types*dp_size) 00725 ALLOCATE (mc_par%rmangle(1:nmol_types),STAT=istat) 00726 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00727 "mc_par%rmangle",nmol_types*dp_size) 00728 ALLOCATE (mc_par%rmdihedral(1:nmol_types),STAT=istat) 00729 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00730 "mc_par%rmdihedral",nmol_types*dp_size) 00731 ALLOCATE (mc_par%rmtrans(1:nmol_types),STAT=istat) 00732 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00733 "mc_par%rmtrans",nmol_types*dp_size) 00734 ALLOCATE (mc_par%rmrot(1:nmol_types),STAT=istat) 00735 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00736 "mc_par%rmrot",nmol_types*dp_size) 00737 00738 ALLOCATE (mc_par%avbmc_atom(1:nmol_types),STAT=istat) 00739 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00740 "mc_par%avbmc_atom",nmol_types*int_size) 00741 ALLOCATE (mc_par%avbmc_rmin(1:nmol_types),STAT=istat) 00742 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00743 "mc_par%avbmc_rmin",nmol_types*dp_size) 00744 ALLOCATE (mc_par%avbmc_rmax(1:nmol_types),STAT=istat) 00745 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00746 "mc_par%avbmc_rmax",nmol_types*dp_size) 00747 ALLOCATE (mc_par%pmavbmc_mol(1:nmol_types),STAT=istat) 00748 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00749 "mc_par%pmavbmc_mol",nmol_types*dp_size) 00750 ALLOCATE (mc_par%pbias(1:nmol_types),STAT=istat) 00751 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00752 "mc_par%pbias",nmol_types*dp_size) 00753 00754 ALLOCATE (mc_par%mc_input_file,STAT=istat) 00755 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00756 "mc_par%mc_input_file",1) 00757 ALLOCATE (mc_par%mc_bias_file,STAT=istat) 00758 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00759 "mc_par%mc_bias_file",1) 00760 00761 END SUBROUTINE mc_sim_par_create 00762 00763 ! ***************************************************************************** 00768 SUBROUTINE mc_sim_par_destroy ( mc_par ) 00769 00770 TYPE(mc_simpar_type), POINTER :: mc_par 00771 00772 CHARACTER(len=*), PARAMETER :: routineN = 'mc_sim_par_destroy', 00773 routineP = moduleN//':'//routineN 00774 00775 INTEGER :: istat 00776 00777 DEALLOCATE (mc_par%mc_input_file,STAT=istat) 00778 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00779 "mc_par%mc_input_file") 00780 DEALLOCATE (mc_par%mc_bias_file,STAT=istat) 00781 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00782 "mc_par%mc_bias_file") 00783 00784 DEALLOCATE (mc_par%pmtraion_mol,STAT=istat) 00785 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00786 "mc_par%pmtraion_mol") 00787 DEALLOCATE (mc_par%pmtrans_mol,STAT=istat) 00788 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00789 "mc_par%pmtrans_mol") 00790 DEALLOCATE (mc_par%pmrot_mol,STAT=istat) 00791 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00792 "mc_par%pmrot_mol") 00793 DEALLOCATE (mc_par%pmswap_mol,STAT=istat) 00794 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00795 "mc_par%pmswap_mol") 00796 00797 DEALLOCATE (mc_par%eta,STAT=istat) 00798 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00799 "mc_par%eta") 00800 00801 DEALLOCATE (mc_par%rmbond,STAT=istat) 00802 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00803 "mc_par%rmbond") 00804 DEALLOCATE (mc_par%rmangle,STAT=istat) 00805 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00806 "mc_par%rmangle") 00807 DEALLOCATE (mc_par%rmdihedral,STAT=istat) 00808 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00809 "mc_par%rmdihedral") 00810 DEALLOCATE (mc_par%rmtrans,STAT=istat) 00811 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00812 "mc_par%rmtrans") 00813 DEALLOCATE (mc_par%rmrot,STAT=istat) 00814 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00815 "mc_par%rmrot") 00816 00817 DEALLOCATE (mc_par%avbmc_atom,STAT=istat) 00818 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00819 "mc_par%avbmc_atom") 00820 DEALLOCATE (mc_par%avbmc_rmin,STAT=istat) 00821 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00822 "mc_par%avbmc_rmin") 00823 DEALLOCATE (mc_par%avbmc_rmax,STAT=istat) 00824 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00825 "mc_par%avbmc_rmax") 00826 DEALLOCATE (mc_par%pmavbmc_mol,STAT=istat) 00827 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00828 "mc_par%pmavbmc_mol") 00829 DEALLOCATE (mc_par%pbias,STAT=istat) 00830 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00831 "mc_par%pbias") 00832 IF(mc_par%ensemble == "VIRIAL") THEN 00833 DEALLOCATE (mc_par%virial_temps,STAT=istat) 00834 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00835 "mc_par%virial_temps") 00836 END IF 00837 00838 DEALLOCATE (mc_par,STAT=istat) 00839 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00840 "mc_par") 00841 NULLIFY(mc_par) 00842 00843 END SUBROUTINE mc_sim_par_destroy 00844 00845 ! ***************************************************************************** 00852 SUBROUTINE mc_input_file_create (mc_input_file,input_file_name,& 00853 mc_molecule_info,empty_coords,lhmc,error) 00854 00855 TYPE(mc_input_file_type), POINTER :: mc_input_file 00856 CHARACTER(LEN=*), INTENT(IN) :: input_file_name 00857 TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info 00858 REAL(dp), DIMENSION(:, :) :: empty_coords 00859 LOGICAL, INTENT(IN) :: lhmc 00860 TYPE(cp_error_type), INTENT(inout) :: error 00861 00862 CHARACTER(len=*), PARAMETER :: routineN = 'mc_input_file_create', 00863 routineP = moduleN//':'//routineN 00864 00865 CHARACTER(default_string_length) :: cdum, line, method_name 00866 CHARACTER(default_string_length), 00867 DIMENSION(:, :), POINTER :: atom_names 00868 INTEGER :: abc_column, abc_row, cell_column, cell_row, idum, iline, 00869 io_stat, irep, istat, iunit, iw, line_number, nlines, nmol_types, 00870 nstart, row_number, unit 00871 INTEGER, DIMENSION(:), POINTER :: nunits 00872 TYPE(cp_logger_type), POINTER :: logger 00873 00874 ! some stuff in case we need to write error messages to the screen 00875 00876 logger=>cp_error_get_logger(error) 00877 iw = cp_logger_get_default_io_unit(logger) 00878 00879 ! allocate the array we'll need in case we have an empty box 00880 CALL get_mc_molecule_info(mc_molecule_info,nmol_types=nmol_types,& 00881 nunits=nunits,atom_names=atom_names) 00882 ALLOCATE (mc_input_file%coordinates_empty(1:3,1:nunits(1)),STAT=istat) 00883 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00884 "mc_input_file%coordinates_empty",nunits(1)*5) 00885 ALLOCATE (mc_input_file%atom_names_empty(1:nunits(1)),STAT=istat) 00886 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00887 "mc_input_file%atom_names_empty",nunits(1)*5) 00888 DO iunit=1,nunits(1) 00889 mc_input_file%atom_names_empty(iunit)=atom_names(iunit,1) 00890 mc_input_file%coordinates_empty(1:3,iunit)=empty_coords(1:3,iunit) 00891 ENDDO 00892 mc_input_file%nunits_empty=nunits(1) 00893 00894 ! allocate some arrays 00895 ALLOCATE (mc_input_file%mol_set_nmol_row(1:nmol_types),STAT=istat) 00896 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00897 "mc_input_file%mol_set_nmol_row",nmol_types*int_size) 00898 ALLOCATE (mc_input_file%mol_set_nmol_column(1:nmol_types),STAT=istat) 00899 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00900 "mc_input_file%mol_set_nmol_column",nmol_types*int_size) 00901 00902 ! now read in and store the input file, first finding out how many lines it is 00903 CALL open_file(file_name=input_file_name,& 00904 unit_number=unit,file_action='READ',file_status='OLD') 00905 00906 nlines=0 00907 DO 00908 READ(unit,'(A)',IOSTAT=io_stat) line 00909 IF(io_stat .NE. 0) EXIT 00910 nlines=nlines+1 00911 ENDDO 00912 00913 ALLOCATE (mc_input_file%text(1:nlines),STAT=istat) 00914 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00915 "mc_input_file%text",nlines*default_string_length) 00916 00917 REWIND(unit) 00918 DO iline=1,nlines 00919 READ(unit,'(A)') mc_input_file%text(iline) 00920 ENDDO 00921 00922 ! close the input file 00923 CALL close_file(unit_number=unit) 00924 00925 ! now we need to find the row and column numbers of a variety of things 00926 ! first, Let's find the first END after GLOBAL 00927 CALL mc_parse_text(mc_input_file%text,1,nlines,"&GLOBAL",.TRUE.,& 00928 mc_input_file%global_row_end,idum) 00929 IF(mc_input_file%global_row_end == 0) THEN 00930 IF (iw>0) THEN 00931 WRITE(iw,*) 00932 WRITE(iw,*) 'File name ',input_file_name 00933 END IF 00934 CALL stop_program(routineN,moduleN,__LINE__,& 00935 'Could not find &END after &GLOBAL (make sure & is in the same column) and last line is blank') 00936 ENDIF 00937 00938 ! Let's find the first END after MOTION...we might need this whole section 00939 ! because of hybrid MD/MC moves, which requires the MD information to always 00940 ! be attached to the force env 00941 CALL mc_parse_text(mc_input_file%text,1,nlines,"&MOTION",.TRUE.,& 00942 mc_input_file%motion_row_end,idum,start_row_number=mc_input_file%motion_row_start) 00943 IF(mc_input_file%motion_row_end == 0) THEN 00944 IF (iw>0) THEN 00945 WRITE(iw,*) 00946 WRITE(iw,*) 'File name ',input_file_name,mc_input_file%motion_row_start 00947 END IF 00948 CALL stop_program(routineN,moduleN,__LINE__,& 00949 'Could not find &END after &MOTION (make sure & is in the same column) and last line is blank') 00950 ENDIF 00951 00952 ! first, let's find the first END after &COORD, and the line of &COORD 00953 CALL mc_parse_text(mc_input_file%text,1,nlines,"&COORD",.TRUE.,& 00954 mc_input_file%coord_row_end,idum) 00955 CALL mc_parse_text(mc_input_file%text,1,nlines,"&COORD",.FALSE.,& 00956 mc_input_file%coord_row_start,idum) 00957 00958 IF(mc_input_file%coord_row_end == 0) THEN 00959 IF (iw>0) THEN 00960 WRITE(iw,*) 00961 WRITE(iw,*) 'File name ',input_file_name 00962 END IF 00963 CALL stop_program(routineN,moduleN,__LINE__,& 00964 'Could not find &END after &COORD (make sure & is the first in the same column after &COORD)') 00965 ENDIF 00966 00967 ! now the RUN_TYPE 00968 CALL mc_parse_text(mc_input_file%text,1,nlines,"RUN_TYPE",.FALSE.,& 00969 mc_input_file%run_type_row, mc_input_file%run_type_column) 00970 mc_input_file%run_type_column=mc_input_file%run_type_column+9 00971 IF(mc_input_file%run_type_row == 0) THEN 00972 IF (iw>0) THEN 00973 WRITE(iw,*) 00974 WRITE(iw,*) 'File name ',input_file_name 00975 END IF 00976 CALL stop_program(routineN,moduleN,__LINE__,& 00977 'Could not find RUN_TYPE in the input file (should be in the &GLOBAL section).') 00978 ENDIF 00979 00980 ! now the CELL stuff..we don't care about REF_CELL since that should 00981 ! never change if it's there 00982 CALL mc_parse_text(mc_input_file%text,1,nlines,"&CELL",.FALSE.,& 00983 mc_input_file%cell_row, mc_input_file%cell_column) 00984 ! now find the ABC input line after CELL 00985 CALL mc_parse_text(mc_input_file%text,mc_input_file%cell_row+1,nlines,& 00986 "ABC",.FALSE.,abc_row,abc_column) 00987 ! is there a &CELL inbetween? If so, that ABC will be for the ref_cell 00988 ! and we need to find the next one 00989 CALL mc_parse_text(mc_input_file%text,mc_input_file%cell_row+1,abc_row,& 00990 "&CELL",.FALSE.,cell_row,cell_column) 00991 IF(cell_row == 0) THEN ! nothing in between...we found the correct ABC 00992 mc_input_file%cell_row=abc_row 00993 mc_input_file%cell_column=abc_column+4 00994 ELSE 00995 CALL mc_parse_text(mc_input_file%text,abc_row+1,nlines,& 00996 "ABC",.FALSE.,mc_input_file%cell_row, mc_input_file%cell_column) 00997 ENDIF 00998 IF(mc_input_file%cell_row == 0) THEN 00999 IF (iw>0) THEN 01000 WRITE(iw,*) 01001 WRITE(iw,*) 'File name ',input_file_name 01002 END IF 01003 CALL stop_program(routineN,moduleN,__LINE__,& 01004 'Could not find &CELL section in the input file. Or could not find the ABC line after it.') 01005 ENDIF 01006 01007 ! now we need all the MOL_SET NMOLS indicies 01008 IF(.NOT. lhmc)THEN 01009 irep=0 01010 nstart=1 01011 DO 01012 CALL mc_parse_text(mc_input_file%text,nstart,nlines,"&MOLECULE",& 01013 .FALSE.,row_number, idum) 01014 IF(row_number == 0) EXIT 01015 nstart=row_number+1 01016 irep=irep+1 01017 CALL mc_parse_text(mc_input_file%text,nstart,nlines,"NMOL",& 01018 .FALSE.,mc_input_file%mol_set_nmol_row(irep), & 01019 mc_input_file%mol_set_nmol_column(irep)) 01020 mc_input_file%mol_set_nmol_column(irep)=mc_input_file%mol_set_nmol_column(irep)+5 01021 01022 ENDDO 01023 IF(irep .NE. nmol_types) THEN 01024 IF (iw>0) THEN 01025 WRITE(iw,*) 01026 WRITE(iw,*) 'File name ',input_file_name 01027 WRITE(iw,*) 'Number of molecule types ',nmol_types 01028 WRITE(iw,*) '&MOLECULE sections found ',irep 01029 END IF 01030 CALL stop_program(routineN,moduleN,__LINE__,& 01031 'Did not find MOLECULE sections for every molecule in the simulation (make sure both input files have all types)') 01032 ENDIF 01033 ENDIF 01034 01035 ! now let's find the type of force_env...right now, I'm only concerned with 01036 ! QS, and Fist, though it should be trivial to add others 01037 CALL mc_parse_text(mc_input_file%text,1,nlines,& 01038 "METHOD",.FALSE.,line_number, idum) 01039 READ(mc_input_file%text(line_number),*) cdum,method_name 01040 CALL uppercase(method_name) 01041 SELECT CASE(TRIM(ADJUSTL(method_name))) 01042 CASE("FIST") 01043 mc_input_file%in_use=use_fist_force 01044 CASE("QS","QUICKSTEP") 01045 mc_input_file%in_use=use_qs_force 01046 CASE default 01047 CALL stop_program(routineN,moduleN,__LINE__,& 01048 'Cannot determine the type of force_env we are using (check METHOD)') 01049 END SELECT 01050 01051 END SUBROUTINE mc_input_file_create 01052 01053 ! ***************************************************************************** 01058 SUBROUTINE mc_input_file_destroy(mc_input_file) 01059 01060 TYPE(mc_input_file_type), POINTER :: mc_input_file 01061 01062 CHARACTER(len=*), PARAMETER :: routineN = 'mc_input_file_destroy', 01063 routineP = moduleN//':'//routineN 01064 01065 INTEGER :: istat 01066 01067 DEALLOCATE (mc_input_file%mol_set_nmol_row,STAT=istat) 01068 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01069 "mc_input_file%mol_set_nmol_row") 01070 DEALLOCATE (mc_input_file%mol_set_nmol_column,STAT=istat) 01071 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01072 "mc_input_file%mol_set_nmol_column") 01073 DEALLOCATE (mc_input_file%text,STAT=istat) 01074 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01075 "mc_input_file%text") 01076 DEALLOCATE (mc_input_file%atom_names_empty,STAT=istat) 01077 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01078 "mc_input_file%atom_names_empty") 01079 DEALLOCATE (mc_input_file%coordinates_empty,STAT=istat) 01080 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01081 "mc_input_file%coordinates_empty") 01082 01083 END SUBROUTINE mc_input_file_destroy 01084 01085 ! ***************************************************************************** 01099 SUBROUTINE mc_parse_text(text,nstart,nend,string_search,lend,& 01100 row_number,column_number,start_row_number) 01101 01102 CHARACTER(LEN=*), DIMENSION(:), 01103 INTENT(IN) :: text 01104 INTEGER, INTENT(IN) :: nstart, nend 01105 CHARACTER(LEN=*), INTENT(IN) :: string_search 01106 LOGICAL, INTENT(IN) :: lend 01107 INTEGER, INTENT(OUT) :: row_number, column_number 01108 INTEGER, INTENT(OUT), OPTIONAL :: start_row_number 01109 01110 CHARACTER(default_string_length) :: text_temp 01111 INTEGER :: iline, index_string, jline, 01112 string_size 01113 01114 ! did not see any string utilities in the code to help, so here I go 01115 01116 row_number=0 01117 column_number=0 01118 01119 ! how long is our string to search for? 01120 string_size=LEN_TRIM(string_search) 01121 whole_file:DO iline=nstart,nend 01122 01123 index_string=-1 01124 index_string=INDEX(text(iline),string_search(1:string_size)) 01125 01126 IF(index_string .GT. 0) THEN ! we found one 01127 IF(.NOT. lend) THEN 01128 row_number=iline 01129 column_number=index_string 01130 ELSE 01131 IF(PRESENT(start_row_number)) start_row_number=iline 01132 column_number=index_string 01133 DO jline=iline+1,nend 01134 ! now we find the &END that matches up with this one... 01135 ! I need proper indentation because I'm not very smart 01136 text_temp=text(jline) 01137 IF(text_temp(column_number:column_number) .EQ. "&") THEN 01138 row_number=jline 01139 EXIT 01140 ENDIF 01141 ENDDO 01142 ENDIF 01143 EXIT whole_file 01144 ENDIF 01145 ENDDO whole_file 01146 01147 END SUBROUTINE mc_parse_text 01148 01149 ! ***************************************************************************** 01159 SUBROUTINE read_mc_section ( mc_par, para_env, globenv, input_file_name, input_file, force_env_section,& 01160 error ) 01161 01162 TYPE(mc_simpar_type), POINTER :: mc_par 01163 TYPE(cp_para_env_type), POINTER :: para_env 01164 TYPE(global_environment_type), POINTER :: globenv 01165 CHARACTER(LEN=*), INTENT(IN) :: input_file_name 01166 TYPE(section_vals_type), POINTER :: input_file, force_env_section 01167 TYPE(cp_error_type), INTENT(inout) :: error 01168 01169 CHARACTER(len=*), PARAMETER :: routineN = 'read_mc_section', 01170 routineP = moduleN//':'//routineN 01171 01172 CHARACTER(LEN=5) :: box_string, molecule_string, 01173 tab_box_string, tab_string, 01174 tab_string_int 01175 CHARACTER(LEN=default_string_length) :: format_box_string, 01176 format_string, 01177 format_string_int 01178 INTEGER :: handle, ia, ie, imol, istat, 01179 itype, iw, method_name_id, 01180 nmol_types, stop_num 01181 INTEGER, DIMENSION(:), POINTER :: temp_i_array 01182 LOGICAL :: failure 01183 REAL(dp), DIMENSION(:), POINTER :: temp_r_array 01184 TYPE(cp_logger_type), POINTER :: logger 01185 TYPE(section_vals_type), POINTER :: mc_section 01186 01187 ! begin the timing of the subroutine 01188 01189 CPPrecondition(ASSOCIATED(input_file),cp_failure_level,routineP,error,failure) 01190 01191 CALL timeset(routineN,handle) 01192 01193 NULLIFY(mc_section) 01194 mc_section => section_vals_get_subs_vals(input_file,& 01195 "MOTION%MC",error=error) 01196 01197 ! need the input file sturcutre that we're reading from for when we make 01198 ! dat files 01199 mc_par % input_file => input_file 01200 01201 ! set the ionode and mepos 01202 mc_par % ionode = para_env%ionode 01203 mc_par % group = para_env%group 01204 mc_par % source = para_env%source 01205 01206 ! for convenience 01207 nmol_types=mc_par%mc_molecule_info%nmol_types 01208 01209 !..defaults...most are set in input_cp2k_motion 01210 mc_par % nstart = 0 01211 CALL section_vals_val_get(force_env_section,"METHOD",i_val=method_name_id,error=error) 01212 SELECT CASE (method_name_id) 01213 CASE (do_fist) 01214 mc_par % iprint = 100 01215 CASE (do_qs) 01216 mc_par % iprint = 1 01217 END SELECT 01218 01219 logger=>cp_error_get_logger(error) 01220 iw = cp_logger_get_default_io_unit(logger) 01221 IF(iw>0) WRITE ( iw, * ) 01222 01223 !..filenames 01224 mc_par % program = input_file_name 01225 CALL xstring ( mc_par % program, ia, ie ) 01226 mc_par%coords_file = mc_par % program(ia:ie) // '.coordinates' 01227 mc_par%molecules_file = mc_par % program(ia:ie) // '.molecules' 01228 mc_par%moves_file = mc_par % program(ia:ie) // '.moves' 01229 mc_par%energy_file = mc_par % program(ia:ie) // '.energy' 01230 mc_par%cell_file = mc_par % program(ia:ie) // '.cell' 01231 mc_par%displacement_file= mc_par % program(ia:ie)& 01232 // '.max_displacements' 01233 mc_par%data_file = mc_par % program(ia:ie) // '.data' 01234 stop_num=ie-3 01235 mc_par%dat_file = mc_par % program(ia:stop_num) // 'dat' 01236 01237 ! set them into the input parameter structure as the new defaults 01238 CALL section_vals_val_set(mc_section,"COORDINATE_FILE_NAME",& 01239 c_val=mc_par%coords_file,error=error) 01240 CALL section_vals_val_set(mc_section,"DATA_FILE_NAME",& 01241 c_val=mc_par%data_file,error=error) 01242 CALL section_vals_val_set(mc_section,"CELL_FILE_NAME",& 01243 c_val=mc_par%cell_file,error=error) 01244 CALL section_vals_val_set(mc_section,"MAX_DISP_FILE_NAME",& 01245 c_val=mc_par%displacement_file,error=error) 01246 CALL section_vals_val_set(mc_section,"MOVES_FILE_NAME",& 01247 c_val=mc_par%moves_file,error=error) 01248 CALL section_vals_val_set(mc_section,"MOLECULES_FILE_NAME",& 01249 c_val=mc_par%molecules_file,error=error) 01250 CALL section_vals_val_set(mc_section,"ENERGY_FILE_NAME",& 01251 c_val=mc_par%energy_file,error=error) 01252 01253 ! grab the FFT library name and print level...this is used for writing the dat file 01254 ! and hopefully will be changed 01255 mc_par % fft_lib = globenv%default_fft_library 01256 01257 ! first, grab all the integer values 01258 CALL section_vals_val_get(mc_section,"NSTEP",i_val=mc_par%nstep,error=error) 01259 CALL section_vals_val_get(mc_section,"NMOVES",i_val=mc_par%nmoves,error=error) 01260 CALL section_vals_val_get(mc_section,"NSWAPMOVES",i_val=mc_par%nswapmoves,error=error) 01261 CALL section_vals_val_get(mc_section,"IUPVOLUME",i_val=mc_par%iupvolume,error=error) 01262 CALL section_vals_val_get(mc_section,"NVIRIAL",i_val=mc_par%nvirial,error=error) 01263 CALL section_vals_val_get(mc_section,"IUPTRANS",i_val=mc_par%iuptrans,error=error) 01264 CALL section_vals_val_get(mc_section,"IPRINT",i_val=mc_par%iprint,error=error) 01265 ! now an integer array 01266 CALL section_vals_val_get(mc_section,"AVBMC_ATOM",i_vals=temp_i_array,error=error) 01267 01268 ! find out if we're only doing HMC moves...we can ignore a lot of extra information 01269 ! then, which would ordinarily be cause for concern 01270 CALL section_vals_val_get(mc_section,"PMHMC",r_val=mc_par%pmhmc,error=error) 01271 CALL section_vals_val_get(mc_section,"PMSWAP",r_val=mc_par%pmswap,error=error) 01272 IF(mc_par%pmhmc - mc_par%pmswap >= 1.0_dp .AND. mc_par%pmhmc == 1.0_dp)THEN 01273 mc_par%lhmc=.TRUE. 01274 ELSE 01275 mc_par%lhmc=.FALSE. 01276 ENDIF 01277 01278 IF( .NOT. mc_par%lhmc)THEN 01279 IF(SIZE(temp_i_array) .NE. nmol_types)THEN 01280 CALL stop_program(routineN,moduleN,__LINE__,& 01281 'AVBMC_ATOM must have as many values as there are molecule types.') 01282 ENDIF 01283 ENDIF 01284 DO imol=1,SIZE(temp_i_array) 01285 mc_par%avbmc_atom(imol)=temp_i_array(imol) 01286 ENDDO 01287 01288 ! now the real values 01289 CALL section_vals_val_get(mc_section,"PRESSURE",r_val=mc_par%pressure,error=error) 01290 CALL section_vals_val_get(mc_section,"TEMPERATURE",r_val=mc_par%temperature,error=error) 01291 CALL section_vals_val_get(mc_section,"PMVOLUME",r_val=mc_par%pmvolume,error=error) 01292 CALL section_vals_val_get(mc_section,"PMHMC_BOX",r_val=mc_par%pmhmc_box,error=error) 01293 CALL section_vals_val_get(mc_section,"PMVOL_BOX",r_val=mc_par%pmvol_box,error=error) 01294 CALL section_vals_val_get(mc_section,"PMAVBMC",r_val=mc_par%pmavbmc,error=error) 01295 CALL section_vals_val_get(mc_section,"PMTRANS",r_val=mc_par%pmtrans,error=error) 01296 CALL section_vals_val_get(mc_section,"PMTRAION",r_val=mc_par%pmtraion,error=error) 01297 CALL section_vals_val_get(mc_section,"DISCRETE_STEP",r_val=mc_par%discrete_step,error=error) 01298 CALL section_vals_val_get(mc_section,"RMVOLUME",r_val=mc_par%rmvolume,error=error) 01299 01300 ! finally the character values 01301 CALL section_vals_val_get(mc_section,"ENSEMBLE",c_val=mc_par%ensemble,error=error) 01302 CALL section_vals_val_get(mc_section,"RESTART_FILE_NAME",c_val=mc_par%restart_file_name,error=error) 01303 CALL section_vals_val_get(mc_section,"COORDINATE_FILE_NAME",c_val=mc_par%coords_file,error=error) 01304 CALL section_vals_val_get(mc_section,"ENERGY_FILE_NAME",c_val=mc_par%energy_file,error=error) 01305 CALL section_vals_val_get(mc_section,"MOVES_FILE_NAME",c_val=mc_par%moves_file,error=error) 01306 CALL section_vals_val_get(mc_section,"MOLECULES_FILE_NAME",c_val=mc_par%molecules_file,error=error) 01307 CALL section_vals_val_get(mc_section,"CELL_FILE_NAME",c_val=mc_par%cell_file,error=error) 01308 CALL section_vals_val_get(mc_section,"DATA_FILE_NAME",c_val=mc_par%data_file,error=error) 01309 CALL section_vals_val_get(mc_section,"MAX_DISP_FILE_NAME",c_val=mc_par%displacement_file,error=error) 01310 CALL section_vals_val_get(mc_section,"BOX2_FILE_NAME",c_val=mc_par%box2_file,error=error) 01311 01312 ! set the values of the arrays...if we just point, we have problems when we start fooling around 01313 ! with releasing force_envs and wonky values get set (despite that these are private) 01314 IF(mc_par%ensemble == "VIRIAL") THEN 01315 CALL section_vals_val_get(mc_section,"VIRIAL_TEMPS",r_vals=temp_r_array,error=error) 01316 ! yes, I'm allocating here...I cannot find a better place to do it, though 01317 ALLOCATE (mc_par%virial_temps(1:SIZE(temp_r_array)),STAT=istat) 01318 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01319 "virial_temps",dp_size*SIZE(temp_r_array)) 01320 DO imol=1,SIZE(temp_r_array) 01321 mc_par%virial_temps(imol)=temp_r_array(imol) 01322 ENDDO 01323 END IF 01324 ! all of these arrays should have one value for each type of molecule...so check that 01325 CALL section_vals_val_get(mc_section,"AVBMC_RMIN",r_vals=temp_r_array,error=error) 01326 IF( .NOT. mc_par%lhmc)THEN 01327 IF(SIZE(temp_r_array) .NE. nmol_types)THEN 01328 CALL stop_program(routineN,moduleN,__LINE__,& 01329 'AVBMC_RMIN must have as many values as there are molecule types.') 01330 ENDIF 01331 ENDIF 01332 DO imol=1,SIZE(temp_r_array) 01333 mc_par%avbmc_rmin(imol)=temp_r_array(imol) 01334 ENDDO 01335 CALL section_vals_val_get(mc_section,"AVBMC_RMAX",r_vals=temp_r_array,error=error) 01336 IF( .NOT. mc_par%lhmc)THEN 01337 IF(SIZE(temp_r_array) .NE. nmol_types)THEN 01338 CALL stop_program(routineN,moduleN,__LINE__,& 01339 'AVBMC_RMAXL must have as many values as there are molecule types.') 01340 ENDIF 01341 ENDIF 01342 DO imol=1,SIZE(temp_r_array) 01343 mc_par%avbmc_rmax(imol)=temp_r_array(imol) 01344 ENDDO 01345 CALL section_vals_val_get(mc_section,"PBIAS",r_vals=temp_r_array,error=error) 01346 IF( .NOT. mc_par%lhmc)THEN 01347 IF(SIZE(temp_r_array) .NE. nmol_types)THEN 01348 CALL stop_program(routineN,moduleN,__LINE__,& 01349 'PBIAS must have as many values as there are molecule types.') 01350 ENDIF 01351 ENDIF 01352 DO imol=1,SIZE(temp_r_array) 01353 mc_par%pbias(imol)=temp_r_array(imol) 01354 ENDDO 01355 CALL section_vals_val_get(mc_section,"PMAVBMC_MOL",r_vals=temp_r_array,error=error) 01356 IF( .NOT. mc_par%lhmc)THEN 01357 IF(SIZE(temp_r_array) .NE. nmol_types)THEN 01358 CALL stop_program(routineN,moduleN,__LINE__,& 01359 'PMAVBMC_MOL must have as many values as there are molecule types.') 01360 ENDIF 01361 ENDIF 01362 DO imol=1,SIZE(temp_r_array) 01363 mc_par%pmavbmc_mol(imol)=temp_r_array(imol) 01364 ENDDO 01365 CALL section_vals_val_get(mc_section,"ETA",r_vals=temp_r_array,error=error) 01366 IF( .NOT. mc_par%lhmc)THEN 01367 IF(SIZE(temp_r_array) .NE. nmol_types)THEN 01368 CALL stop_program(routineN,moduleN,__LINE__,& 01369 'ETA must have as many values as there are molecule types.') 01370 ENDIF 01371 ENDIF 01372 DO imol=1,SIZE(temp_r_array) 01373 mc_par%eta(imol)=temp_r_array(imol) 01374 ENDDO 01375 CALL section_vals_val_get(mc_section,"RMBOND",r_vals=temp_r_array,error=error) 01376 IF( .NOT. mc_par%lhmc)THEN 01377 IF(SIZE(temp_r_array) .NE. nmol_types)THEN 01378 CALL stop_program(routineN,moduleN,__LINE__,& 01379 'RMBOND must have as many values as there are molecule types.') 01380 ENDIF 01381 ENDIF 01382 DO imol=1,SIZE(temp_r_array) 01383 mc_par%rmbond(imol)=temp_r_array(imol) 01384 ENDDO 01385 CALL section_vals_val_get(mc_section,"RMANGLE",r_vals=temp_r_array,error=error) 01386 IF( .NOT. mc_par%lhmc)THEN 01387 IF(SIZE(temp_r_array) .NE. nmol_types)THEN 01388 CALL stop_program(routineN,moduleN,__LINE__,& 01389 'RMANGLE must have as many values as there are molecule types.') 01390 ENDIF 01391 ENDIF 01392 DO imol=1,SIZE(temp_r_array) 01393 mc_par%rmangle(imol)=temp_r_array(imol) 01394 ENDDO 01395 CALL section_vals_val_get(mc_section,"RMDIHEDRAL",r_vals=temp_r_array,error=error) 01396 IF( .NOT. mc_par%lhmc)THEN 01397 IF(SIZE(temp_r_array) .NE. nmol_types)THEN 01398 CALL stop_program(routineN,moduleN,__LINE__,& 01399 'RMDIHEDRAL must have as many values as there are molecule types.') 01400 ENDIF 01401 ENDIF 01402 DO imol=1,SIZE(temp_r_array) 01403 mc_par%rmdihedral(imol)=temp_r_array(imol) 01404 ENDDO 01405 CALL section_vals_val_get(mc_section,"RMROT",r_vals=temp_r_array,error=error) 01406 IF( .NOT. mc_par%lhmc)THEN 01407 IF(SIZE(temp_r_array) .NE. nmol_types)THEN 01408 CALL stop_program(routineN,moduleN,__LINE__,& 01409 'RMROT must have as many values as there are molecule types.') 01410 ENDIF 01411 ENDIF 01412 DO imol=1,SIZE(temp_r_array) 01413 mc_par%rmrot(imol)=temp_r_array(imol) 01414 ENDDO 01415 CALL section_vals_val_get(mc_section,"RMTRANS",r_vals=temp_r_array,error=error) 01416 IF( .NOT. mc_par%lhmc)THEN 01417 IF(SIZE(temp_r_array) .NE. nmol_types)THEN 01418 CALL stop_program(routineN,moduleN,__LINE__,& 01419 'RMTRANS must have as many values as there are molecule types.') 01420 ENDIF 01421 ENDIF 01422 DO imol=1,SIZE(temp_r_array) 01423 mc_par%rmtrans(imol)=temp_r_array(imol) 01424 ENDDO 01425 01426 CALL section_vals_val_get(mc_section,"PMTRAION_MOL",r_vals=temp_r_array,error=error) 01427 IF( .NOT. mc_par%lhmc)THEN 01428 IF(SIZE(temp_r_array) .NE. nmol_types)THEN 01429 CALL stop_program(routineN,moduleN,__LINE__,& 01430 'PMTRAION_MOL must have as many values as there are molecule types.') 01431 ENDIF 01432 ENDIF 01433 DO imol=1,SIZE(temp_r_array) 01434 mc_par%pmtraion_mol(imol)=temp_r_array(imol) 01435 ENDDO 01436 01437 CALL section_vals_val_get(mc_section,"PMTRANS_MOL",r_vals=temp_r_array,error=error) 01438 IF( .NOT. mc_par%lhmc)THEN 01439 IF(SIZE(temp_r_array) .NE. nmol_types)THEN 01440 CALL stop_program(routineN,moduleN,__LINE__,& 01441 'PMTRANS_MOL must have as many values as there are molecule types.') 01442 ENDIF 01443 ENDIF 01444 DO imol=1,SIZE(temp_r_array) 01445 mc_par%pmtrans_mol(imol)=temp_r_array(imol) 01446 ENDDO 01447 01448 CALL section_vals_val_get(mc_section,"PMROT_MOL",r_vals=temp_r_array,error=error) 01449 IF( .NOT. mc_par%lhmc)THEN 01450 IF(SIZE(temp_r_array) .NE. nmol_types)THEN 01451 CALL stop_program(routineN,moduleN,__LINE__,& 01452 'PMROT_MOL must have as many values as there are molecule types.') 01453 ENDIF 01454 ENDIF 01455 DO imol=1,SIZE(temp_r_array) 01456 mc_par%pmrot_mol(imol)=temp_r_array(imol) 01457 ENDDO 01458 01459 CALL section_vals_val_get(mc_section,"PMSWAP_MOL",r_vals=temp_r_array,error=error) 01460 IF( .NOT. mc_par%lhmc)THEN 01461 IF(SIZE(temp_r_array) .NE. nmol_types)THEN 01462 CALL stop_program(routineN,moduleN,__LINE__,& 01463 'PMSWAP_MOL must have as many values as there are molecule types.') 01464 ENDIF 01465 ENDIF 01466 DO imol=1,SIZE(temp_r_array) 01467 mc_par%pmswap_mol(imol)=temp_r_array(imol) 01468 ENDDO 01469 01470 ! now some logical values 01471 CALL section_vals_val_get(mc_section,"LBIAS",l_val=mc_par%lbias,error=error) 01472 CALL section_vals_val_get(mc_section,"LDISCRETE",l_val=mc_par%ldiscrete,error=error) 01473 CALL section_vals_val_get(mc_section,"LSTOP",l_val=mc_par%lstop,error=error) 01474 CALL section_vals_val_get(mc_section,"RESTART",l_val=mc_par%lrestart,error=error) 01475 CALL section_vals_val_get(mc_section,"LBIAS",l_val=mc_par%lbias,error=error) 01476 01477 !..end of parsing the input section 01478 01479 !..write some information to output 01480 IF (iw>0) THEN 01481 WRITE(box_string,'(I2)') mc_par%mc_molecule_info%nboxes 01482 WRITE(molecule_string,'(I2)') mc_par%mc_molecule_info%nmol_types 01483 WRITE(tab_string,'(I4)') 81-10*mc_par%mc_molecule_info%nmol_types 01484 WRITE(tab_box_string,'(I4)') 81-10*mc_par%mc_molecule_info%nboxes 01485 WRITE(tab_string_int,'(I4)') 81-5*mc_par%mc_molecule_info%nmol_types 01486 format_string="(A,T" // TRIM(ADJUSTL(tab_string)) // "," // TRIM(ADJUSTL(molecule_string)) // "(2X,F8.4))" 01487 format_box_string="(A,T" // TRIM(ADJUSTL(tab_box_string)) // "," // TRIM(ADJUSTL(box_string)) // "(2X,F8.4))" 01488 format_string_int="(A,T" // TRIM(ADJUSTL(tab_string)) // "," // TRIM(ADJUSTL(molecule_string)) // "(I3,2X))" 01489 WRITE ( iw, '( A )' ) ' MC| Monte Carlo Protocol ' 01490 WRITE ( iw, '( A,T71,I10 )' ) ' MC| total number of steps ', & 01491 mc_par % nstep 01492 WRITE ( iw, '( A,T71,F10.3 )' ) ' MC| pmvolume ', & 01493 mc_par % pmvolume 01494 WRITE ( iw, '( A,T71,F10.3 )' ) ' MC| pmvol_box ', & 01495 mc_par % pmvol_box 01496 WRITE ( iw, '( A,T71,F10.3 )' ) ' MC| pmhmc ', & 01497 mc_par % pmhmc 01498 WRITE ( iw, '( A,T71,F10.3 )' ) ' MC| pmhmc_box ', & 01499 mc_par % pmhmc_box 01500 WRITE ( iw, '( A,T71,F10.3 )' ) ' MC| pmswap ', & 01501 mc_par % pmswap 01502 WRITE ( iw, '( A,T71,F10.3 )' ) ' MC| pmavbmc ', & 01503 mc_par % pmavbmc 01504 WRITE ( iw, '( A,T71,F10.3 )' ) ' MC| pmtraion ', & 01505 mc_par % pmtraion 01506 WRITE ( iw, '( A,T71,F10.3 )' ) ' MC| pmtrans ', & 01507 mc_par % pmtrans 01508 WRITE ( iw, '( A,T71,I10 )' ) ' MC| iupvolume ', & 01509 mc_par % iupvolume 01510 WRITE ( iw, '( A,T71,I10 )' ) ' MC| nvirial ', & 01511 mc_par % nvirial 01512 WRITE ( iw, '( A,T71,I10 )' ) ' MC| iuptrans ', & 01513 mc_par % iuptrans 01514 WRITE ( iw, '( A,T71,I10 )' ) ' MC| iprint ', & 01515 mc_par % iprint 01516 WRITE ( iw, '( A,T61,A20 )' ) ' MC| ensemble ', & 01517 ADJUSTR(mc_par % ensemble) 01518 ! format string may not be valid if all the molecules are atomic, 01519 ! like in HMC 01520 IF(.NOT. mc_par % lhmc)THEN 01521 WRITE ( iw, format_string ) ' MC| pmswap_mol ', & 01522 mc_par % pmswap_mol 01523 WRITE ( iw, format_string ) ' MC| pmavbmc_mol ', & 01524 mc_par % pmavbmc_mol 01525 WRITE ( iw, format_string ) ' MC| pbias ', & 01526 mc_par % pbias 01527 WRITE ( iw, format_string ) ' MC| pmtraion_mol ', & 01528 mc_par % pmtraion_mol 01529 WRITE ( iw, format_string ) ' MC| pmtrans_mol ', & 01530 mc_par % pmtrans_mol 01531 WRITE ( iw, format_string ) ' MC| pmrot_mol ', & 01532 mc_par % pmrot_mol 01533 WRITE ( iw, format_string ) ' MC| eta [K]', & 01534 mc_par % eta 01535 WRITE ( iw, format_string ) ' MC| rmbond [angstroms]', & 01536 mc_par % rmbond 01537 WRITE ( iw, format_string ) ' MC| rmangle [degrees]', & 01538 mc_par % rmangle 01539 WRITE ( iw, format_string ) ' MC| rmdihedral [degrees]', & 01540 mc_par % rmdihedral 01541 WRITE ( iw, format_string ) ' MC| rmtrans [angstroms]', & 01542 mc_par % rmtrans 01543 WRITE ( iw, format_string ) ' MC| rmrot [degrees]', & 01544 mc_par % rmrot 01545 WRITE ( iw, format_string_int ) ' MC| AVBMC target atom ', & 01546 mc_par % avbmc_atom 01547 WRITE ( iw, format_string ) ' MC| AVBMC inner cutoff [ang]', & 01548 mc_par % avbmc_rmin 01549 WRITE ( iw, format_string ) ' MC| AVBMC outer cutoff [ang]', & 01550 mc_par % avbmc_rmax 01551 ENDIF 01552 IF (mc_par%ensemble .EQ. 'GEMC-NVT') THEN 01553 WRITE ( iw, '( A,T58,A)' ) ' MC| Box 2 file', & 01554 TRIM(mc_par % box2_file) 01555 ENDIF 01556 WRITE ( iw, '( A,T58,A )' ) ' MC| Name of restart file:',& 01557 TRIM(mc_par % restart_file_name) 01558 WRITE ( iw, '( A,A,T44,A )' ) ' MC| Name of output ',& 01559 'coordinate file:',& 01560 TRIM(mc_par % coords_file) 01561 WRITE ( iw, '( A,T44,A )' ) ' MC| Name of output data file:',& 01562 TRIM(mc_par % data_file) 01563 WRITE ( iw, '( A,A,T44,A )' ) ' MC| Name of output ',& 01564 'molecules file:',& 01565 TRIM(mc_par %molecules_file) 01566 WRITE ( iw, '( A,T44,A )' ) ' MC| Name of output moves file:',& 01567 TRIM(mc_par % moves_file) 01568 WRITE ( iw, '( A,T44,A )') ' MC| Name of output energy file:',& 01569 TRIM(mc_par % energy_file) 01570 WRITE ( iw, '( A,T44,A )' ) ' MC| Name of output cell file:',& 01571 TRIM(mc_par % cell_file) 01572 WRITE ( iw, '( A,A,T44,A )' ) ' MC| Name of output',& 01573 ' displacement file:',& 01574 TRIM(mc_par % displacement_file) 01575 IF(mc_par % ldiscrete) THEN 01576 WRITE ( iw, '(A,A,T71,F10.3)' ) ' MC| discrete step size',& 01577 '[angstroms]', & 01578 mc_par % discrete_step 01579 ELSE 01580 WRITE ( iw, '( A,A,T71,F10.3 )' ) ' MC| rmvolume ',& 01581 '[cubic angstroms]', & 01582 mc_par % rmvolume 01583 ENDIF 01584 WRITE ( iw, '( A,T71,F10.2 )' ) ' MC| Temperature [K] ', & 01585 mc_par % temperature 01586 WRITE ( iw, '( A,T71,F10.5 )' ) ' MC| Pressure [bar] ', & 01587 mc_par % pressure 01588 IF ( mc_par % lrestart ) THEN 01589 WRITE ( iw, '(A,A)') ' MC| Initial data will be read from a',& 01590 ' restart file.' 01591 ENDIF 01592 IF ( mc_par % lbias ) THEN 01593 WRITE ( iw, '(A)') ' MC| The moves will be biased.' 01594 ELSE 01595 WRITE ( iw, '(A)') ' MC| The moves will not be biased,' 01596 ENDIF 01597 IF (mc_par%nmoves .EQ. 1) THEN 01598 WRITE(iw,'(A,A)') ' MC| A full energy calculation ',& 01599 'will be done at every step.' 01600 ELSE 01601 WRITE( iw, '(A,I4,A,A)' ) ' MC| ',mc_par%nmoves,& 01602 ' moves will be attempted ',& 01603 'before a Quickstep energy calculation' 01604 WRITE( iw, '(A)' ) ' MC| takes place.' 01605 ENDIF 01606 01607 WRITE( iw, '(A,I4,A,A)') ' MC| ',mc_par%nswapmoves,& 01608 ' swap insertions will be attempted ',& 01609 'per molecular swap move' 01610 01611 IF(mc_par % lhmc)THEN 01612 WRITE(iw,'(A)') '********************************************************' 01613 WRITE(iw,'(A)') '************** ONLY DOING HYBRID MONTE CARLO MOVES **************************' 01614 WRITE(iw,'(A)') '********************************************************' 01615 01616 ENDIF 01617 END IF 01618 01619 ! figure out what beta (1/kT) is in atomic units (1/Hartree) 01620 mc_par % BETA = 1 / mc_par%temperature / boltzmann * joule 01621 01622 ! 0.9_dp is to avoid overflow or underflow 01623 mc_par % exp_max_val = 0.9_dp*LOG(HUGE(0.0_dp)) 01624 mc_par % exp_min_val = 0.9_dp*LOG(TINY(0.0_dp)) 01625 mc_par % max_val = EXP(mc_par % exp_max_val) 01626 mc_par % min_val = 0.0_dp 01627 01628 ! convert from bar to a.u. 01629 mc_par%pressure = cp_unit_to_cp2k(mc_par%pressure,"bar",error=error) 01630 ! convert from angstrom to a.u. 01631 DO itype=1,mc_par%mc_molecule_info%nmol_types 01632 ! convert from Kelvin to a.u. 01633 ! mc_par % eta = mc_par%eta(itype) * boltzmann / joule 01634 ! convert from degrees to radians 01635 mc_par%rmrot(itype) = mc_par%rmrot(itype)/180.0e0_dp*pi 01636 mc_par%rmangle(itype) = mc_par%rmangle(itype)/180.0e0_dp*pi 01637 mc_par%rmdihedral(itype) = mc_par%rmdihedral(itype)/180.0e0_dp*pi 01638 mc_par%rmtrans(itype) = cp_unit_to_cp2k(mc_par%rmtrans(itype),"angstrom",error=error) 01639 mc_par%rmbond(itype) = cp_unit_to_cp2k(mc_par%rmbond(itype),"angstrom",error=error) 01640 mc_par%eta(itype) = cp_unit_to_cp2k(mc_par%eta(itype),"K_e",error=error) 01641 mc_par%avbmc_rmin(itype) = cp_unit_to_cp2k(mc_par%avbmc_rmin(itype),"angstrom",error=error) 01642 mc_par%avbmc_rmax(itype) = cp_unit_to_cp2k(mc_par%avbmc_rmax(itype),"angstrom",error=error) 01643 ENDDO 01644 mc_par%rmvolume = cp_unit_to_cp2k(mc_par%rmvolume,"angstrom^3",error=error) 01645 mc_par%discrete_step = cp_unit_to_cp2k(mc_par%discrete_step,"angstrom",error=error) 01646 01647 ! end the timing 01648 CALL timestop(handle) 01649 01650 END SUBROUTINE read_mc_section 01651 01652 ! ***************************************************************************** 01664 SUBROUTINE find_mc_rcut ( mc_par, force_env, lterminate ) 01665 01666 TYPE(mc_simpar_type), INTENT(INOUT) :: mc_par 01667 TYPE(force_env_type), POINTER :: force_env 01668 LOGICAL, INTENT(OUT) :: lterminate 01669 01670 CHARACTER(len=*), PARAMETER :: routineN = 'find_mc_rcut', 01671 routineP = moduleN//':'//routineN 01672 01673 INTEGER :: itype, jtype 01674 REAL(KIND=dp) :: rcutsq_max 01675 REAL(KIND=dp), DIMENSION(1:3) :: abc 01676 TYPE(cell_type), POINTER :: cell 01677 TYPE(cp_error_type) :: error 01678 TYPE(fist_environment_type), POINTER :: fist_env 01679 TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env 01680 TYPE(pair_potential_pp_type), POINTER :: potparm 01681 01682 NULLIFY(cell,potparm,fist_nonbond_env,fist_env) 01683 01684 lterminate=.FALSE. 01685 CALL force_env_get(force_env,fist_env=fist_env,error=error) 01686 CALL fist_env_get(fist_env,cell=cell,& 01687 fist_nonbond_env=fist_nonbond_env,error=error) 01688 CALL fist_nonbond_env_get (fist_nonbond_env, potparm=potparm,error=error) 01689 CALL get_cell(cell,abc=abc) 01690 01691 ! find the largest value of rcutsq 01692 rcutsq_max=0.0e0_dp 01693 DO itype=1,SIZE(potparm%pot,1) 01694 DO jtype=itype,SIZE(potparm%pot,2) 01695 IF(potparm%pot(itype,jtype)%pot%rcutsq .GT. rcutsq_max) & 01696 rcutsq_max=potparm%pot(itype,jtype)%pot%rcutsq 01697 ENDDO 01698 ENDDO 01699 01700 ! check to make sure all box dimensions are greater than two times this 01701 ! value 01702 mc_par % rcut=rcutsq_max**0.5_dp 01703 DO itype=1,3 01704 IF(abc(itype) .LT. 2.0_dp*mc_par % rcut) THEN 01705 lterminate=.TRUE. 01706 ENDIF 01707 ENDDO 01708 01709 END SUBROUTINE find_mc_rcut 01710 01711 ! ***************************************************************************** 01724 SUBROUTINE mc_determine_molecule_info(force_env,mc_molecule_info,error,box_number,& 01725 coordinates_empty) 01726 01727 TYPE(force_env_p_type), DIMENSION(:), 01728 POINTER :: force_env 01729 TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info 01730 TYPE(cp_error_type), INTENT(inout) :: error 01731 INTEGER, INTENT(IN), OPTIONAL :: box_number 01732 REAL(dp), DIMENSION(:, :), OPTIONAL, 01733 POINTER :: coordinates_empty 01734 01735 CHARACTER(LEN=*), PARAMETER :: routineN = 'mc_determine_molecule_info', 01736 routineP = moduleN//':'//routineN 01737 01738 CHARACTER(LEN=default_string_length) :: name 01739 CHARACTER(LEN=default_string_length), 01740 ALLOCATABLE, DIMENSION(:) :: names_init 01741 INTEGER :: first_mol, iatom, ibox, imol, imolecule, istat, itype, 01742 iunique, iunit, jtype, last_mol, natom, natoms_large, nbend, nbond, 01743 nboxes, nmol_types, nmolecule, ntorsion, ntypes, skip_box, 01744 this_molecule, total 01745 INTEGER, DIMENSION(:), POINTER :: molecule_list 01746 LOGICAL :: lnew_type 01747 TYPE(atom_type), DIMENSION(:), POINTER :: atom_list 01748 TYPE(atomic_kind_type), POINTER :: atomic_kind 01749 TYPE(cp_subsys_type), POINTER :: subsys 01750 TYPE(mol_kind_new_list_p_type), 01751 DIMENSION(:), POINTER :: molecule_kinds_new 01752 TYPE(molecule_kind_type), POINTER :: molecule_kind 01753 TYPE(particle_list_p_type), 01754 DIMENSION(:), POINTER :: particles 01755 01756 ! how many simulation boxes do we have? 01757 01758 nboxes=SIZE(force_env(:)) 01759 01760 ! allocate the pointer 01761 ALLOCATE(mc_molecule_info) 01762 mc_molecule_info%nboxes=nboxes 01763 01764 ! if box_number is present, that box is supposed to be empty, 01765 ! so we don't count it in anything 01766 IF(PRESENT(box_number)) THEN 01767 skip_box=box_number 01768 ELSE 01769 skip_box=0 01770 ENDIF 01771 01772 ! we need to figure out how many different kinds of molecules we have... 01773 ! the fun stuff comes in when box 1 is missing a molecule type...I'll 01774 ! distinguish molecules based on names from the .psf files... 01775 ! this is only really a problem when we have more than one simulation 01776 ! box 01777 NULLIFY(subsys,molecule_kinds_new,molecule_kind) 01778 01779 ALLOCATE(particles(1:nboxes),STAT=istat) 01780 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01781 "particles",nboxes) 01782 ALLOCATE (molecule_kinds_new(1:nboxes),STAT=istat) 01783 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01784 "molecule_kinds_new",nboxes) 01785 01786 ! find out how many types we have 01787 ntypes=0 01788 DO ibox=1,nboxes 01789 IF(ibox == skip_box) CYCLE 01790 CALL force_env_get(force_env(ibox)%force_env,& 01791 subsys=subsys,error=error) 01792 CALL cp_subsys_get(subsys, & 01793 molecule_kinds_new=molecule_kinds_new(ibox)%list,& 01794 particles=particles(ibox)%list,error=error) 01795 ntypes=ntypes+SIZE(molecule_kinds_new(ibox)%list%els(:)) 01796 ENDDO 01797 01798 ALLOCATE (names_init(1:ntypes),STAT=istat) 01799 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01800 "names_init",ntypes*default_string_length) 01801 01802 ! now get the names of all types so we can figure out how many distinct 01803 ! types we have 01804 itype=1 01805 natoms_large=0 01806 DO ibox=1,nboxes 01807 IF(ibox == skip_box) CYCLE 01808 DO imolecule=1,SIZE(molecule_kinds_new(ibox)%list%els(:)) 01809 molecule_kind=> molecule_kinds_new(ibox)%list%els(imolecule) 01810 CALL get_molecule_kind(molecule_kind,name=names_init(itype),& 01811 natom=natom) 01812 IF(natom .GT. natoms_large) natoms_large=natom 01813 itype=itype+1 01814 ENDDO 01815 ENDDO 01816 01817 nmol_types=0 01818 DO itype=1,ntypes 01819 lnew_type=.TRUE. 01820 DO jtype=1,itype-1 01821 IF(TRIM(names_init(itype)) .EQ. TRIM(names_init(jtype)))& 01822 lnew_type=.FALSE. 01823 ENDDO 01824 IF(lnew_type) THEN 01825 nmol_types=nmol_types+1 01826 ELSE 01827 names_init(itype)='' 01828 ENDIF 01829 ENDDO 01830 01831 ! allocate arrays for the names of the molecules, the number of atoms, and 01832 ! the conformational probabilities 01833 mc_molecule_info%nmol_types=nmol_types 01834 ALLOCATE(mc_molecule_info%names(1:nmol_types),STAT=istat) 01835 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01836 "mc_molecule_info%names",nmol_types*default_string_length) 01837 ALLOCATE(mc_molecule_info%nunits(1:nmol_types),STAT=istat) 01838 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01839 "mc_molecule_info%nunits",nmol_types*int_size) 01840 ALLOCATE(mc_molecule_info%conf_prob(1:3,1:nmol_types),STAT=istat) 01841 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01842 "mc_molecule_info%conf_prob",3*nmol_types*dp_size) 01843 ALLOCATE(mc_molecule_info%nchains(1:nmol_types,1:nboxes),STAT=istat) 01844 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01845 "mc_molecule_info%nchains",nboxes*nmol_types*dp_size) 01846 ALLOCATE(mc_molecule_info%nunits_tot(1:nboxes),STAT=istat) 01847 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01848 "mc_molecule_info%nunits_tot",nboxes*int_size) 01849 ALLOCATE(mc_molecule_info%atom_names(1:natoms_large,1:nmol_types),STAT=istat) 01850 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01851 "mc_molecule_info%atom_names",natoms_large*nmol_types*default_string_length) 01852 ALLOCATE(mc_molecule_info%mass(1:natoms_large,1:nmol_types),STAT=istat) 01853 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01854 "mc_molecule_info%mass",natoms_large*nmol_types*dp_size) 01855 01856 ! now record all the information for every molecule type 01857 iunique=0 01858 itype=0 01859 DO ibox=1,nboxes 01860 IF(ibox == skip_box) CYCLE 01861 DO imolecule=1,SIZE(molecule_kinds_new(ibox)%list%els(:)) 01862 itype=itype+1 01863 IF(names_init(itype) .EQ. '') CYCLE 01864 iunique=iunique+1 01865 mc_molecule_info%names(iunique)=names_init(itype) 01866 molecule_kind=> molecule_kinds_new(ibox)%list%els(imolecule) 01867 CALL get_molecule_kind(molecule_kind,natom=mc_molecule_info%nunits(iunique),& 01868 nbond=nbond,nbend=nbend,ntorsion=ntorsion,atom_list=atom_list) 01869 01870 DO iatom=1,mc_molecule_info%nunits(iunique) 01871 atomic_kind => atom_list(iatom)%atomic_kind 01872 CALL get_atomic_kind(atomic_kind=atomic_kind,& 01873 name=mc_molecule_info%atom_names(iatom,iunique),& 01874 mass=mc_molecule_info%mass(iatom,iunique)) 01875 ENDDO 01876 01877 ! compute the probabilities of doing any particular kind of conformation change 01878 total=nbond+nbend+ntorsion 01879 01880 IF(total == 0) THEN 01881 mc_molecule_info%conf_prob(:,iunique)=0.0e0_dp 01882 ELSE 01883 mc_molecule_info%conf_prob(1,iunique)=REAL(nbond,dp)/REAL(total,dp) 01884 mc_molecule_info%conf_prob(2,iunique)=REAL(nbend,dp)/REAL(total,dp) 01885 mc_molecule_info%conf_prob(3,iunique)=REAL(ntorsion,dp)/REAL(total,dp) 01886 ENDIF 01887 01888 ENDDO 01889 ENDDO 01890 01891 ! figure out the empty coordinates 01892 IF(PRESENT(coordinates_empty)) THEN 01893 ALLOCATE(coordinates_empty(1:3,1:mc_molecule_info%nunits(1)),STAT=istat) 01894 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01895 "coordinates_empty",3*mc_molecule_info%nunits(1)*dp_size) 01896 DO iunit=1,mc_molecule_info%nunits(1) 01897 coordinates_empty(1:3,iunit)=particles(1)%list%els(iunit)%r(1:3) 01898 ENDDO 01899 ENDIF 01900 01901 ! now we need to figure out the total number of molecules of each kind in 01902 ! each box, and the total number of interaction sites in each box 01903 mc_molecule_info%nchains(:,:)=0 01904 DO iunique=1,nmol_types 01905 DO ibox=1,nboxes 01906 IF(ibox == skip_box) CYCLE 01907 DO imolecule=1,SIZE(molecule_kinds_new(ibox)%list%els(:)) 01908 molecule_kind=> molecule_kinds_new(ibox)%list%els(imolecule) 01909 CALL get_molecule_kind(molecule_kind,nmolecule=nmolecule,& 01910 name=name) 01911 IF(TRIM(name) .NE. mc_molecule_info%names(iunique) ) CYCLE 01912 mc_molecule_info%nchains(iunique,ibox)=mc_molecule_info%nchains(iunique,ibox)+nmolecule 01913 ENDDO 01914 ENDDO 01915 ENDDO 01916 mc_molecule_info%nchain_total=0 01917 DO ibox=1,nboxes 01918 mc_molecule_info%nunits_tot(ibox)=0 01919 IF(ibox == skip_box) CYCLE 01920 DO iunique=1,nmol_types 01921 mc_molecule_info%nunits_tot(ibox)=mc_molecule_info%nunits_tot(ibox)+& 01922 mc_molecule_info%nunits(iunique)*mc_molecule_info%nchains(iunique,ibox) 01923 ENDDO 01924 mc_molecule_info%nchain_total=mc_molecule_info%nchain_total+SUM(mc_molecule_info%nchains(:,ibox)) 01925 ENDDO 01926 01927 ! now we need to figure out which type every molecule is, 01928 ! and which force_env (box) it's in 01929 ALLOCATE(mc_molecule_info%mol_type(1:mc_molecule_info%nchain_total),STAT=istat) 01930 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01931 "mol_type",mc_molecule_info%nchain_total*int_size) 01932 ALLOCATE(mc_molecule_info%in_box(1:mc_molecule_info%nchain_total),STAT=istat) 01933 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01934 "in_box",mc_molecule_info%nchain_total*int_size) 01935 01936 last_mol=0 01937 DO ibox=1,nboxes 01938 IF(ibox == skip_box) CYCLE 01939 first_mol=last_mol+1 01940 last_mol=first_mol+SUM(mc_molecule_info%nchains(:,ibox))-1 01941 mc_molecule_info%in_box(first_mol:last_mol)=ibox 01942 DO imolecule=1,SIZE(molecule_kinds_new(ibox)%list%els(:)) 01943 molecule_kind=> molecule_kinds_new(ibox)%list%els(imolecule) 01944 CALL get_molecule_kind(molecule_kind,nmolecule=nmolecule,& 01945 name=name,molecule_list=molecule_list) 01946 ! figure out which molecule number this is 01947 this_molecule=0 01948 DO iunique=1,nmol_types 01949 IF(TRIM(name) .EQ. TRIM(mc_molecule_info%names(iunique))) THEN 01950 this_molecule=iunique 01951 ENDIF 01952 ENDDO 01953 DO imol=1,SIZE(molecule_list(:)) 01954 mc_molecule_info%mol_type(first_mol+molecule_list(imol)-1)=this_molecule 01955 ENDDO 01956 ENDDO 01957 ENDDO 01958 01959 ! get rid of stuff we don't need 01960 DEALLOCATE(names_init,stat=istat) 01961 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"names_init") 01962 DEALLOCATE(molecule_kinds_new,stat=istat) 01963 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"molecule_kinds_new") 01964 DEALLOCATE(particles,stat=istat) 01965 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"particles") 01966 01967 END SUBROUTINE MC_DETERMINE_MOLECULE_INFO 01968 01969 ! ***************************************************************************** 01976 SUBROUTINE mc_molecule_info_destroy(mc_molecule_info) 01977 01978 TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info 01979 01980 CHARACTER(LEN=*), PARAMETER :: routineN = 'mc_molecule_info_destroy', 01981 routineP = moduleN//':'//routineN 01982 01983 INTEGER :: istat 01984 01985 DEALLOCATE(mc_molecule_info%nchains,stat=istat) 01986 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"nchains") 01987 DEALLOCATE(mc_molecule_info%nunits,stat=istat) 01988 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"nunits") 01989 DEALLOCATE(mc_molecule_info%mol_type,stat=istat) 01990 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"mol_type") 01991 DEALLOCATE(mc_molecule_info%in_box,stat=istat) 01992 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"in_box") 01993 DEALLOCATE(mc_molecule_info%names,stat=istat) 01994 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"names") 01995 DEALLOCATE(mc_molecule_info%atom_names,stat=istat) 01996 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"atom_names") 01997 DEALLOCATE(mc_molecule_info%conf_prob,stat=istat) 01998 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"conf_prob") 01999 DEALLOCATE(mc_molecule_info%nunits_tot,stat=istat) 02000 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"nunits_tot") 02001 DEALLOCATE(mc_molecule_info%mass,stat=istat) 02002 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"mass") 02003 02004 DEALLOCATE(mc_molecule_info,stat=istat) 02005 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"mc_molecule_info") 02006 NULLIFY(mc_molecule_info) 02007 02008 END SUBROUTINE mc_molecule_info_destroy 02009 02010 ! ***************************************************************************** 02011 SUBROUTINE mc_input_parameters_check ( mc_par) 02012 02013 ! ***************************************************************************** 02020 TYPE(mc_simpar_type), POINTER :: mc_par 02021 02022 CHARACTER(len=*), PARAMETER :: routineN = 'mc_input_parameters_check', 02023 routineP = moduleN//':'//routineN 02024 02025 INTEGER :: imol_type, nboxes, 02026 nchain_total, nmol_types 02027 INTEGER, DIMENSION(:), POINTER :: nunits 02028 LOGICAL :: lhmc 02029 TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info 02030 02031 CALL get_mc_par(mc_par,mc_molecule_info=mc_molecule_info,lhmc=lhmc) 02032 02033 CALL get_mc_molecule_info(mc_molecule_info,nmol_types=nmol_types,& 02034 nboxes=nboxes,nunits=nunits,nchain_total=nchain_total) 02035 02036 ! if we're only doing HMC, we can skip these checks 02037 IF(lhmc) RETURN 02038 02039 ! the final value of these arrays needs to be 1.0, otherwise there is 02040 ! a chance we won't select a molecule type for a move 02041 IF(mc_par%pmavbmc_mol(nmol_types) .LT. 0.9999E0_dp) THEN 02042 CALL stop_program(routineN,moduleN,__LINE__,& 02043 'The last value of PMAVBMC_MOL needs to be 1.0') 02044 ENDIF 02045 IF(mc_par%pmswap_mol(nmol_types) .LT. 0.9999E0_dp) THEN 02046 CALL stop_program(routineN,moduleN,__LINE__,& 02047 'The last value of PMSWAP_MOL needs to be 1.0') 02048 ENDIF 02049 IF(mc_par%pmtraion_mol(nmol_types) .LT. 0.9999E0_dp) THEN 02050 CALL stop_program(routineN,moduleN,__LINE__,& 02051 'The last value of PMTRAION_MOL needs to be 1.0') 02052 ENDIF 02053 IF(mc_par%pmtrans_mol(nmol_types) .LT. 0.9999E0_dp) THEN 02054 CALL stop_program(routineN,moduleN,__LINE__,& 02055 'The last value of PMTRANS_MOL needs to be 1.0') 02056 ENDIF 02057 IF(mc_par%pmrot_mol(nmol_types) .LT. 0.9999E0_dp) THEN 02058 CALL stop_program(routineN,moduleN,__LINE__,& 02059 'The last value of PMROT_MOL needs to be 1.0') 02060 ENDIF 02061 02062 ! check to make sure we have all the desired values for various 02063 02064 ! some ensembles require multiple boxes or molecule types 02065 SELECT CASE(mc_par%ensemble) 02066 CASE("GEMC_NPT") 02067 IF(nmol_types .LE. 1) & 02068 CALL stop_program(routineN,moduleN,__LINE__,& 02069 'Cannot have GEMC-NPT simulation with only one molecule type') 02070 IF(nboxes .LE. 1) & 02071 CALL stop_program(routineN,moduleN,__LINE__,& 02072 'Cannot have GEMC-NPT simulation with only one box') 02073 CASE("GEMC_NVT") 02074 IF(nboxes .LE. 1) & 02075 CALL stop_program(routineN,moduleN,__LINE__,& 02076 'Cannot have GEMC-NVT simulation with only one box') 02077 CASE("TRADITIONAL") 02078 IF(mc_par%pmswap .GT. 0.0E0_dp)& 02079 CALL stop_program(routineN,moduleN,__LINE__,& 02080 'You cannot do swap moves in a system with only one box') 02081 CASE("VIRIAL") 02082 IF(nchain_total .NE. 2)& 02083 CALL stop_program(routineN,moduleN,__LINE__,& 02084 'You need exactly two molecules in the box to compute the second virial.') 02085 END SELECT 02086 02087 ! can't choose an AVBMC target atom number higher than the number 02088 ! of units in that molecule 02089 DO imol_type=1,nmol_types 02090 IF(mc_par%avbmc_atom(imol_type) .GT. nunits(imol_type)) THEN 02091 CALL stop_program(routineN,moduleN,__LINE__,& 02092 'Cannot have avbmc_atom higher than the number of atoms for that molecule type!') 02093 ENDIF 02094 ENDDO 02095 02096 END SUBROUTINE mc_input_parameters_check 02097 02098 END MODULE mc_types 02099
1.7.3