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