Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2023 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Perform a molecular dynamics (MD) run using QUICKSTEP
10 : !> \par History
11 : !> - Added support for Langevin regions (2014/02/05, LT)
12 : !> \author Matthias Krack (07.11.2002)
13 : ! **************************************************************************************************
14 : MODULE md_run
15 : USE atomic_kind_list_types, ONLY: atomic_kind_list_type
16 : USE averages_types, ONLY: average_quantities_type
17 : USE barostat_types, ONLY: barostat_type,&
18 : create_barostat_type
19 : USE cell_types, ONLY: cell_type
20 : USE cp_external_control, ONLY: external_control
21 : USE cp_log_handling, ONLY: cp_get_default_logger,&
22 : cp_logger_type
23 : USE cp_output_handling, ONLY: cp_add_iter_level,&
24 : cp_iterate,&
25 : cp_p_file,&
26 : cp_print_key_finished_output,&
27 : cp_print_key_should_output,&
28 : cp_print_key_unit_nr,&
29 : cp_rm_iter_level
30 : USE cp_subsys_types, ONLY: cp_subsys_get,&
31 : cp_subsys_type
32 : USE distribution_1d_types, ONLY: distribution_1d_type
33 : USE force_env_methods, ONLY: force_env_calc_energy_force
34 : USE force_env_types, ONLY: force_env_get,&
35 : force_env_type
36 : USE free_energy_methods, ONLY: free_energy_evaluate
37 : USE free_energy_types, ONLY: fe_env_create,&
38 : free_energy_type
39 : USE global_types, ONLY: global_environment_type
40 : USE input_constants, ONLY: &
41 : ehrenfest, langevin_ensemble, npe_f_ensemble, npe_i_ensemble, &
42 : nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, npt_f_ensemble, npt_i_ensemble, &
43 : npt_ia_ensemble, reftraj_ensemble
44 : USE input_cp2k_check, ONLY: remove_restart_info
45 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
46 : section_vals_remove_values,&
47 : section_vals_type,&
48 : section_vals_val_get
49 : USE kinds, ONLY: default_string_length,&
50 : dp
51 : USE machine, ONLY: m_walltime
52 : USE md_ener_types, ONLY: create_md_ener,&
53 : md_ener_type
54 : USE md_energies, ONLY: initialize_md_ener,&
55 : md_ener_reftraj,&
56 : md_energy,&
57 : md_write_output
58 : USE md_environment_types, ONLY: get_md_env,&
59 : md_env_create,&
60 : md_env_release,&
61 : md_environment_type,&
62 : need_per_atom_wiener_process,&
63 : set_md_env
64 : USE md_util, ONLY: md_output
65 : USE md_vel_utils, ONLY: angvel_control,&
66 : comvel_control,&
67 : setup_velocities,&
68 : temperature_control
69 : USE mdctrl_methods, ONLY: mdctrl_callback
70 : USE mdctrl_types, ONLY: mdctrl_type
71 : USE message_passing, ONLY: mp_para_env_type
72 : USE metadynamics, ONLY: metadyn_finalise_plumed,&
73 : metadyn_forces,&
74 : metadyn_initialise_plumed,&
75 : metadyn_write_colvar
76 : USE metadynamics_types, ONLY: set_meta_env
77 : USE particle_list_types, ONLY: particle_list_type
78 : USE qs_environment_methods, ONLY: qs_env_time_update
79 : USE reftraj_types, ONLY: create_reftraj,&
80 : reftraj_type
81 : USE reftraj_util, ONLY: initialize_reftraj,&
82 : write_output_reftraj
83 : USE rt_propagation, ONLY: rt_prop_setup
84 : USE simpar_methods, ONLY: read_md_section
85 : USE simpar_types, ONLY: create_simpar_type,&
86 : release_simpar_type,&
87 : simpar_type
88 : USE thermal_region_types, ONLY: thermal_regions_type
89 : USE thermal_region_utils, ONLY: create_thermal_regions,&
90 : print_thermal_regions_langevin
91 : USE thermostat_methods, ONLY: create_thermostats
92 : USE thermostat_types, ONLY: thermostats_type
93 : USE velocity_verlet_control, ONLY: velocity_verlet
94 : USE virial_methods, ONLY: virial_evaluate
95 : USE virial_types, ONLY: virial_type
96 : USE wiener_process, ONLY: create_wiener_process,&
97 : create_wiener_process_cv
98 : #include "../base/base_uses.f90"
99 :
100 : IMPLICIT NONE
101 :
102 : PRIVATE
103 :
104 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'md_run'
105 :
106 : PUBLIC :: qs_mol_dyn
107 :
108 : CONTAINS
109 :
110 : ! **************************************************************************************************
111 : !> \brief Main driver module for Molecular Dynamics
112 : !> \param force_env ...
113 : !> \param globenv ...
114 : !> \param averages ...
115 : !> \param rm_restart_info ...
116 : !> \param hmc_e_initial ...
117 : !> \param hmc_e_final ...
118 : !> \param mdctrl ...
119 : ! **************************************************************************************************
120 1727 : SUBROUTINE qs_mol_dyn(force_env, globenv, averages, rm_restart_info, hmc_e_initial, hmc_e_final, mdctrl)
121 :
122 : TYPE(force_env_type), POINTER :: force_env
123 : TYPE(global_environment_type), POINTER :: globenv
124 : TYPE(average_quantities_type), OPTIONAL, POINTER :: averages
125 : LOGICAL, INTENT(IN), OPTIONAL :: rm_restart_info
126 : REAL(KIND=dp), OPTIONAL :: hmc_e_initial, hmc_e_final
127 : TYPE(mdctrl_type), OPTIONAL, POINTER :: mdctrl
128 :
129 : LOGICAL :: my_rm_restart_info
130 : TYPE(md_environment_type), POINTER :: md_env
131 : TYPE(mp_para_env_type), POINTER :: para_env
132 : TYPE(section_vals_type), POINTER :: md_section, motion_section
133 :
134 1727 : my_rm_restart_info = .TRUE.
135 1727 : IF (PRESENT(rm_restart_info)) my_rm_restart_info = rm_restart_info
136 1727 : NULLIFY (md_env, para_env)
137 1727 : para_env => force_env%para_env
138 1727 : motion_section => section_vals_get_subs_vals(force_env%root_section, "MOTION")
139 1727 : md_section => section_vals_get_subs_vals(motion_section, "MD")
140 :
141 : ! Real call to MD driver - Low Level
142 1727 : ALLOCATE (md_env)
143 1727 : CALL md_env_create(md_env, md_section, para_env, force_env=force_env)
144 1727 : CALL set_md_env(md_env, averages=averages)
145 1727 : IF (PRESENT(hmc_e_initial) .AND. PRESENT(hmc_e_final)) THEN
146 : CALL qs_mol_dyn_low(md_env, md_section, motion_section, force_env, globenv, &
147 28 : hmc_e_initial=hmc_e_initial, hmc_e_final=hmc_e_final)
148 : ELSE
149 1699 : CALL qs_mol_dyn_low(md_env, md_section, motion_section, force_env, globenv, mdctrl=mdctrl)
150 : END IF
151 1727 : CALL md_env_release(md_env)
152 1727 : DEALLOCATE (md_env)
153 :
154 : ! Clean restartable sections..
155 1727 : IF (my_rm_restart_info) CALL remove_restart_info(force_env%root_section)
156 1727 : END SUBROUTINE qs_mol_dyn
157 :
158 : ! **************************************************************************************************
159 : !> \brief Purpose: Driver routine for MD run using QUICKSTEP.
160 : !> \param md_env ...
161 : !> \param md_section ...
162 : !> \param motion_section ...
163 : !> \param force_env ...
164 : !> \param globenv ...
165 : !> \param hmc_e_initial ...
166 : !> \param hmc_e_final ...
167 : !> \param mdctrl ...
168 : !> \par History
169 : !> - Cleaning (09.2007) Teodoro Laino [tlaino] - University of Zurich
170 : !> - Added lines to print out langevin regions (2014/02/04, LT)
171 : !> \author Creation (07.11.2002,MK)
172 : ! **************************************************************************************************
173 6908 : SUBROUTINE qs_mol_dyn_low(md_env, md_section, motion_section, force_env, globenv, hmc_e_initial, hmc_e_final, mdctrl)
174 :
175 : TYPE(md_environment_type), POINTER :: md_env
176 : TYPE(section_vals_type), POINTER :: md_section, motion_section
177 : TYPE(force_env_type), POINTER :: force_env
178 : TYPE(global_environment_type), POINTER :: globenv
179 : REAL(KIND=dp), OPTIONAL :: hmc_e_initial, hmc_e_final
180 : TYPE(mdctrl_type), OPTIONAL, POINTER :: mdctrl
181 :
182 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_mol_dyn_low'
183 :
184 : CHARACTER(LEN=default_string_length) :: my_act, my_pos
185 : INTEGER :: handle, i, istep, md_stride, run_type_id
186 : INTEGER, POINTER :: itimes
187 : LOGICAL :: check, ehrenfest_md, save_mem, &
188 : should_stop, write_binary_restart_file
189 : REAL(KIND=dp) :: dummy, time_iter_start, time_iter_stop
190 : REAL(KIND=dp), POINTER :: constant, time, used_time
191 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
192 : TYPE(barostat_type), POINTER :: barostat
193 : TYPE(cell_type), POINTER :: cell
194 : TYPE(cp_logger_type), POINTER :: logger
195 : TYPE(cp_subsys_type), POINTER :: subsys, subsys_i
196 : TYPE(distribution_1d_type), POINTER :: local_particles
197 : TYPE(free_energy_type), POINTER :: fe_env
198 : TYPE(md_ener_type), POINTER :: md_ener
199 : TYPE(mp_para_env_type), POINTER :: para_env
200 : TYPE(particle_list_type), POINTER :: particles
201 : TYPE(reftraj_type), POINTER :: reftraj
202 : TYPE(section_vals_type), POINTER :: constraint_section, force_env_section, &
203 : free_energy_section, global_section, reftraj_section, subsys_section, work_section
204 : TYPE(simpar_type), POINTER :: simpar
205 : TYPE(thermal_regions_type), POINTER :: thermal_regions
206 : TYPE(thermostats_type), POINTER :: thermostats
207 : TYPE(virial_type), POINTER :: virial
208 :
209 1727 : CALL timeset(routineN, handle)
210 1727 : CPASSERT(ASSOCIATED(globenv))
211 1727 : CPASSERT(ASSOCIATED(force_env))
212 :
213 1727 : NULLIFY (particles, cell, simpar, itimes, used_time, subsys, &
214 1727 : md_ener, thermostats, barostat, reftraj, force_env_section, &
215 1727 : reftraj_section, work_section, atomic_kinds, &
216 1727 : local_particles, time, fe_env, free_energy_section, &
217 1727 : constraint_section, thermal_regions, virial, subsys_i)
218 1727 : logger => cp_get_default_logger()
219 1727 : para_env => force_env%para_env
220 :
221 1727 : global_section => section_vals_get_subs_vals(force_env%root_section, "GLOBAL")
222 1727 : free_energy_section => section_vals_get_subs_vals(motion_section, "FREE_ENERGY")
223 1727 : constraint_section => section_vals_get_subs_vals(motion_section, "CONSTRAINT")
224 1727 : CALL section_vals_val_get(global_section, "SAVE_MEM", l_val=save_mem)
225 :
226 1727 : CALL section_vals_val_get(global_section, "RUN_TYPE", i_val=run_type_id)
227 1727 : IF (run_type_id == ehrenfest) CALL set_md_env(md_env, ehrenfest_md=.TRUE.)
228 :
229 1727 : CALL create_simpar_type(simpar)
230 1727 : force_env_section => force_env%force_env_section
231 1727 : subsys_section => section_vals_get_subs_vals(force_env_section, "SUBSYS")
232 1727 : CALL cp_add_iter_level(logger%iter_info, "MD")
233 1727 : CALL cp_iterate(logger%iter_info, iter_nr=0)
234 : ! Read MD section
235 1727 : CALL read_md_section(simpar, motion_section, md_section)
236 : ! Setup print_keys
237 : simpar%info_constraint = cp_print_key_unit_nr(logger, constraint_section, &
238 1727 : "CONSTRAINT_INFO", extension=".shakeLog", log_filename=.FALSE.)
239 : simpar%lagrange_multipliers = cp_print_key_unit_nr(logger, constraint_section, &
240 1727 : "LAGRANGE_MULTIPLIERS", extension=".LagrangeMultLog", log_filename=.FALSE.)
241 : simpar%dump_lm = BTEST(cp_print_key_should_output(logger%iter_info, constraint_section, &
242 1727 : "LAGRANGE_MULTIPLIERS"), cp_p_file)
243 :
244 : ! Create the structure for the md energies
245 1727 : ALLOCATE (md_ener)
246 1727 : CALL create_md_ener(md_ener)
247 1727 : CALL set_md_env(md_env, md_ener=md_ener)
248 1727 : NULLIFY (md_ener)
249 :
250 : ! If requested setup Thermostats
251 : CALL create_thermostats(thermostats, md_section, force_env, simpar, para_env, &
252 1727 : globenv, global_section)
253 :
254 : ! If requested setup Barostat
255 1727 : CALL create_barostat_type(barostat, md_section, force_env, simpar, globenv)
256 :
257 : ! If requested setup different thermal regions
258 1727 : CALL create_thermal_regions(thermal_regions, md_section, simpar, force_env)
259 :
260 : ! If doing langevin_ensemble, then print out langevin_regions information upon request
261 1727 : IF (simpar%ensemble == langevin_ensemble) THEN
262 42 : my_pos = "REWIND"
263 42 : my_act = "WRITE"
264 : CALL print_thermal_regions_langevin(thermal_regions, simpar, &
265 42 : pos=my_pos, act=my_act)
266 : END IF
267 :
268 1727 : CALL set_md_env(md_env, thermostats=thermostats, barostat=barostat, thermal_regions=thermal_regions)
269 :
270 1727 : CALL get_md_env(md_env, ehrenfest_md=ehrenfest_md)
271 :
272 : !If requested set up the REFTRAJ run
273 1727 : IF (simpar%ensemble == reftraj_ensemble .AND. ehrenfest_md) &
274 0 : CPABORT("Ehrenfest MD does not support reftraj ensemble ")
275 1727 : IF (simpar%ensemble == reftraj_ensemble) THEN
276 32 : reftraj_section => section_vals_get_subs_vals(md_section, "REFTRAJ")
277 32 : ALLOCATE (reftraj)
278 32 : CALL create_reftraj(reftraj, reftraj_section, para_env)
279 32 : CALL set_md_env(md_env, reftraj=reftraj)
280 : END IF
281 :
282 : CALL force_env_get(force_env, subsys=subsys, cell=cell, &
283 1727 : force_env_section=force_env_section)
284 1727 : CALL cp_subsys_get(subsys, virial=virial)
285 :
286 : ! Set V0 if needed
287 1727 : IF (simpar%ensemble == nph_uniaxial_ensemble .OR. simpar%ensemble == nph_uniaxial_damped_ensemble) THEN
288 6 : IF (simpar%v0 == 0._dp) simpar%v0 = cell%deth
289 : END IF
290 :
291 : ! Initialize velocities possibly applying constraints at the zeroth MD step
292 : CALL section_vals_val_get(motion_section, "PRINT%RESTART%SPLIT_RESTART_FILE", &
293 1727 : l_val=write_binary_restart_file)
294 : CALL setup_velocities(force_env, simpar, globenv, md_env, md_section, constraint_section, &
295 1727 : write_binary_restart_file)
296 :
297 : ! Setup Free Energy Calculation (if required)
298 1727 : CALL fe_env_create(fe_env, free_energy_section)
299 :
300 : CALL set_md_env(md_env=md_env, simpar=simpar, fe_env=fe_env, cell=cell, &
301 1727 : force_env=force_env)
302 :
303 : ! Possibly initialize Wiener processes
304 : ![NB] Tested again within create_wiener_process. Why??
305 1727 : IF (need_per_atom_wiener_process(md_env)) CALL create_wiener_process(md_env)
306 :
307 1727 : time_iter_start = m_walltime()
308 :
309 : CALL get_md_env(md_env, force_env=force_env, itimes=itimes, constant=constant, &
310 1727 : md_ener=md_ener, t=time, used_time=used_time)
311 :
312 : ! Attach the time counter of the meta_env to the one of the MD
313 1727 : CALL set_meta_env(force_env%meta_env, time=time)
314 :
315 : ! Initialize the md_ener structure
316 1727 : CALL initialize_md_ener(md_ener, force_env, simpar)
317 :
318 : ! Check for ensembles requiring the stress tensor - takes into account the possibility for
319 : ! multiple force_evals
320 : IF ((simpar%ensemble == npt_i_ensemble) .OR. &
321 : (simpar%ensemble == npt_ia_ensemble) .OR. &
322 : (simpar%ensemble == npt_f_ensemble) .OR. &
323 : (simpar%ensemble == npe_f_ensemble) .OR. &
324 : (simpar%ensemble == npe_i_ensemble) .OR. &
325 1727 : (simpar%ensemble == nph_uniaxial_ensemble) .OR. &
326 : (simpar%ensemble == nph_uniaxial_damped_ensemble)) THEN
327 160 : check = virial%pv_availability
328 160 : IF (.NOT. check) &
329 : CALL cp_abort(__LOCATION__, &
330 : "Virial evaluation not requested for this run in the input file! "// &
331 : "You may consider to switch on the virial evaluation with the keyword: STRESS_TENSOR."// &
332 0 : "Be sure the method you are using can compute the virial!")
333 160 : IF (ASSOCIATED(force_env%sub_force_env)) THEN
334 26 : DO i = 1, SIZE(force_env%sub_force_env)
335 26 : IF (ASSOCIATED(force_env%sub_force_env(i)%force_env)) THEN
336 10 : CALL force_env_get(force_env%sub_force_env(i)%force_env, subsys=subsys_i)
337 10 : CALL cp_subsys_get(subsys_i, virial=virial)
338 10 : check = check .AND. virial%pv_availability
339 : END IF
340 : END DO
341 : END IF
342 160 : IF (.NOT. check) &
343 : CALL cp_abort(__LOCATION__, &
344 : "Virial evaluation not requested for all the force_eval sections present in"// &
345 : " the input file! You have to switch on the virial evaluation with the keyword: STRESS_TENSOR "// &
346 0 : " in each force_eval section. Be sure the method you are using can compute the virial!")
347 : END IF
348 :
349 : ! Computing Forces at zero MD step
350 1727 : IF (simpar%ensemble /= reftraj_ensemble) THEN
351 1695 : CALL section_vals_val_get(md_section, "STEP_START_VAL", i_val=itimes)
352 1695 : CALL section_vals_val_get(md_section, "TIME_START_VAL", r_val=time)
353 1695 : CALL section_vals_val_get(md_section, "ECONS_START_VAL", r_val=constant)
354 1695 : CALL cp_iterate(logger%iter_info, iter_nr=itimes)
355 1695 : IF (save_mem) THEN
356 2 : work_section => section_vals_get_subs_vals(subsys_section, "VELOCITY")
357 2 : CALL section_vals_remove_values(work_section)
358 2 : work_section => section_vals_get_subs_vals(subsys_section, "SHELL_VELOCITY")
359 2 : CALL section_vals_remove_values(work_section)
360 2 : work_section => section_vals_get_subs_vals(subsys_section, "CORE_VELOCITY")
361 2 : CALL section_vals_remove_values(work_section)
362 : END IF
363 :
364 1695 : IF (ehrenfest_md) THEN
365 66 : CALL rt_prop_setup(force_env)
366 66 : force_env%qs_env%rtp%dt = simpar%dt
367 : ELSE
368 : ![NB] Lets let all methods, even ones without consistent energies, succeed here.
369 : ! They'll fail in actual integrator if needed
370 : ! consistent_energies=.FALSE. by default
371 1629 : CALL force_env_calc_energy_force(force_env, calc_force=.TRUE.)
372 : END IF
373 :
374 1695 : IF (ASSOCIATED(force_env%qs_env)) THEN
375 584 : CALL qs_env_time_update(force_env%qs_env, time, itimes)
376 : !deb force_env%qs_env%sim_time = time
377 : !deb force_env%qs_env%sim_step = itimes
378 : END IF
379 : ! Warm-up engines for metadynamics
380 1695 : IF (ASSOCIATED(force_env%meta_env)) THEN
381 : ! Setup stuff for plumed if needed
382 144 : IF (force_env%meta_env%use_plumed .EQV. .TRUE.) THEN
383 2 : CALL metadyn_initialise_plumed(force_env, simpar, itimes)
384 : ELSE
385 142 : IF (force_env%meta_env%langevin) THEN
386 4 : CALL create_wiener_process_cv(force_env%meta_env)
387 : END IF
388 142 : IF (force_env%meta_env%well_tempered) THEN
389 2 : force_env%meta_env%wttemperature = simpar%temp_ext
390 2 : IF (force_env%meta_env%wtgamma > EPSILON(1._dp)) THEN
391 0 : dummy = force_env%meta_env%wttemperature*(force_env%meta_env%wtgamma - 1._dp)
392 0 : IF (force_env%meta_env%delta_t > EPSILON(1._dp)) THEN
393 0 : check = ABS(force_env%meta_env%delta_t - dummy) < 1.E+3_dp*EPSILON(1._dp)
394 0 : IF (.NOT. check) &
395 : CALL cp_abort(__LOCATION__, &
396 : "Inconsistency between DELTA_T and WTGAMMA (both specified):"// &
397 0 : " please, verify that DELTA_T=(WTGAMMA-1)*TEMPERATURE")
398 : ELSE
399 0 : force_env%meta_env%delta_t = dummy
400 : END IF
401 : ELSE
402 : force_env%meta_env%wtgamma = 1._dp &
403 2 : + force_env%meta_env%delta_t/force_env%meta_env%wttemperature
404 : END IF
405 2 : force_env%meta_env%invdt = 1._dp/force_env%meta_env%delta_t
406 : END IF
407 142 : CALL metadyn_forces(force_env)
408 142 : CALL metadyn_write_colvar(force_env)
409 : END IF
410 : END IF
411 :
412 1695 : IF (simpar%do_respa) THEN
413 : CALL force_env_calc_energy_force(force_env%sub_force_env(1)%force_env, &
414 6 : calc_force=.TRUE.)
415 : END IF
416 :
417 1695 : CALL force_env_get(force_env, subsys=subsys)
418 :
419 : CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, &
420 1695 : particles=particles, virial=virial)
421 :
422 : CALL virial_evaluate(atomic_kinds%els, particles%els, local_particles, &
423 1695 : virial, force_env%para_env)
424 :
425 1695 : CALL md_energy(md_env, md_ener)
426 1695 : CALL md_write_output(md_env) !inits the print env at itimes == 0 also writes trajectories
427 1695 : md_stride = 1
428 : ELSE
429 32 : CALL get_md_env(md_env, reftraj=reftraj)
430 32 : CALL initialize_reftraj(reftraj, reftraj_section, md_env)
431 32 : itimes = reftraj%info%first_snapshot - 1
432 32 : md_stride = reftraj%info%stride
433 32 : IF (ASSOCIATED(force_env%meta_env)) THEN
434 4 : IF (force_env%meta_env%use_plumed .EQV. .TRUE.) THEN
435 0 : CALL metadyn_initialise_plumed(force_env, simpar, itimes)
436 : END IF
437 : END IF
438 : END IF
439 :
440 : CALL cp_print_key_finished_output(simpar%info_constraint, logger, &
441 1727 : constraint_section, "CONSTRAINT_INFO")
442 : CALL cp_print_key_finished_output(simpar%lagrange_multipliers, logger, &
443 1727 : constraint_section, "LAGRANGE_MULTIPLIERS")
444 :
445 : ! if we need the initial kinetic energy for Hybrid Monte Carlo
446 1727 : IF (PRESENT(hmc_e_initial)) hmc_e_initial = md_ener%ekin
447 :
448 1727 : IF (itimes >= simpar%max_steps) CALL cp_abort(__LOCATION__, &
449 0 : "maximum step number smaller than initial step value")
450 :
451 : ! Real MD Loop
452 41377 : DO istep = 1, simpar%nsteps, md_stride
453 : ! Increase counters
454 39683 : itimes = itimes + 1
455 39683 : time = time + simpar%dt
456 : !needed when electric field fields are applied
457 39683 : IF (ASSOCIATED(force_env%qs_env)) THEN
458 3038 : CALL qs_env_time_update(force_env%qs_env, time, itimes)
459 : !deb force_env%qs_env%sim_time = time
460 : !deb force_env%qs_env%sim_step = itimes
461 : END IF
462 39683 : IF (ehrenfest_md) force_env%qs_env%rtp%istep = istep
463 :
464 39683 : IF (.NOT. logger%iter_info%last_iter(logger%iter_info%n_rlevel)) THEN
465 39681 : CALL cp_iterate(logger%iter_info, last=(istep == simpar%nsteps), iter_nr=itimes)
466 : ELSE
467 2 : CALL cp_iterate(logger%iter_info, last=.TRUE., iter_nr=itimes)
468 : END IF
469 :
470 : ! Open possible Shake output units
471 : simpar%info_constraint = cp_print_key_unit_nr(logger, constraint_section, "CONSTRAINT_INFO", &
472 39683 : extension=".shakeLog", log_filename=.FALSE.)
473 : simpar%lagrange_multipliers = cp_print_key_unit_nr( &
474 : logger, constraint_section, &
475 39683 : "LAGRANGE_MULTIPLIERS", extension=".LagrangeMultLog", log_filename=.FALSE.)
476 : simpar%dump_lm = BTEST(cp_print_key_should_output(logger%iter_info, constraint_section, &
477 39683 : "LAGRANGE_MULTIPLIERS"), cp_p_file)
478 :
479 : ! Velocity Verlet Integrator
480 39683 : CALL velocity_verlet(md_env, globenv)
481 :
482 : ! Close Shake output if requested...
483 : CALL cp_print_key_finished_output(simpar%info_constraint, logger, &
484 39683 : constraint_section, "CONSTRAINT_INFO")
485 : CALL cp_print_key_finished_output(simpar%lagrange_multipliers, logger, &
486 39683 : constraint_section, "LAGRANGE_MULTIPLIERS")
487 :
488 : ! Free Energy calculation
489 39683 : CALL free_energy_evaluate(md_env, should_stop, free_energy_section)
490 :
491 39683 : IF (should_stop) EXIT
492 :
493 : ! Test for <PROJECT_NAME>.EXIT_MD or for WALL_TIME to exit
494 : ! Default:
495 : ! IF so we don't overwrite the restart or append to the trajectory
496 : ! because the execution could in principle stop inside the SCF where energy
497 : ! and forces are not converged.
498 : ! But:
499 : ! You can force to print the last step (for example if the method used
500 : ! to compute energy and forces is not SCF based) activating the print_key
501 : ! MOTION%MD%PRINT%FORCE_LAST.
502 39683 : CALL external_control(should_stop, "MD", globenv=globenv)
503 :
504 : !check if upper bound of total steps has been reached
505 39683 : IF (.NOT. (istep == simpar%nsteps) .AND. logger%iter_info%last_iter(logger%iter_info%n_rlevel)) should_stop = .TRUE.
506 39683 : IF (itimes >= simpar%max_steps) should_stop = .TRUE.
507 :
508 : ! call external hook e.g. from global optimization
509 39683 : IF (PRESENT(mdctrl)) &
510 2633 : CALL mdctrl_callback(mdctrl, md_env, should_stop)
511 :
512 39683 : IF (should_stop) THEN
513 33 : CALL cp_iterate(logger%iter_info, last=.TRUE., iter_nr=itimes)
514 : !In Ehrenfest molecular dynamics the external control is only checked after a converged propagation
515 : !The restart needs to be written in order to be consistent with the mos/density matrix for the restart
516 33 : IF (run_type_id == ehrenfest) THEN
517 2 : CALL md_output(md_env, md_section, force_env%root_section, .FALSE.)
518 : ELSE
519 31 : CALL md_output(md_env, md_section, force_env%root_section, should_stop)
520 : END IF
521 : EXIT
522 : END IF
523 :
524 39650 : IF (simpar%ensemble /= reftraj_ensemble) THEN
525 39388 : CALL md_energy(md_env, md_ener)
526 39388 : CALL temperature_control(simpar, md_env, md_ener, force_env, logger)
527 39388 : CALL comvel_control(md_ener, force_env, md_section, logger)
528 39388 : CALL angvel_control(md_ener, force_env, md_section, logger)
529 : ELSE
530 262 : CALL md_ener_reftraj(md_env, md_ener)
531 : END IF
532 :
533 39650 : time_iter_stop = m_walltime()
534 39650 : used_time = time_iter_stop - time_iter_start
535 39650 : time_iter_start = time_iter_stop
536 :
537 39650 : CALL md_output(md_env, md_section, force_env%root_section, should_stop)
538 120710 : IF (simpar%ensemble == reftraj_ensemble) THEN
539 262 : CALL write_output_reftraj(md_env)
540 : END IF
541 : END DO
542 :
543 : ! if we need the final kinetic energy for Hybrid Monte Carlo
544 1727 : IF (PRESENT(hmc_e_final)) hmc_e_final = md_ener%ekin
545 :
546 : ! Remove the iteration level
547 1727 : CALL cp_rm_iter_level(logger%iter_info, "MD")
548 :
549 : ! Clean up PLUMED
550 1727 : IF (ASSOCIATED(force_env%meta_env)) THEN
551 148 : IF (force_env%meta_env%use_plumed .EQV. .TRUE.) THEN
552 2 : CALL metadyn_finalise_plumed()
553 : END IF
554 : END IF
555 :
556 : ! Deallocate Thermostats and Barostats
557 1727 : CALL release_simpar_type(simpar)
558 1727 : CALL timestop(handle)
559 :
560 1727 : END SUBROUTINE qs_mol_dyn_low
561 :
562 : END MODULE md_run
|
