1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Methods for storing MD parameters type 8!> \author CJM 9!> \author Teodoro Laino [tlaino] - University of Zurich - 10.2008 10!> reorganization of the original routines/modules 11! ************************************************************************************************** 12MODULE simpar_methods 13 USE bibliography, ONLY: Evans1983,& 14 Kuhne2007,& 15 Minary2003,& 16 Ricci2003,& 17 cite_reference 18 USE cp_log_handling, ONLY: cp_get_default_logger,& 19 cp_logger_type 20 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 21 cp_print_key_generate_filename,& 22 cp_print_key_unit_nr 23 USE cp_units, ONLY: cp_unit_from_cp2k 24 USE input_constants, ONLY: & 25 isokin_ensemble, langevin_ensemble, npe_f_ensemble, npe_i_ensemble, & 26 nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, npt_f_ensemble, npt_i_ensemble, & 27 nvt_ensemble, reftraj_ensemble 28 USE input_cp2k_md, ONLY: create_md_section 29 USE input_enumeration_types, ONLY: enum_i2c,& 30 enumeration_type 31 USE input_keyword_types, ONLY: keyword_get,& 32 keyword_type 33 USE input_section_types, ONLY: section_get_keyword,& 34 section_release,& 35 section_type,& 36 section_vals_get,& 37 section_vals_get_subs_vals,& 38 section_vals_type,& 39 section_vals_val_get 40 USE kinds, ONLY: default_path_length,& 41 dp 42 USE simpar_types, ONLY: simpar_type 43#include "../base/base_uses.f90" 44 45 IMPLICIT NONE 46 47 PRIVATE 48 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'simpar_methods' 49 PUBLIC :: read_md_section 50 51CONTAINS 52 53! ************************************************************************************************** 54!> \brief Reads the MD section and setup the simulation parameters type 55!> \param simpar ... 56!> \param motion_section ... 57!> \param md_section ... 58!> \author Teodoro Laino 59! ************************************************************************************************** 60 SUBROUTINE read_md_section(simpar, motion_section, md_section) 61 TYPE(simpar_type), POINTER :: simpar 62 TYPE(section_vals_type), POINTER :: motion_section, md_section 63 64 CHARACTER(len=*), PARAMETER :: routineN = 'read_md_section', & 65 routineP = moduleN//':'//routineN 66 67 CHARACTER(LEN=default_path_length) :: filename 68 INTEGER :: iprint, iw 69 REAL(kind=dp) :: tmp_r1, tmp_r2, tmp_r3 70 TYPE(cp_logger_type), POINTER :: logger 71 TYPE(enumeration_type), POINTER :: enum 72 TYPE(keyword_type), POINTER :: keyword 73 TYPE(section_type), POINTER :: section 74 TYPE(section_vals_type), POINTER :: print_key 75 76 NULLIFY (logger, print_key, enum, keyword, section) 77 logger => cp_get_default_logger() 78 iw = cp_print_key_unit_nr(logger, md_section, "PRINT%PROGRAM_RUN_INFO", extension=".log") 79 80 CALL read_md_low(simpar, motion_section, md_section) 81 IF (iw > 0) WRITE (iw, *) 82 83 ! Begin setup Langevin dynamics 84 IF (simpar%ensemble == langevin_ensemble) THEN 85 CALL cite_reference(Ricci2003) 86 IF (simpar%noisy_gamma > 0.0_dp) CALL cite_reference(Kuhne2007) 87 ! Normalization factor using a normal Gaussian random number distribution 88 simpar%var_w = 2.0_dp*simpar%temp_ext*simpar%dt*(simpar%gamma + simpar%noisy_gamma) 89 IF (iw > 0) THEN 90 tmp_r1 = cp_unit_from_cp2k(simpar%gamma, "fs^-1") 91 tmp_r2 = cp_unit_from_cp2k(simpar%noisy_gamma, "fs^-1") 92 tmp_r3 = cp_unit_from_cp2k(simpar%shadow_gamma, "fs^-1") 93 WRITE (UNIT=iw, FMT="(T2,A,T71,ES10.3)") & 94 "LD| Gamma [1/fs] ", tmp_r1, & 95 "LD| Noisy Gamma [1/fs]", tmp_r2, & 96 "LD| Shadow Gamma [1/fs]", tmp_r3, & 97 "LD| Variance [a.u.]", simpar%var_w 98 END IF 99 END IF 100 101 ! Create section for output enumeration infos 102 CALL create_md_section(section) 103 keyword => section_get_keyword(section, "ENSEMBLE") 104 CALL keyword_get(keyword, enum=enum) 105 106 !..write some information to output 107 IF (iw > 0) THEN 108 WRITE (iw, '( A )') ' MD| Molecular Dynamics Protocol ' 109 WRITE (iw, '( A,T61,A20)') ' MD| Ensemble Type ', & 110 ADJUSTR(TRIM(enum_i2c(enum, simpar%ensemble))) 111 WRITE (iw, '( A,T71,I10 )') ' MD| Number of Time Steps ', & 112 simpar%nsteps 113 IF (simpar%variable_dt) THEN 114 WRITE (iw, '( A )') ' MD| Variable Time Step is activated ' 115 tmp_r1 = cp_unit_from_cp2k(simpar%dt, "fs") 116 WRITE (iw, '( A,A2,A,T71,F10.2 )') ' MD| Max. Time Step [', 'fs', '] ', tmp_r1 117 tmp_r1 = cp_unit_from_cp2k(simpar%dr_tol, "angstrom") 118 WRITE (iw, '( A, T71, F10.4 )') ' MD| Max. atomic displacement permitted [A] ', tmp_r1 119 ELSE 120 tmp_r1 = cp_unit_from_cp2k(simpar%dt, "fs") 121 WRITE (iw, '( A,A2,A,T71,F10.2 )') ' MD| Time Step [', 'fs', '] ', tmp_r1 122 END IF 123 tmp_r1 = cp_unit_from_cp2k(simpar%temp_ext, "K") 124 WRITE (iw, '( A,T71,F10.2 )') ' MD| Temperature [K] ', tmp_r1 125 tmp_r1 = cp_unit_from_cp2k(simpar%temp_tol, "K") 126 WRITE (iw, '( A,T71,F10.2 )') ' MD| Temperature tolerance [K] ', tmp_r1 127 128 IF (simpar%annealing) & 129 WRITE (iw, '( A,T71,F10.2 )') ' MD| Annealing ion factor ', & 130 simpar%f_annealing 131 IF ((simpar%ensemble == npe_f_ensemble .OR. simpar%ensemble == npe_i_ensemble) & 132 .AND. simpar%annealing_cell) & 133 WRITE (iw, '( A,T71,F10.2 )') ' MD| Annealing cell factor ', & 134 simpar%f_annealing_cell 135 IF (simpar%ensemble == npt_i_ensemble .OR. & 136 simpar%ensemble == npt_f_ensemble .OR. & 137 simpar%ensemble == npe_i_ensemble .OR. & 138 simpar%ensemble == npe_f_ensemble) THEN 139 tmp_r1 = cp_unit_from_cp2k(simpar%p_ext, "bar") 140 WRITE (iw, '( A,A3,A, T71, F10.2 )') & 141 ' MD| Pressure [', 'Bar', '] ', tmp_r1 142 tmp_r1 = cp_unit_from_cp2k(simpar%tau_cell, "fs") 143 WRITE (iw, '( A,A4,A, T71, F10.2 )') & 144 ' MD| Barostat time constant [', 'fs', '] ', tmp_r1 145 END IF 146 IF (simpar%ensemble == isokin_ensemble) THEN 147 CALL cite_reference(Evans1983) 148 CALL cite_reference(Minary2003) 149 WRITE (iw, '( A )') ' MD| Simulation in the isokinetic ensemble' 150 END IF 151 IF (simpar%constraint) THEN 152 WRITE (iw, '( A )') ' MD| Constraints activated ' 153 WRITE (iw, '( A,T71,G10.4 )') ' MD| Tolerance for shake ', & 154 simpar%shake_tol 155 END IF 156 157 print_key => section_vals_get_subs_vals(motion_section, "MD%PRINT%PROGRAM_RUN_INFO") 158 CALL section_vals_val_get(print_key, "EACH%MD", i_val=iprint) 159 WRITE (iw, '( A,T63,i10,A )') ' MD| Print MD information every', iprint, ' step(s)' 160 WRITE (iw, '( A,T20,A,T71,A10 )') ' MD| File type', 'Print frequency[steps]', 'File names' 161 162 print_key => section_vals_get_subs_vals(motion_section, "PRINT%TRAJECTORY") 163 CALL section_vals_val_get(print_key, "EACH%MD", i_val=iprint) 164 filename = cp_print_key_generate_filename(logger, print_key, & 165 extension=".xyz", middle_name="pos", my_local=.FALSE.) 166 WRITE (iw, '( A,T20,i10,T31,A50 )') ' MD| Coordinates', iprint, & 167 ADJUSTR(TRIM(filename)) 168 169 IF (simpar%ensemble == nph_uniaxial_ensemble .OR. & 170 simpar%ensemble == nph_uniaxial_damped_ensemble .OR. & 171 simpar%ensemble == npt_i_ensemble .OR. & 172 simpar%ensemble == npt_f_ensemble .OR. & 173 simpar%ensemble == npe_i_ensemble .OR. & 174 simpar%ensemble == npe_f_ensemble) THEN 175 176 print_key => section_vals_get_subs_vals(motion_section, "PRINT%CELL") 177 CALL section_vals_val_get(print_key, "EACH%MD", i_val=iprint) 178 filename = cp_print_key_generate_filename(logger, print_key, & 179 extension=".cell", my_local=.FALSE.) 180 WRITE (iw, '( A,T20,i10,T31,A50 )') ' MD| Simulation Cell', iprint, & 181 ADJUSTR(TRIM(filename)) 182 END IF 183 184 print_key => section_vals_get_subs_vals(motion_section, "PRINT%VELOCITIES") 185 CALL section_vals_val_get(print_key, "EACH%MD", i_val=iprint) 186 filename = cp_print_key_generate_filename(logger, print_key, & 187 extension=".xyz", middle_name="vel", my_local=.FALSE.) 188 WRITE (iw, '( A,T20,i10,T31,A50 )') ' MD| Velocities', iprint, & 189 ADJUSTR(TRIM(filename)) 190 191 print_key => section_vals_get_subs_vals(motion_section, "MD%PRINT%ENERGY") 192 CALL section_vals_val_get(print_key, "EACH%MD", i_val=iprint) 193 filename = cp_print_key_generate_filename(logger, print_key, & 194 extension=".ener", my_local=.FALSE.) 195 WRITE (iw, '( A,T20,i10,T31,A50 )') ' MD| Energies', iprint, & 196 ADJUSTR(TRIM(filename)) 197 198 print_key => section_vals_get_subs_vals(motion_section, "PRINT%RESTART") 199 CALL section_vals_val_get(print_key, "EACH%MD", i_val=iprint) 200 filename = cp_print_key_generate_filename(logger, print_key, & 201 extension=".restart", my_local=.FALSE.) 202 WRITE (iw, '( A,T20,i10,T31,A50 )') ' MD| Dump', iprint, & 203 ADJUSTR(TRIM(filename)) 204 205 WRITE (iw, *) 206 IF (simpar%ensemble == nph_uniaxial_ensemble .OR. & 207 simpar%ensemble == nph_uniaxial_damped_ensemble) THEN 208 WRITE (iw, '( A )') ' SHOCK| Uniaxial Shock Parameters: ' 209 tmp_r1 = cp_unit_from_cp2k(simpar%v_shock, "m*s^-1") 210 WRITE (iw, '( A,A4,A,T71,F10.4 )') & 211 ' SHOCK| Shock Velocity [', 'm/s', '] ', tmp_r1 212 tmp_r1 = cp_unit_from_cp2k(simpar%gamma_nph, "fs^-1") 213 WRITE (iw, '( A,A4,A,T71,F10.4 )') & 214 ' SHOCK| Damping Coefficient [', '1/fs', '] ', tmp_r1 215 tmp_r1 = cp_unit_from_cp2k(simpar%p0, "bar") 216 WRITE (iw, '( A,A3,A, T71, F10.2 )') & 217 ' SHOCK| Pressure [', 'Bar', '] ', tmp_r1 218 WRITE (iw, '( A,A4,A, T71, E10.4 )') & 219 ' SHOCK| Barostat Mass [', 'au', '] ', simpar%cmass 220 END IF 221 ! Print warning for temp_tol 222 IF (simpar%temp_tol > 0.0_dp) THEN 223 CALL cp_warn(__LOCATION__, & 224 "A temperature tolerance (TEMP_TOL) is used during the MD. "// & 225 "Due to the velocity rescaling algorithm jumps may appear in the conserved quantity.") 226 END IF 227 ! Print warning for annealing 228 IF (simpar%annealing) THEN 229 IF (simpar%ensemble == nvt_ensemble .OR. & 230 simpar%ensemble == npt_i_ensemble .OR. & 231 simpar%ensemble == npt_f_ensemble) THEN 232 CALL cp_abort(__LOCATION__, & 233 "Annealing of the ions has been required "// & 234 "even if the thermostat is active (nvt or npt_i or npt_f) "// & 235 "These two methods to control the temperature act one against the other.") 236 END IF 237 END IF 238 ! Print warning for variable time step 239 IF (simpar%variable_dt) THEN 240 IF (simpar%ensemble == langevin_ensemble .OR. & 241 simpar%ensemble == reftraj_ensemble .OR. & 242 simpar%do_respa) THEN 243 CALL cp_warn( & 244 __LOCATION__, & 245 "The variable timestep has been required, however "// & 246 "this option is not available either with the Langevin ensemble or with the multiple timestep schme. "// & 247 "The run will proceed with constant timestep, as read from input.") 248 END IF 249 END IF 250 END IF 251 CALL section_release(section) 252 CALL cp_print_key_finished_output(iw, logger, md_section, & 253 "PRINT%PROGRAM_RUN_INFO") 254 255 END SUBROUTINE read_md_section 256 257! ************************************************************************************************** 258!> \brief Low Level: Parses the MD input section 259!> \param simpar ... 260!> \param motion_section ... 261!> \param md_section ... 262!> \author teo 263! ************************************************************************************************** 264 SUBROUTINE read_md_low(simpar, motion_section, md_section) 265 TYPE(simpar_type), POINTER :: simpar 266 TYPE(section_vals_type), POINTER :: motion_section, md_section 267 268 CHARACTER(len=*), PARAMETER :: routineN = 'read_md_low', routineP = moduleN//':'//routineN 269 270 LOGICAL :: explicit 271 TYPE(section_vals_type), POINTER :: tmp_section 272 273 NULLIFY (tmp_section) 274 CALL section_vals_val_get(md_section, "ENSEMBLE", i_val=simpar%ensemble) 275 CALL section_vals_val_get(md_section, "STEPS", i_val=simpar%nsteps) 276 CALL section_vals_val_get(md_section, "MAX_STEPS", i_val=simpar%max_steps) 277 CALL section_vals_val_get(md_section, "TEMPERATURE", r_val=simpar%temp_ext) 278 CALL section_vals_val_get(md_section, "TEMP_TOL", r_val=simpar%temp_tol) 279 CALL section_vals_val_get(md_section, "ANGVEL_ZERO", l_val=simpar%angvel_zero) 280 CALL section_vals_val_get(md_section, "TEMP_KIND", l_val=simpar%temperature_per_kind) 281 CALL section_vals_val_get(md_section, "SCALE_TEMP_KIND", l_val=simpar%scale_temperature_per_kind) 282 CALL section_vals_val_get(md_section, "ANNEALING", r_val=simpar%f_annealing, explicit=simpar%annealing) 283 CALL section_vals_val_get(md_section, "ANNEALING_CELL", r_val=simpar%f_annealing_cell, & 284 explicit=simpar%annealing_cell) 285 CALL section_vals_val_get(md_section, "TEMPERATURE_ANNEALING", r_val=simpar%f_temperature_annealing, & 286 explicit=simpar%temperature_annealing) 287 CALL section_vals_val_get(md_section, "DISPLACEMENT_TOL", r_val=simpar%dr_tol, & 288 explicit=simpar%variable_dt) 289 CALL section_vals_val_get(md_section, "TIMESTEP", r_val=simpar%dt) 290 CALL section_vals_val_get(md_section, "INITIALIZATION_METHOD", & 291 i_val=simpar%initialization_method) 292 ! Initialize dt_fact to 1.0 293 simpar%dt_fact = 1.0_dp 294 295 IF (simpar%ensemble == langevin_ensemble) THEN 296 CALL section_vals_val_get(md_section, "LANGEVIN%GAMMA", r_val=simpar%gamma) 297 CALL section_vals_val_get(md_section, "LANGEVIN%NOISY_GAMMA", r_val=simpar%noisy_gamma) 298 CALL section_vals_val_get(md_section, "LANGEVIN%SHADOW_GAMMA", r_val=simpar%shadow_gamma) 299 END IF 300 301 tmp_section => section_vals_get_subs_vals(motion_section, "CONSTRAINT") 302 CALL section_vals_get(tmp_section, explicit=simpar%constraint) 303 IF (simpar%constraint) THEN 304 CALL section_vals_val_get(tmp_section, "SHAKE_TOLERANCE", r_val=simpar%shake_tol) 305 IF (simpar%shake_tol <= EPSILON(0.0_dp)*1000.0_dp) & 306 CALL cp_warn(__LOCATION__, & 307 "Shake tolerance lower than 1000*EPSILON, where EPSILON is the machine precision. "// & 308 "This may lead to numerical problems. Setting up shake_tol to 1000*EPSILON!") 309 simpar%shake_tol = MAX(EPSILON(0.0_dp)*1000.0_dp, simpar%shake_tol) 310 311 CALL section_vals_val_get(tmp_section, "ROLL_TOLERANCE", r_val=simpar%roll_tol) 312 IF (simpar%roll_tol <= EPSILON(0.0_dp)*1000.0_dp) & 313 CALL cp_warn(__LOCATION__, & 314 "Roll tolerance lower than 1000*EPSILON, where EPSILON is the machine precision. "// & 315 "This may lead to numerical problems. Setting up roll_tol to 1000*EPSILON!") 316 simpar%roll_tol = MAX(EPSILON(0.0_dp)*1000.0_dp, simpar%roll_tol) 317 END IF 318 319 IF (simpar%ensemble == nph_uniaxial_ensemble .OR. simpar%ensemble == nph_uniaxial_damped_ensemble) THEN 320 tmp_section => section_vals_get_subs_vals(md_section, "MSST") 321 CALL section_vals_val_get(tmp_section, "PRESSURE", r_val=simpar%p0) 322 CALL section_vals_val_get(tmp_section, "ENERGY", r_val=simpar%e0) 323 CALL section_vals_val_get(tmp_section, "VOLUME", r_val=simpar%v0) 324 CALL section_vals_val_get(tmp_section, "GAMMA", r_val=simpar%gamma_nph) 325 IF (simpar%gamma_nph /= 0.0_dp) simpar%ensemble = nph_uniaxial_damped_ensemble 326 CALL section_vals_val_get(tmp_section, "CMASS", r_val=simpar%cmass) 327 CALL section_vals_val_get(tmp_section, "VSHOCK", r_val=simpar%v_shock) 328 END IF 329 330 SELECT CASE (simpar%ensemble) 331 CASE (nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, & 332 npt_f_ensemble, npt_i_ensemble, npe_f_ensemble, npe_i_ensemble) 333 tmp_section => section_vals_get_subs_vals(md_section, "BAROSTAT") 334 CALL section_vals_val_get(tmp_section, "PRESSURE", r_val=simpar%p_ext) 335 CALL section_vals_val_get(tmp_section, "TIMECON", r_val=simpar%tau_cell) 336 END SELECT 337 338 ! RESPA 339 tmp_section => section_vals_get_subs_vals(md_section, "RESPA") 340 CALL section_vals_get(tmp_section, explicit=simpar%do_respa) 341 CALL section_vals_val_get(tmp_section, "FREQUENCY", i_val=simpar%n_time_steps) 342 simpar%multi_time_switch = simpar%do_respa 343 344 ! CORE-SHELL MODEL 345 tmp_section => section_vals_get_subs_vals(md_section, "SHELL") 346 CALL section_vals_val_get(tmp_section, "TEMPERATURE", r_val=simpar%temp_sh_ext) 347 CALL section_vals_val_get(tmp_section, "TEMP_TOL", r_val=simpar%temp_sh_tol) 348 349 CALL section_vals_val_get(tmp_section, "DISPLACEMENT_SHELL_TOL", r_val=simpar%dsc_tol, & 350 explicit=explicit) 351 simpar%variable_dt = simpar%variable_dt .OR. explicit 352 ! ADIABATIC DYNAMICS 353 tmp_section => section_vals_get_subs_vals(md_section, "ADIABATIC_DYNAMICS") 354 CALL section_vals_val_get(tmp_section, "TEMP_FAST", r_val=simpar%temp_fast) 355 CALL section_vals_val_get(tmp_section, "TEMP_SLOW", r_val=simpar%temp_slow) 356 CALL section_vals_val_get(tmp_section, "TEMP_TOL_FAST", r_val=simpar%temp_tol_fast) 357 CALL section_vals_val_get(tmp_section, "TEMP_TOL_SLOW", r_val=simpar%temp_tol_slow) 358 CALL section_vals_val_get(tmp_section, "N_RESP_FAST", i_val=simpar%n_resp_fast) 359 360 ! VELOCITY SOFTENING 361 tmp_section => section_vals_get_subs_vals(md_section, "VELOCITY_SOFTENING") 362 CALL section_vals_val_get(tmp_section, "STEPS", i_val=simpar%soften_nsteps) 363 CALL section_vals_val_get(tmp_section, "ALPHA", r_val=simpar%soften_alpha) 364 CALL section_vals_val_get(tmp_section, "DELTA", r_val=simpar%soften_delta) 365 END SUBROUTINE read_md_low 366 367END MODULE simpar_methods 368