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