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 temperature accelarated hybrid monte carlo (TAHMC) run using QUICKSTEP
10 : !> \par History
11 : !> none
12 : !> \author Alin M Elena
13 : ! **************************************************************************************************
14 : MODULE tamc_run
15 :
16 : USE atomic_kind_list_types, ONLY: atomic_kind_list_type
17 : USE atomic_kind_types, ONLY: atomic_kind_type
18 : USE averages_types, ONLY: average_quantities_type
19 : USE barostat_types, ONLY: barostat_type,&
20 : create_barostat_type
21 : USE bibliography, ONLY: VandenCic2006
22 : USE cell_types, ONLY: cell_type
23 : USE colvar_methods, ONLY: colvar_eval_glob_f
24 : USE colvar_types, ONLY: HBP_colvar_id,&
25 : WC_colvar_id,&
26 : colvar_p_type
27 : USE constraint_fxd, ONLY: fix_atom_control
28 : USE cp_external_control, ONLY: external_control
29 : USE cp_log_handling, ONLY: cp_get_default_logger,&
30 : cp_logger_get_default_io_unit,&
31 : cp_logger_type
32 : USE cp_output_handling, ONLY: cp_add_iter_level,&
33 : cp_iterate,&
34 : cp_p_file,&
35 : cp_print_key_finished_output,&
36 : cp_print_key_should_output,&
37 : cp_print_key_unit_nr,&
38 : cp_rm_iter_level
39 : USE cp_subsys_types, ONLY: cp_subsys_get,&
40 : cp_subsys_type
41 : USE cp_units, ONLY: cp_unit_from_cp2k
42 : USE distribution_1d_types, ONLY: distribution_1d_type
43 : USE force_env_methods, ONLY: force_env_calc_energy_force
44 : USE force_env_types, ONLY: force_env_get,&
45 : force_env_type
46 : USE free_energy_types, ONLY: fe_env_create,&
47 : free_energy_type
48 : USE global_types, ONLY: global_environment_type
49 : USE input_constants, ONLY: &
50 : langevin_ensemble, npe_f_ensemble, npe_i_ensemble, nph_uniaxial_damped_ensemble, &
51 : nph_uniaxial_ensemble, npt_f_ensemble, npt_i_ensemble, npt_ia_ensemble, reftraj_ensemble
52 : USE input_cp2k_check, ONLY: remove_restart_info
53 : USE input_cp2k_restarts, ONLY: write_restart
54 : USE input_section_types, ONLY: section_vals_get,&
55 : section_vals_get_subs_vals,&
56 : section_vals_remove_values,&
57 : section_vals_type,&
58 : section_vals_val_get,&
59 : section_vals_val_set
60 : USE kinds, ONLY: dp
61 : USE machine, ONLY: m_walltime
62 : USE mc_environment_types, ONLY: get_mc_env,&
63 : mc_env_create,&
64 : mc_env_release,&
65 : mc_environment_type,&
66 : set_mc_env
67 : USE mc_misc, ONLY: mc_averages_create,&
68 : mc_averages_release
69 : USE mc_move_control, ONLY: init_mc_moves,&
70 : mc_moves_release
71 : USE mc_types, ONLY: get_mc_par,&
72 : mc_averages_type,&
73 : mc_ekin_type,&
74 : mc_moves_type,&
75 : mc_simpar_type,&
76 : set_mc_par
77 : USE md_ener_types, ONLY: create_md_ener,&
78 : md_ener_type
79 : USE md_energies, ONLY: initialize_md_ener,&
80 : md_energy
81 : USE md_environment_types, ONLY: get_md_env,&
82 : md_env_create,&
83 : md_env_release,&
84 : md_environment_type,&
85 : set_md_env
86 : USE md_run, ONLY: qs_mol_dyn
87 : USE message_passing, ONLY: mp_comm_type,&
88 : mp_para_env_type
89 : USE metadynamics_types, ONLY: meta_env_type,&
90 : metavar_type,&
91 : set_meta_env
92 : USE molecule_kind_list_types, ONLY: molecule_kind_list_type
93 : USE molecule_kind_types, ONLY: molecule_kind_type
94 : USE molecule_list_types, ONLY: molecule_list_type
95 : USE molecule_types, ONLY: global_constraint_type,&
96 : molecule_type
97 : USE parallel_rng_types, ONLY: UNIFORM,&
98 : rng_stream_type
99 : USE particle_list_types, ONLY: particle_list_type
100 : USE particle_types, ONLY: particle_type
101 : USE physcon, ONLY: boltzmann,&
102 : femtoseconds,&
103 : joule,&
104 : kelvin
105 : USE qmmm_util, ONLY: apply_qmmm_walls_reflective
106 : USE qs_environment_types, ONLY: get_qs_env
107 : USE qs_scf_post_gpw, ONLY: scf_post_calculation_gpw
108 : USE reference_manager, ONLY: cite_reference
109 : USE reftraj_types, ONLY: create_reftraj,&
110 : reftraj_type
111 : USE reftraj_util, ONLY: initialize_reftraj
112 : USE simpar_methods, ONLY: read_md_section
113 : USE simpar_types, ONLY: create_simpar_type,&
114 : release_simpar_type,&
115 : simpar_type
116 : USE string_utilities, ONLY: str_comp
117 : USE thermal_region_types, ONLY: thermal_regions_type
118 : USE thermal_region_utils, ONLY: create_thermal_regions
119 : USE thermostat_methods, ONLY: create_thermostats
120 : USE thermostat_types, ONLY: thermostats_type
121 : USE virial_methods, ONLY: virial_evaluate
122 : USE virial_types, ONLY: virial_type
123 : USE wiener_process, ONLY: create_wiener_process,&
124 : create_wiener_process_cv
125 : !!!!! monte carlo part
126 : #include "../../base/base_uses.f90"
127 :
128 : IMPLICIT NONE
129 :
130 : PRIVATE
131 :
132 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'tamc_run'
133 :
134 : PUBLIC :: qs_tamc
135 :
136 : CONTAINS
137 :
138 : ! **************************************************************************************************
139 : !> \brief Driver routine for TAHMC
140 : !> \param force_env ...
141 : !> \param globenv ...
142 : !> \param averages ...
143 : !> \author Alin M Elena
144 : !> \note it computes the forces using QuickStep.
145 : ! **************************************************************************************************
146 2 : SUBROUTINE qs_tamc(force_env, globenv, averages)
147 :
148 : TYPE(force_env_type), POINTER :: force_env
149 : TYPE(global_environment_type), POINTER :: globenv
150 : TYPE(average_quantities_type), OPTIONAL, POINTER :: averages
151 :
152 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_tamc'
153 :
154 : CHARACTER(LEN=20) :: ensemble
155 : INTEGER :: handle, i, initialStep, iprint, isos, &
156 : istep, j, md_stride, nmccycles, &
157 : output_unit, rand2skip, run_type_id
158 : INTEGER, POINTER :: itimes
159 : LOGICAL :: check, explicit, my_rm_restart_info, &
160 : save_mem, should_stop
161 : REAL(KIND=dp) :: auxRandom, inittime, rval, temp, &
162 : time_iter_start, time_iter_stop
163 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: An, fz, xieta, zbuff
164 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: r
165 : REAL(KIND=dp), POINTER :: constant, time, used_time
166 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
167 : TYPE(barostat_type), POINTER :: barostat
168 : TYPE(cell_type), POINTER :: cell
169 : TYPE(cp_logger_type), POINTER :: logger
170 : TYPE(cp_subsys_type), POINTER :: subsys, subsys_i
171 : TYPE(distribution_1d_type), POINTER :: local_particles
172 : TYPE(free_energy_type), POINTER :: fe_env
173 : TYPE(mc_averages_type), POINTER :: MCaverages
174 : TYPE(mc_environment_type), POINTER :: mc_env
175 : TYPE(mc_moves_type), POINTER :: gmoves, moves
176 : TYPE(mc_simpar_type), POINTER :: mc_par
177 : TYPE(md_ener_type), POINTER :: md_ener
178 : TYPE(md_environment_type), POINTER :: md_env
179 : TYPE(meta_env_type), POINTER :: meta_env_saved
180 : TYPE(mp_para_env_type), POINTER :: para_env
181 : TYPE(particle_list_type), POINTER :: particles
182 : TYPE(reftraj_type), POINTER :: reftraj
183 : TYPE(rng_stream_type) :: rng_stream_mc
184 : TYPE(section_vals_type), POINTER :: constraint_section, force_env_section, &
185 : free_energy_section, fs_section, global_section, mc_section, md_section, motion_section, &
186 : reftraj_section, subsys_section, work_section
187 : TYPE(simpar_type), POINTER :: simpar
188 : TYPE(thermal_regions_type), POINTER :: thermal_regions
189 : TYPE(thermostats_type), POINTER :: thermostats
190 : TYPE(virial_type), POINTER :: virial
191 :
192 2 : initialStep = 0
193 2 : inittime = 0.0_dp
194 :
195 2 : CALL timeset(routineN, handle)
196 2 : my_rm_restart_info = .TRUE.
197 2 : NULLIFY (para_env, fs_section, virial)
198 2 : para_env => force_env%para_env
199 2 : motion_section => section_vals_get_subs_vals(force_env%root_section, "MOTION")
200 2 : md_section => section_vals_get_subs_vals(motion_section, "MD")
201 :
202 : ! Real call to MD driver - Low Level
203 2 : ALLOCATE (md_env)
204 2 : CALL md_env_create(md_env, md_section, para_env, force_env=force_env)
205 2 : IF (PRESENT(averages)) CALL set_md_env(md_env, averages=averages)
206 :
207 2 : CPASSERT(ASSOCIATED(globenv))
208 2 : CPASSERT(ASSOCIATED(force_env))
209 :
210 2 : NULLIFY (particles, cell, simpar, itimes, used_time, subsys, &
211 2 : md_ener, thermostats, barostat, reftraj, force_env_section, &
212 2 : reftraj_section, work_section, atomic_kinds, &
213 2 : local_particles, time, fe_env, free_energy_section, &
214 2 : constraint_section, thermal_regions, subsys_i)
215 2 : logger => cp_get_default_logger()
216 2 : para_env => force_env%para_env
217 :
218 2 : global_section => section_vals_get_subs_vals(force_env%root_section, "GLOBAL")
219 2 : free_energy_section => section_vals_get_subs_vals(motion_section, "FREE_ENERGY")
220 2 : constraint_section => section_vals_get_subs_vals(motion_section, "CONSTRAINT")
221 2 : CALL section_vals_val_get(global_section, "SAVE_MEM", l_val=save_mem)
222 :
223 2 : CALL section_vals_val_get(global_section, "RUN_TYPE", i_val=run_type_id)
224 :
225 2 : CALL create_simpar_type(simpar)
226 2 : force_env_section => force_env%force_env_section
227 2 : subsys_section => section_vals_get_subs_vals(force_env_section, "SUBSYS")
228 2 : CALL cp_add_iter_level(logger%iter_info, "MD")
229 2 : CALL cp_iterate(logger%iter_info, iter_nr=initialStep)
230 : ! Read MD section
231 2 : CALL read_md_section(simpar, motion_section, md_section)
232 : ! Setup print_keys
233 : simpar%info_constraint = cp_print_key_unit_nr(logger, constraint_section, &
234 2 : "CONSTRAINT_INFO", extension=".shakeLog", log_filename=.FALSE.)
235 : simpar%lagrange_multipliers = cp_print_key_unit_nr(logger, constraint_section, &
236 2 : "LAGRANGE_MULTIPLIERS", extension=".LagrangeMultLog", log_filename=.FALSE.)
237 : simpar%dump_lm = BTEST(cp_print_key_should_output(logger%iter_info, constraint_section, &
238 2 : "LAGRANGE_MULTIPLIERS"), cp_p_file)
239 :
240 : ! Create the structure for the md energies
241 8 : ALLOCATE (md_ener)
242 2 : CALL create_md_ener(md_ener)
243 2 : CALL set_md_env(md_env, md_ener=md_ener)
244 :
245 : ! If requested setup Thermostats
246 : CALL create_thermostats(thermostats, md_section, force_env, simpar, para_env, &
247 2 : globenv, global_section)
248 :
249 : ! If requested setup Barostat
250 2 : CALL create_barostat_type(barostat, md_section, force_env, simpar, globenv)
251 :
252 : ! If requested setup different thermal regions
253 2 : CALL create_thermal_regions(thermal_regions, md_section, simpar, force_env)
254 :
255 2 : CALL set_md_env(md_env, thermostats=thermostats, barostat=barostat, thermal_regions=thermal_regions)
256 :
257 2 : IF (simpar%ensemble == reftraj_ensemble) THEN
258 0 : reftraj_section => section_vals_get_subs_vals(md_section, "REFTRAJ")
259 0 : ALLOCATE (reftraj)
260 0 : CALL create_reftraj(reftraj, reftraj_section, para_env)
261 0 : CALL set_md_env(md_env, reftraj=reftraj)
262 : END IF
263 :
264 : CALL force_env_get(force_env, subsys=subsys, cell=cell, &
265 2 : force_env_section=force_env_section)
266 :
267 : ! Set V0 if needed
268 2 : IF (simpar%ensemble == nph_uniaxial_ensemble .OR. simpar%ensemble == nph_uniaxial_damped_ensemble) THEN
269 0 : IF (simpar%v0 == 0._dp) simpar%v0 = cell%deth
270 : END IF
271 :
272 : ! Setup Free Energy Calculation (if required)
273 2 : CALL fe_env_create(fe_env, free_energy_section)
274 : CALL set_md_env(md_env=md_env, simpar=simpar, fe_env=fe_env, cell=cell, &
275 2 : force_env=force_env)
276 :
277 : ! Possibly initialize Wiener processes
278 2 : IF (simpar%ensemble == langevin_ensemble) CALL create_wiener_process(md_env)
279 2 : time_iter_start = m_walltime()
280 :
281 : CALL get_md_env(md_env, force_env=force_env, itimes=itimes, constant=constant, &
282 2 : md_ener=md_ener, t=time, used_time=used_time)
283 :
284 : ! Attach the time counter of the meta_env to the one of the MD
285 2 : CALL set_meta_env(force_env%meta_env, time=time)
286 : ! Initialize the md_ener structure
287 :
288 2 : force_env%meta_env%dt = force_env%meta_env%zdt
289 2 : CALL initialize_md_ener(md_ener, force_env, simpar)
290 :
291 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! MC setup up
292 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
293 :
294 2 : NULLIFY (mc_env, mc_par, MCaverages)
295 :
296 2 : CALL section_vals_get(force_env_section, n_repetition=isos)
297 2 : CPASSERT(isos == 1)
298 : ! set some values...will use get_globenv if that ever comes around
299 :
300 : ! initialize the random numbers
301 : rng_stream_mc = rng_stream_type(name="Random numbers for monte carlo acc/rej", &
302 2 : distribution_type=UNIFORM)
303 :
304 2 : NULLIFY (mc_section)
305 2 : ALLOCATE (mc_par)
306 :
307 : mc_section => section_vals_get_subs_vals(force_env%root_section, &
308 2 : "MOTION%MC")
309 : CALL section_vals_val_get(mc_section, "ENSEMBLE", &
310 2 : c_val=ensemble)
311 2 : CPASSERT(str_comp(ensemble, "TRADITIONAL"))
312 : CALL section_vals_val_get(mc_section, "NSTEP", &
313 2 : i_val=nmccycles)
314 2 : CPASSERT(nmccycles > 0)
315 : CALL section_vals_val_get(mc_section, "IPRINT", &
316 2 : i_val=iprint)
317 2 : CALL section_vals_val_get(mc_section, "RANDOMTOSKIP", i_val=rand2skip)
318 2 : CPASSERT(rand2skip >= 0)
319 2 : temp = cp_unit_from_cp2k(simpar%temp_ext, "K")
320 :
321 : CALL set_mc_par(mc_par, ensemble=ensemble, nstep=nmccycles, iprint=iprint, temperature=temp, &
322 : beta=1.0_dp/temp/boltzmann*joule, exp_max_val=0.9_dp*LOG(HUGE(0.0_dp)), &
323 : exp_min_val=0.9_dp*LOG(TINY(0.0_dp)), max_val=HUGE(0.0_dp), min_val=0.0_dp, &
324 2 : source=para_env%source, group=para_env, ionode=para_env%is_source(), rand2skip=rand2skip)
325 :
326 2 : output_unit = cp_logger_get_default_io_unit(logger)
327 2 : IF (output_unit > 0) THEN
328 1 : WRITE (output_unit, '(a,a)') "HMC| Hybrid Monte Carlo Scheme "
329 1 : WRITE (output_unit, '(a,a)') "HMC| Ensemble ", ADJUSTL(ensemble)
330 1 : WRITE (output_unit, '(a,i0)') "HMC| MC Cycles ", nmccycles
331 1 : WRITE (output_unit, '(a,i0,a)') "HMC| Print every ", iprint, " cycles"
332 1 : WRITE (output_unit, '(a,i0)') "HMC| Number of random numbers to skip ", rand2skip
333 1 : WRITE (output_unit, '(a,f16.8,a)') "HMC| Temperature ", temp, "K"
334 : END IF
335 :
336 2 : CALL force_env_get(force_env, subsys=subsys)
337 :
338 : CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, &
339 2 : particles=particles, virial=virial)
340 :
341 2 : DO i = 1, rand2skip
342 0 : auxRandom = rng_stream_mc%next()
343 2 : DO j = 1, 3*SIZE(particles%els)
344 0 : auxRandom = globenv%gaussian_rng_stream%next()
345 : END DO
346 : END DO
347 :
348 2 : ALLOCATE (mc_env)
349 2 : CALL mc_env_create(mc_env)
350 2 : CALL set_mc_env(mc_env, mc_par=mc_par, force_env=force_env)
351 : !!!!!!!end mc setup
352 :
353 : ! Check for ensembles requiring the stress tensor - takes into account the possibility for
354 : ! multiple force_evals
355 : IF ((simpar%ensemble == npt_i_ensemble) .OR. &
356 : (simpar%ensemble == npt_ia_ensemble) .OR. &
357 : (simpar%ensemble == npt_f_ensemble) .OR. &
358 : (simpar%ensemble == npe_f_ensemble) .OR. &
359 : (simpar%ensemble == npe_i_ensemble) .OR. &
360 2 : (simpar%ensemble == nph_uniaxial_ensemble) .OR. &
361 : (simpar%ensemble == nph_uniaxial_damped_ensemble)) THEN
362 0 : check = virial%pv_availability
363 0 : IF (.NOT. check) THEN
364 : CALL cp_abort(__LOCATION__, &
365 : "Virial evaluation not requested for this run in the input file! "// &
366 : "You may consider to switch on the virial evaluation with the keyword: STRESS_TENSOR. "// &
367 0 : "Be sure the method you are using can compute the virial!")
368 : END IF
369 0 : IF (ASSOCIATED(force_env%sub_force_env)) THEN
370 0 : DO i = 1, SIZE(force_env%sub_force_env)
371 0 : IF (ASSOCIATED(force_env%sub_force_env(i)%force_env)) THEN
372 0 : CALL force_env_get(force_env%sub_force_env(i)%force_env, subsys=subsys_i)
373 0 : CALL cp_subsys_get(subsys_i, virial=virial)
374 0 : check = check .AND. virial%pv_availability
375 : END IF
376 : END DO
377 : END IF
378 0 : IF (.NOT. check) THEN
379 : CALL cp_abort(__LOCATION__, &
380 : "Virial evaluation not requested for all the force_eval sections present in"// &
381 : " the input file! You have to switch on the virial evaluation with the keyword: STRESS_TENSOR"// &
382 0 : " in each force_eval section. Be sure the method you are using can compute the virial!")
383 : END IF
384 : END IF
385 :
386 : ! Computing Forces at zero MD step
387 2 : IF (simpar%ensemble /= reftraj_ensemble) THEN
388 2 : CALL section_vals_val_get(md_section, "STEP_START_VAL", i_val=itimes)
389 2 : CALL section_vals_val_get(md_section, "TIME_START_VAL", r_val=time)
390 2 : CALL section_vals_val_get(md_section, "ECONS_START_VAL", r_val=constant)
391 2 : CALL section_vals_val_set(md_section, "STEP_START_VAL", i_val=initialStep)
392 2 : CALL section_vals_val_set(md_section, "TIME_START_VAL", r_val=inittime)
393 2 : initialStep = itimes
394 2 : CALL cp_iterate(logger%iter_info, iter_nr=itimes)
395 2 : IF (save_mem) THEN
396 0 : work_section => section_vals_get_subs_vals(subsys_section, "VELOCITY")
397 0 : CALL section_vals_remove_values(work_section)
398 0 : work_section => section_vals_get_subs_vals(subsys_section, "SHELL_VELOCITY")
399 0 : CALL section_vals_remove_values(work_section)
400 0 : work_section => section_vals_get_subs_vals(subsys_section, "CORE_VELOCITY")
401 0 : CALL section_vals_remove_values(work_section)
402 : END IF
403 :
404 2 : meta_env_saved => force_env%meta_env
405 2 : NULLIFY (force_env%meta_env)
406 2 : CALL force_env_calc_energy_force(force_env, calc_force=.FALSE.)
407 2 : force_env%meta_env => meta_env_saved
408 :
409 2 : IF (ASSOCIATED(force_env%qs_env)) THEN
410 2 : force_env%qs_env%sim_time = 0.0_dp
411 2 : force_env%qs_env%sim_step = 0
412 : END IF
413 : ! Warm-up engines for metadynamics
414 2 : IF (ASSOCIATED(force_env%meta_env)) THEN
415 2 : IF (force_env%meta_env%langevin) THEN
416 2 : CALL create_wiener_process_cv(force_env%meta_env)
417 2 : DO j = 1, (rand2skip - 1)/nmccycles
418 2 : DO i = 1, force_env%meta_env%n_colvar
419 0 : auxRandom = force_env%meta_env%rng(i)%next()
420 0 : auxRandom = force_env%meta_env%rng(i)%next()
421 : END DO
422 : END DO
423 : END IF
424 2 : CALL tamc_force(force_env)
425 : END IF
426 :
427 2 : IF (simpar%do_respa) THEN
428 : CALL force_env_calc_energy_force(force_env%sub_force_env(1)%force_env, &
429 0 : calc_force=.TRUE.)
430 : END IF
431 :
432 : CALL virial_evaluate(atomic_kinds%els, particles%els, local_particles, &
433 2 : virial, force_env%para_env)
434 :
435 2 : CALL md_energy(md_env, md_ener)
436 2 : md_stride = 1
437 : ELSE
438 0 : CALL get_md_env(md_env, reftraj=reftraj)
439 0 : CALL initialize_reftraj(reftraj, reftraj_section, md_env)
440 0 : itimes = reftraj%info%first_snapshot - 1
441 0 : md_stride = reftraj%info%stride
442 : END IF
443 :
444 : CALL cp_print_key_finished_output(simpar%info_constraint, logger, &
445 2 : constraint_section, "CONSTRAINT_INFO")
446 : CALL cp_print_key_finished_output(simpar%lagrange_multipliers, logger, &
447 2 : constraint_section, "LAGRANGE_MULTIPLIERS")
448 2 : CALL init_mc_moves(moves)
449 2 : CALL init_mc_moves(gmoves)
450 6 : ALLOCATE (r(1:3, SIZE(particles%els)))
451 2 : CALL mc_averages_create(MCaverages)
452 : !!!!! some more buffers
453 : ! Allocate random number for Langevin Thermostat acting on COLVARS
454 6 : ALLOCATE (xieta(2*force_env%meta_env%n_colvar))
455 2 : xieta(:) = 0.0_dp
456 6 : ALLOCATE (An(force_env%meta_env%n_colvar))
457 2 : An(:) = 0.0_dp
458 4 : ALLOCATE (fz(force_env%meta_env%n_colvar))
459 2 : fz(:) = 0.0_dp
460 4 : ALLOCATE (zbuff(2*force_env%meta_env%n_colvar))
461 2 : zbuff(:) = 0.0_dp
462 :
463 2 : IF (output_unit > 0) THEN
464 1 : WRITE (output_unit, '(a)') "HMC|==== Initial average forces"
465 : END IF
466 2 : CALL metadyn_write_colvar_header(force_env)
467 2 : moves%hmc%attempts = 0
468 2 : moves%hmc%successes = 0
469 2 : gmoves%hmc%attempts = 0
470 2 : gmoves%hmc%successes = 0
471 2 : IF (initialStep == 0) THEN
472 : !!! if we come from a restart we shall properly compute the average force
473 : !!! read the average force up to now
474 4 : DO i = 1, force_env%meta_env%n_colvar
475 2 : fs_section => section_vals_get_subs_vals(force_env%meta_env%metadyn_section, "EXT_LAGRANGE_FS")
476 2 : CALL section_vals_get(fs_section, explicit=explicit)
477 4 : IF (explicit) THEN
478 : CALL section_vals_val_get(fs_section, "_DEFAULT_KEYWORD_", &
479 0 : i_rep_val=i, r_val=rval)
480 0 : fz(i) = rval*rand2skip
481 : END IF
482 : END DO
483 :
484 : CALL HMCsampler(globenv, force_env, MCaverages, r, mc_par, moves, gmoves, rng_stream_mc, output_unit, &
485 2 : fz, zbuff, nskip=rand2skip)
486 2 : CALL cp_iterate(logger%iter_info, last=.FALSE., iter_nr=0)
487 2 : CALL section_vals_val_set(mc_section, "RANDOMTOSKIP", i_val=rand2skip + nmccycles)
488 2 : CALL write_restart(md_env=md_env, root_section=force_env%root_section)
489 : END IF
490 2 : IF (output_unit > 0) THEN
491 1 : WRITE (output_unit, '(a)') "HMC|==== end initial average forces"
492 : END IF
493 :
494 2 : CALL metadyn_write_colvar(force_env)
495 :
496 4 : DO istep = 1, force_env%meta_env%TAMCSteps
497 : ! Increase counters
498 2 : itimes = itimes + 1
499 2 : time = time + force_env%meta_env%dt
500 2 : IF (output_unit > 0) THEN
501 1 : WRITE (output_unit, '(a)') "HMC|==================================="
502 1 : WRITE (output_unit, '(a,1x,i0)') "HMC| on z step ", istep
503 : END IF
504 : !needed when electric field fields are applied
505 2 : IF (ASSOCIATED(force_env%qs_env)) THEN
506 2 : force_env%qs_env%sim_time = time
507 2 : force_env%qs_env%sim_step = itimes
508 2 : force_env%meta_env%time = force_env%qs_env%sim_time
509 : END IF
510 :
511 2 : CALL cp_iterate(logger%iter_info, last=(istep == force_env%meta_env%TAMCSteps), iter_nr=itimes)
512 : ! Open possible Shake output units
513 : simpar%info_constraint = cp_print_key_unit_nr(logger, constraint_section, "CONSTRAINT_INFO", &
514 2 : extension=".shakeLog", log_filename=.FALSE.)
515 : simpar%lagrange_multipliers = cp_print_key_unit_nr( &
516 : logger, constraint_section, &
517 2 : "LAGRANGE_MULTIPLIERS", extension=".LagrangeMultLog", log_filename=.FALSE.)
518 : simpar%dump_lm = BTEST(cp_print_key_should_output(logger%iter_info, constraint_section, &
519 2 : "LAGRANGE_MULTIPLIERS"), cp_p_file)
520 :
521 : ! Velocity Verlet Integrator
522 :
523 2 : moves%hmc%attempts = 0
524 2 : moves%hmc%successes = 0
525 : CALL langevinVEC(md_env, globenv, mc_env, moves, gmoves, r, &
526 2 : rng_stream_mc, xieta, An, fz, MCaverages, zbuff)
527 :
528 : ! Close Shake output if requested...
529 : CALL cp_print_key_finished_output(simpar%info_constraint, logger, &
530 2 : constraint_section, "CONSTRAINT_INFO")
531 : CALL cp_print_key_finished_output(simpar%lagrange_multipliers, logger, &
532 2 : constraint_section, "LAGRANGE_MULTIPLIERS")
533 2 : CALL cp_iterate(logger%iter_info, iter_nr=initialStep)
534 2 : CALL metadyn_write_colvar(force_env)
535 : ! Free Energy calculation
536 : ! CALL free_energy_evaluate(md_env,should_stop,free_energy_section)
537 :
538 : ! Test for <PROJECT_NAME>.EXIT_MD or for WALL_TIME to exit
539 : ! Default:
540 : ! IF so we don't overwrite the restart or append to the trajectory
541 : ! because the execution could in principle stop inside the SCF where energy
542 : ! and forces are not converged.
543 : ! But:
544 : ! You can force to print the last step (for example if the method used
545 : ! to compute energy and forces is not SCF based) activating the print_key
546 : ! MOTION%MD%PRINT%FORCE_LAST.
547 2 : CALL external_control(should_stop, "MD", globenv=globenv)
548 2 : IF (should_stop) THEN
549 0 : CALL cp_iterate(logger%iter_info, last=.TRUE., iter_nr=itimes)
550 0 : EXIT
551 : END IF
552 :
553 2 : time_iter_stop = m_walltime()
554 2 : used_time = time_iter_stop - time_iter_start
555 2 : time_iter_start = time_iter_stop
556 :
557 6 : IF (output_unit > 0) THEN
558 1 : WRITE (output_unit, '(a,1x,i0)') "HMC| end z step ", istep
559 1 : WRITE (output_unit, '(a)') "HMC|==================================="
560 : END IF
561 : END DO
562 2 : CALL cp_iterate(logger%iter_info, last=.TRUE., iter_nr=itimes)
563 2 : force_env%qs_env%sim_time = 0.0_dp
564 2 : force_env%qs_env%sim_step = 0
565 2 : rand2skip = rand2skip + nmccycles*force_env%meta_env%TAMCSteps
566 2 : IF (initialStep == 0) rand2skip = rand2skip + nmccycles
567 2 : CALL section_vals_val_set(mc_section, "RANDOMTOSKIP", i_val=rand2skip)
568 :
569 2 : CALL write_restart(md_env=md_env, root_section=force_env%root_section)
570 : ! if we need the final kinetic energy for Hybrid Monte Carlo
571 :
572 : ! Remove the iteration level
573 2 : CALL cp_rm_iter_level(logger%iter_info, "MD")
574 :
575 : ! Deallocate Thermostats and Barostats
576 2 : CALL release_simpar_type(simpar)
577 :
578 2 : CALL md_env_release(md_env)
579 2 : DEALLOCATE (md_env)
580 : ! Clean restartable sections..
581 2 : IF (my_rm_restart_info) CALL remove_restart_info(force_env%root_section)
582 2 : CALL MC_ENV_RELEASE(mc_env)
583 2 : DEALLOCATE (mc_env)
584 2 : DEALLOCATE (mc_par)
585 2 : CALL MC_MOVES_RELEASE(moves)
586 2 : CALL MC_MOVES_RELEASE(gmoves)
587 2 : DEALLOCATE (r)
588 2 : DEALLOCATE (xieta)
589 2 : DEALLOCATE (An)
590 2 : DEALLOCATE (fz)
591 2 : DEALLOCATE (zbuff)
592 2 : CALL mc_averages_release(MCaverages)
593 2 : CALL timestop(handle)
594 :
595 64 : END SUBROUTINE qs_tamc
596 :
597 : ! **************************************************************************************************
598 : !> \brief Propagates velocities for z half a step
599 : !> \param force_env ...
600 : !> \param An ...
601 : !> \author Alin M Elena
602 : !> \note Vanden-Eijnden Ciccotti C.Phys.Letter 429 (2006) 310-316
603 : ! **************************************************************************************************
604 4 : SUBROUTINE tamc_velocities_colvar(force_env, An)
605 : TYPE(force_env_type), POINTER :: force_env
606 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: An
607 :
608 : CHARACTER(len=*), PARAMETER :: routineN = 'tamc_velocities_colvar'
609 :
610 : INTEGER :: handle, i_c
611 : REAL(kind=dp) :: dt, fft, sigma
612 : TYPE(cp_logger_type), POINTER :: logger
613 : TYPE(meta_env_type), POINTER :: meta_env
614 : TYPE(metavar_type), POINTER :: cv
615 :
616 4 : NULLIFY (logger, meta_env, cv)
617 4 : meta_env => force_env%meta_env
618 4 : CALL timeset(routineN, handle)
619 4 : logger => cp_get_default_logger()
620 : ! Add citation
621 4 : IF (meta_env%langevin) CALL cite_reference(VandenCic2006)
622 4 : dt = meta_env%dt
623 :
624 : ! Evolve Velocities
625 4 : meta_env%epot_walls = 0.0_dp
626 8 : DO i_c = 1, meta_env%n_colvar
627 4 : cv => meta_env%metavar(i_c)
628 4 : fft = cv%ff_s + cv%ff_hills
629 4 : sigma = SQRT((meta_env%temp_wanted*kelvin)*2.0_dp*(boltzmann/joule)*cv%gamma/cv%mass)
630 4 : cv%vvp = cv%vvp + 0.5_dp*dt*(fft/cv%mass - cv%gamma*cv%vvp)*(1.0_dp - 0.25_dp*dt*cv%gamma) + An(i_c)
631 8 : meta_env%epot_walls = meta_env%epot_walls + cv%epot_walls
632 : END DO
633 4 : CALL timestop(handle)
634 4 : END SUBROUTINE tamc_velocities_colvar
635 :
636 : ! **************************************************************************************************
637 : !> \brief propagates z one step
638 : !> \param force_env ...
639 : !> \param xieta ...
640 : !> \author Alin M Elena
641 : !> \note Vanden-Eijnden Ciccotti C.Phys.Letter 429 (2006) 310-316
642 : ! **************************************************************************************************
643 2 : SUBROUTINE tamc_position_colvar(force_env, xieta)
644 : TYPE(force_env_type), POINTER :: force_env
645 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: xieta
646 :
647 : CHARACTER(len=*), PARAMETER :: routineN = 'tamc_position_colvar'
648 :
649 : INTEGER :: handle, i_c
650 : REAL(kind=dp) :: dt, sigma
651 : TYPE(cp_logger_type), POINTER :: logger
652 : TYPE(meta_env_type), POINTER :: meta_env
653 : TYPE(metavar_type), POINTER :: cv
654 :
655 2 : NULLIFY (logger, meta_env, cv)
656 2 : meta_env => force_env%meta_env
657 : ! IF (.NOT.ASSOCIATED(meta_env)) RETURN
658 :
659 2 : CALL timeset(routineN, handle)
660 2 : logger => cp_get_default_logger()
661 :
662 : ! Add citation
663 2 : IF (meta_env%langevin) CALL cite_reference(VandenCic2006)
664 2 : dt = meta_env%dt
665 :
666 : ! Update of ss0
667 4 : DO i_c = 1, meta_env%n_colvar
668 2 : cv => meta_env%metavar(i_c)
669 2 : sigma = SQRT((meta_env%temp_wanted*kelvin)*2.0_dp*(boltzmann/joule)*cv%gamma/cv%mass)
670 : ! cv%ss0 =cv%ss0 +dt*cv%vvp
671 2 : cv%ss0 = cv%ss0 + dt*cv%vvp + dt*SQRT(dt/12.0_dp)*sigma*xieta(i_c + meta_env%n_colvar)
672 4 : IF (cv%periodic) THEN
673 : ! A periodic COLVAR is always within [-pi,pi]
674 0 : cv%ss0 = SIGN(1.0_dp, ASIN(SIN(cv%ss0)))*ACOS(COS(cv%ss0))
675 : END IF
676 : END DO
677 2 : CALL timestop(handle)
678 :
679 2 : END SUBROUTINE tamc_position_colvar
680 :
681 : ! **************************************************************************************************
682 : !> \brief Computes forces on z
683 : !> #details also can be used to get the potenzial evergy of z
684 : !> \param force_env ...
685 : !> \param zpot ...
686 : !> \author Alin M Elena
687 : ! **************************************************************************************************
688 10 : SUBROUTINE tamc_force(force_env, zpot)
689 : TYPE(force_env_type), POINTER :: force_env
690 : REAL(KIND=dp), INTENT(inout), OPTIONAL :: zpot
691 :
692 : CHARACTER(len=*), PARAMETER :: routineN = 'tamc_force'
693 :
694 : INTEGER :: handle, i, i_c, icolvar, ii
695 : LOGICAL :: explicit
696 : REAL(kind=dp) :: diff_ss, dt, rval
697 10 : TYPE(colvar_p_type), DIMENSION(:), POINTER :: colvar_p
698 : TYPE(cp_logger_type), POINTER :: logger
699 : TYPE(cp_subsys_type), POINTER :: subsys
700 : TYPE(meta_env_type), POINTER :: meta_env
701 : TYPE(metavar_type), POINTER :: cv
702 : TYPE(particle_list_type), POINTER :: particles
703 : TYPE(section_vals_type), POINTER :: ss0_section, ss_section, vvp_section
704 :
705 10 : NULLIFY (logger, meta_env)
706 10 : meta_env => force_env%meta_env
707 : ! IF (.NOT.ASSOCIATED(meta_env)) RETURN
708 :
709 10 : CALL timeset(routineN, handle)
710 10 : logger => cp_get_default_logger()
711 10 : NULLIFY (colvar_p, subsys, cv, ss0_section, vvp_section, ss_section)
712 10 : CALL force_env_get(force_env, subsys=subsys)
713 :
714 10 : dt = meta_env%dt
715 10 : IF (.NOT. meta_env%restart) meta_env%n_steps = meta_env%n_steps + 1
716 : ! compute ss and the derivative of ss with respect to the atomic positions
717 20 : DO i_c = 1, meta_env%n_colvar
718 10 : cv => meta_env%metavar(i_c)
719 10 : icolvar = cv%icolvar
720 10 : CALL colvar_eval_glob_f(icolvar, force_env)
721 10 : cv%ss = subsys%colvar_p(icolvar)%colvar%ss
722 : ! Restart for Extended Lagrangian Metadynamics
723 20 : IF (meta_env%restart) THEN
724 : ! Initialize the position of the collective variable in the extended lagrange
725 2 : ss0_section => section_vals_get_subs_vals(meta_env%metadyn_section, "EXT_LAGRANGE_SS0")
726 2 : CALL section_vals_get(ss0_section, explicit=explicit)
727 2 : IF (explicit) THEN
728 : CALL section_vals_val_get(ss0_section, "_DEFAULT_KEYWORD_", &
729 2 : i_rep_val=i_c, r_val=rval)
730 2 : cv%ss0 = rval
731 : ELSE
732 0 : cv%ss0 = cv%ss
733 : END IF
734 2 : vvp_section => section_vals_get_subs_vals(meta_env%metadyn_section, "EXT_LAGRANGE_VVP")
735 2 : CALL section_vals_get(vvp_section, explicit=explicit)
736 2 : IF (explicit) THEN
737 0 : CALL setup_velocities_z(force_env)
738 : CALL section_vals_val_get(vvp_section, "_DEFAULT_KEYWORD_", &
739 0 : i_rep_val=i_c, r_val=rval)
740 0 : cv%vvp = rval
741 : ELSE
742 2 : CALL setup_velocities_z(force_env)
743 : END IF
744 2 : ss_section => section_vals_get_subs_vals(meta_env%metadyn_section, "EXT_LAGRANGE_SS")
745 2 : CALL section_vals_get(ss_section, explicit=explicit)
746 2 : IF (explicit) THEN
747 : CALL section_vals_val_get(ss_section, "_DEFAULT_KEYWORD_", &
748 0 : i_rep_val=i_c, r_val=rval)
749 0 : cv%ss = rval
750 : END IF
751 : END IF
752 : !
753 : END DO
754 : ! forces on the atoms
755 10 : NULLIFY (particles)
756 : CALL cp_subsys_get(subsys, colvar_p=colvar_p, &
757 10 : particles=particles)
758 :
759 10 : meta_env%restart = .FALSE.
760 10 : meta_env%epot_s = 0.0_dp
761 10 : meta_env%epot_walls = 0.0_dp
762 20 : DO i_c = 1, meta_env%n_colvar
763 10 : cv => meta_env%metavar(i_c)
764 10 : diff_ss = cv%ss - cv%ss0
765 10 : IF (cv%periodic) THEN
766 : ! The difference of a periodic COLVAR is always within [-pi,pi]
767 0 : diff_ss = SIGN(1.0_dp, ASIN(SIN(diff_ss)))*ACOS(COS(diff_ss))
768 : END IF
769 10 : cv%epot_s = 0.5_dp*cv%lambda*diff_ss*diff_ss
770 10 : cv%ff_s = cv%lambda*(diff_ss)
771 10 : meta_env%epot_s = meta_env%epot_s + cv%epot_s
772 10 : icolvar = cv%icolvar
773 :
774 50 : DO ii = 1, colvar_p(icolvar)%colvar%n_atom_s
775 30 : i = colvar_p(icolvar)%colvar%i_atom(ii)
776 220 : particles%els(i)%f = particles%els(i)%f - cv%ff_s*colvar_p(icolvar)%colvar%dsdr(:, ii)
777 : END DO
778 :
779 : END DO
780 10 : IF (PRESENT(zpot)) zpot = meta_env%epot_s
781 10 : CALL fix_atom_control(force_env)
782 :
783 10 : CALL timestop(handle)
784 10 : END SUBROUTINE tamc_force
785 :
786 : ! **************************************************************************************************
787 : !> \brief propagates one time step both z systems and samples the x system
788 : !> \param md_env ...
789 : !> \param globenv ...
790 : !> \param mc_env ...
791 : !> \param moves ...
792 : !> \param gmoves ...
793 : !> \param r ...
794 : !> \param rng_stream_mc ...
795 : !> \param xieta ...
796 : !> \param An ...
797 : !> \param fz ...
798 : !> \param averages ...
799 : !> \param zbuff ...
800 : !> \author Alin M Elena
801 : ! **************************************************************************************************
802 8 : SUBROUTINE langevinVEC(md_env, globenv, mc_env, moves, gmoves, r, &
803 2 : rng_stream_mc, xieta, An, fz, averages, zbuff)
804 :
805 : TYPE(md_environment_type), POINTER :: md_env
806 : TYPE(global_environment_type), POINTER :: globenv
807 : TYPE(mc_environment_type), POINTER :: mc_env
808 : TYPE(mc_moves_type), POINTER :: moves, gmoves
809 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: r
810 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream_mc
811 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: xieta, An, fz
812 : TYPE(mc_averages_type), INTENT(INOUT), POINTER :: averages
813 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: zbuff
814 :
815 : INTEGER :: iprint, ivar, nparticle, nparticle_kind, &
816 : nstep, output_unit
817 : REAL(KIND=dp) :: dt, gamma, mass, sigma
818 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
819 2 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
820 : TYPE(cell_type), POINTER :: cell
821 : TYPE(cp_logger_type), POINTER :: logger
822 : TYPE(cp_subsys_type), POINTER :: subsys
823 : TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles
824 : TYPE(force_env_type), POINTER :: force_env
825 : TYPE(global_constraint_type), POINTER :: gci
826 : TYPE(mc_simpar_type), POINTER :: mc_par
827 : TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
828 2 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
829 : TYPE(molecule_list_type), POINTER :: molecules
830 2 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
831 : TYPE(mp_para_env_type), POINTER :: para_env
832 : TYPE(particle_list_type), POINTER :: particles
833 2 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
834 : TYPE(simpar_type), POINTER :: simpar
835 : TYPE(virial_type), POINTER :: virial
836 :
837 2 : NULLIFY (logger, mc_par)
838 4 : logger => cp_get_default_logger()
839 2 : output_unit = cp_logger_get_default_io_unit(logger)
840 :
841 : ! quantitites to be nullified for the get_md_env
842 2 : NULLIFY (simpar, force_env, para_env)
843 : ! quantities to be nullified for the force_env_get environment
844 2 : NULLIFY (subsys, cell)
845 : ! quantitites to be nullified for the cp_subsys_get
846 2 : NULLIFY (atomic_kinds, local_particles, particles, local_molecules, molecules, molecule_kinds, gci)
847 :
848 : CALL get_md_env(md_env=md_env, simpar=simpar, force_env=force_env, &
849 2 : para_env=para_env)
850 2 : CALL get_mc_env(mc_env, mc_par=mc_par)
851 2 : CALL get_mc_par(mc_par, nstep=nstep, iprint=iprint)
852 :
853 2 : dt = simpar%dt
854 2 : CALL force_env_get(force_env=force_env, subsys=subsys, cell=cell)
855 :
856 : !!!! this bit should vanish once I understand what the hell is with it
857 :
858 : ! ! Do some checks on coordinates and box
859 2 : CALL apply_qmmm_walls_reflective(force_env)
860 :
861 : CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, &
862 : particles=particles, local_molecules=local_molecules, molecules=molecules, &
863 2 : molecule_kinds=molecule_kinds, gci=gci, virial=virial)
864 :
865 2 : nparticle_kind = atomic_kinds%n_els
866 2 : atomic_kind_set => atomic_kinds%els
867 2 : molecule_kind_set => molecule_kinds%els
868 :
869 2 : nparticle = particles%n_els
870 2 : particle_set => particles%els
871 2 : molecule_set => molecules%els
872 2 : CPASSERT(ASSOCIATED(force_env%meta_env))
873 2 : CPASSERT(force_env%meta_env%langevin)
874 : ! *** Velocity Verlet for Langevin *** v(t)--> v(t+1/2)
875 : !!!!!! noise xi is in the first half, eta in the second half
876 4 : DO ivar = 1, force_env%meta_env%n_colvar
877 2 : xieta(ivar) = force_env%meta_env%rng(ivar)%next()
878 2 : xieta(ivar + force_env%meta_env%n_colvar) = force_env%meta_env%rng(ivar)%next()
879 2 : gamma = force_env%meta_env%metavar(ivar)%gamma
880 2 : mass = force_env%meta_env%metavar(ivar)%mass
881 2 : sigma = SQRT((force_env%meta_env%temp_wanted*kelvin)*2.0_dp*(boltzmann/joule)*gamma/mass)
882 : An(ivar) = 0.5_dp*SQRT(dt)*sigma*(xieta(ivar)*(1.0_dp - 0.5_dp*dt*gamma) - &
883 4 : dt*gamma*xieta(ivar + force_env%meta_env%n_colvar)/SQRT(12.0_dp))
884 : END DO
885 : ! *** Velocity Verlet for Langeving *** v(t)--> v(t+1/2)
886 2 : CALL tamc_velocities_colvar(force_env, An)
887 : ! *** Velocity Verlet for Langevin S(t)->S(t+1)
888 2 : CALL tamc_position_colvar(force_env, xieta)
889 : !!!!! start zHMC sampler
890 2 : CALL HMCsampler(globenv, force_env, averages, r, mc_par, moves, gmoves, rng_stream_mc, output_unit, fz, zbuff)
891 :
892 : ! CALL final_mc_write(mc_par,tmp_moves,&
893 : ! output_unit,energy_check,&
894 : ! initial_energy,final_energy,&
895 : ! averages)
896 :
897 : !!!!!!!!!!!!!!!!!!!! end zHMC sampler
898 : ! *** Velocity Verlet for Langeving *** v(t+1/2)--> v(t+1)
899 2 : CALL tamc_velocities_colvar(force_env, An)
900 : ! CALL virial_evaluate ( atomic_kind_set, particle_set, &
901 : ! local_particles, virial, para_env)
902 :
903 2 : END SUBROUTINE langevinVEC
904 :
905 : ! **************************************************************************************************
906 : !> \brief Driver routin for the canonical sampler using modified HMC
907 : !> \param globenv ...
908 : !> \param force_env ...
909 : !> \param averages ...
910 : !> \param r ...
911 : !> \param mc_par ...
912 : !> \param moves ...
913 : !> \param gmoves ...
914 : !> \param rng_stream_mc ...
915 : !> \param output_unit ...
916 : !> \param fz ...
917 : !> \param zbuff ...
918 : !> \param nskip ...
919 : !> \author Alin M Elena
920 : !> \note at the end of this routine %ff_s will contain mean force
921 : ! **************************************************************************************************
922 :
923 20 : SUBROUTINE HMCsampler(globenv, force_env, averages, r, mc_par, moves, gmoves, rng_stream_mc, output_unit, &
924 4 : fz, zbuff, nskip)
925 : TYPE(global_environment_type), POINTER :: globenv
926 : TYPE(force_env_type), POINTER :: force_env
927 : TYPE(mc_averages_type), POINTER :: averages
928 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: r
929 : TYPE(mc_simpar_type), POINTER :: mc_par
930 : TYPE(mc_moves_type), POINTER :: moves, gmoves
931 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream_mc
932 : INTEGER, INTENT(IN) :: output_unit
933 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: fz, zbuff
934 : INTEGER, INTENT(IN), OPTIONAL :: nskip
935 :
936 : INTEGER :: i, iprint, ishift, it1, j, nsamples, &
937 : nstep
938 : REAL(KIND=dp) :: energy_check, old_epx, old_epz, t1
939 : TYPE(meta_env_type), POINTER :: meta_env_saved
940 :
941 4 : IF (PRESENT(nskip)) THEN
942 2 : nsamples = nskip
943 2 : ishift = nskip
944 : ELSE
945 4 : nsamples = 0
946 4 : fz = 0.0_dp
947 : ishift = 0
948 : END IF
949 4 : CALL get_mc_par(mc_par, nstep=nstep, iprint=iprint)
950 4 : meta_env_saved => force_env%meta_env
951 4 : NULLIFY (force_env%meta_env)
952 4 : CALL force_env_get(force_env, potential_energy=old_epx)
953 4 : force_env%meta_env => meta_env_saved
954 :
955 4 : old_epz = force_env%meta_env%epot_s
956 : !!! average energy will be wrong on restarts
957 4 : averages%ave_energy = 0.0_dp
958 4 : t1 = force_env%qs_env%sim_time
959 4 : it1 = force_env%qs_env%sim_step
960 4 : IF (output_unit > 0) THEN
961 2 : WRITE (output_unit, '(a,l4)') "HMC|restart? ", force_env%meta_env%restart
962 : WRITE (output_unit, '(a,3(f16.8,1x))') &
963 2 : "HMC|Ep, Epx, Epz ", old_epx + force_env%meta_env%epot_s, old_epx, force_env%meta_env%epot_s
964 2 : WRITE (output_unit, '(a)') "#HMC| No | z.. | theta.. | ff_z... | ff_z/n |"
965 : END IF
966 12 : DO i = 1, nstep
967 8 : IF (MOD(i, iprint) == 0 .AND. (output_unit > 0)) THEN
968 4 : WRITE (output_unit, '(a,1x,i0)') "HMC|========== On Monte Carlo cycle ", i + ishift
969 4 : WRITE (output_unit, '(a)') "HMC| Attempting a minitrajectory move"
970 4 : WRITE (output_unit, '(a,1x,i0)') "HMC| start mini-trajectory", i + ishift
971 4 : WRITE (output_unit, '(a,1x,i0,1x)', advance="no") "#HMC|0 ", i + ishift
972 8 : DO j = 1, force_env%meta_env%n_colvar
973 4 : WRITE (output_unit, '(f16.8,1x,f16.8,1x,f16.8)', advance="no") force_env%meta_env%metavar(j)%ss0, &
974 4 : force_env%meta_env%metavar(j)%ss, &
975 12 : force_env%meta_env%metavar(j)%ff_s !,fz(j)/real(i+ishift,dp)
976 : END DO
977 4 : WRITE (output_unit, *)
978 : END IF
979 :
980 : CALL mc_hmc_move(mc_par, force_env, globenv, moves, gmoves, old_epx, old_epz, energy_check, &
981 8 : r, output_unit, rng_stream_mc, zbuff)
982 : ! check averages...
983 : ! force average for z needed too...
984 : averages%ave_energy = averages%ave_energy*REAL(i - 1, dp)/REAL(i, dp) + &
985 8 : old_epx/REAL(i, dp)
986 16 : DO j = 1, force_env%meta_env%n_colvar
987 16 : fz(j) = fz(j) + force_env%meta_env%metavar(j)%ff_s
988 : END DO
989 8 : IF (output_unit > 0) THEN
990 4 : WRITE (output_unit, '(a,1x,i0)') "HMC|end mini-trajectory", i + ishift
991 : !!!!!!!! this prints z and theta(x) --ss0,ss-- needed to determine an acceptable k then
992 : ! the instanteneous force and some instanteneous average for force
993 4 : WRITE (output_unit, '(a,1x,i0,1x)', advance="no") "#HMC|1 ", i + ishift
994 8 : DO j = 1, force_env%meta_env%n_colvar
995 4 : WRITE (output_unit, '(f16.8,1x,f16.8,1x,f16.8,1x,f16.8)', advance="no") force_env%meta_env%metavar(j)%ss0, &
996 4 : force_env%meta_env%metavar(j)%ss, &
997 12 : force_env%meta_env%metavar(j)%ff_s, fz(j)/REAL(i + ishift, dp)
998 : END DO
999 4 : WRITE (output_unit, *)
1000 : END IF
1001 8 : nsamples = nsamples + 1
1002 12 : IF (MOD(i, iprint) == 0 .AND. (output_unit > 0)) THEN
1003 4 : WRITE (output_unit, '(a,f16.8)') "HMC| Running average for potential energy ", averages%ave_energy
1004 4 : WRITE (output_unit, '(a,1x,i0)') "HMC|======== End Monte Carlo cycle ", i + ishift
1005 : END IF
1006 : END DO
1007 4 : force_env%qs_env%sim_time = t1
1008 4 : force_env%qs_env%sim_step = it1
1009 4 : IF (output_unit > 0) THEN
1010 2 : WRITE (output_unit, '(a,i0,a,i0,a,f16.8)') "HMC| local acceptance ratio: ", moves%hmc%successes, "/", &
1011 4 : moves%hmc%attempts, "=", REAL(moves%hmc%successes, dp)/REAL(moves%hmc%attempts, dp)
1012 2 : WRITE (output_unit, '(a,i0,a,i0,a,f16.8)') "HMC| global acceptance ratio: ", gmoves%hmc%successes, "/", &
1013 4 : gmoves%hmc%attempts, "=", REAL(gmoves%hmc%successes, dp)/REAL(gmoves%hmc%attempts, dp)
1014 : END IF
1015 : !average force
1016 8 : DO j = 1, force_env%meta_env%n_colvar
1017 8 : force_env%meta_env%metavar(j)%ff_s = fz(j)/nsamples
1018 : END DO
1019 4 : END SUBROUTINE HMCsampler
1020 :
1021 : ! **************************************************************************************************
1022 : !> \brief performs a hybrid Monte Carlo move
1023 : !> \param mc_par ...
1024 : !> \param force_env ...
1025 : !> \param globenv ...
1026 : !> \param moves ...
1027 : !> \param gmoves ...
1028 : !> \param old_epx ...
1029 : !> \param old_epz ...
1030 : !> \param energy_check ...
1031 : !> \param r ...
1032 : !> \param output_unit ...
1033 : !> \param rng_stream ...
1034 : !> \param zbuff ...
1035 : !> \author Alin M Elena
1036 : !> \note It runs a NVE trajectory, followed by localisation and accepts rejects
1037 : !> using the biased Hamiltonian, rather than the traditional guiding Hamiltonian
1038 : ! **************************************************************************************************
1039 32 : SUBROUTINE mc_hmc_move(mc_par, force_env, globenv, moves, gmoves, old_epx, old_epz, &
1040 8 : energy_check, r, output_unit, rng_stream, zbuff)
1041 :
1042 : TYPE(mc_simpar_type), POINTER :: mc_par
1043 : TYPE(force_env_type), POINTER :: force_env
1044 : TYPE(global_environment_type), POINTER :: globenv
1045 : TYPE(mc_moves_type), POINTER :: moves, gmoves
1046 : REAL(KIND=dp), INTENT(INOUT) :: old_epx, old_epz, energy_check
1047 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: r
1048 : INTEGER, INTENT(IN) :: output_unit
1049 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
1050 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: zbuff
1051 :
1052 : CHARACTER(LEN=*), PARAMETER :: routineN = 'mc_hmc_move'
1053 :
1054 : INTEGER :: handle, iatom, j, nAtoms, source
1055 : LOGICAL :: ionode, localise
1056 : REAL(KIND=dp) :: BETA, energy_term, exp_max_val, &
1057 : exp_min_val, new_energy, new_epx, &
1058 : new_epz, rand, value, w
1059 : TYPE(cp_subsys_type), POINTER :: oldsys
1060 : TYPE(mc_ekin_type), POINTER :: hmc_ekin
1061 : TYPE(meta_env_type), POINTER :: meta_env_saved
1062 : TYPE(mp_comm_type) :: group
1063 : TYPE(particle_list_type), POINTER :: particles_set
1064 : TYPE(section_vals_type), POINTER :: dft_section, input
1065 :
1066 : ! begin the timing of the subroutine
1067 :
1068 8 : CALL timeset(routineN, handle)
1069 :
1070 8 : CALL get_qs_env(force_env%qs_env, input=input)
1071 8 : dft_section => section_vals_get_subs_vals(input, "DFT")
1072 :
1073 : ! get a bunch of stuff from mc_par
1074 : CALL get_mc_par(mc_par, ionode=ionode, &
1075 : BETA=BETA, exp_max_val=exp_max_val, &
1076 8 : exp_min_val=exp_min_val, source=source, group=group)
1077 :
1078 : ! nullify some pointers
1079 : ! NULLIFY(particles_set,oldsys,hmc_ekin)
1080 8 : NULLIFY (particles_set, oldsys, meta_env_saved, hmc_ekin)
1081 : ! now let's grab the particle positions
1082 8 : CALL force_env_get(force_env, subsys=oldsys)
1083 8 : CALL cp_subsys_get(oldsys, particles=particles_set)
1084 8 : nAtoms = SIZE(particles_set%els)
1085 : ! do some allocation
1086 :
1087 8 : ALLOCATE (hmc_ekin)
1088 :
1089 : ! record the attempt
1090 8 : moves%hmc%attempts = moves%hmc%attempts + 1
1091 8 : gmoves%hmc%attempts = gmoves%hmc%attempts + 1
1092 :
1093 : ! save the old coordinates just in case we need to go back
1094 56 : DO iatom = 1, nAtoms
1095 200 : r(1:3, iatom) = particles_set%els(iatom)%r(1:3)
1096 : END DO
1097 8 : localise = .TRUE.
1098 : ! the same for collective variables data should be made,ss first half and ff_s the last half
1099 16 : DO j = 1, force_env%meta_env%n_colvar
1100 8 : zbuff(j) = force_env%meta_env%metavar(j)%ss
1101 8 : zbuff(j + force_env%meta_env%n_colvar) = force_env%meta_env%metavar(j)%ff_s
1102 8 : IF ((oldsys%colvar_p(force_env%meta_env%metavar(j)%icolvar)%colvar%type_id == HBP_colvar_id) .OR. &
1103 8 : (oldsys%colvar_p(force_env%meta_env%metavar(j)%icolvar)%colvar%type_id == WC_colvar_id)) THEN
1104 8 : localise = .FALSE.
1105 : END IF
1106 : END DO
1107 :
1108 : ! now run the MD simulation
1109 8 : meta_env_saved => force_env%meta_env
1110 8 : NULLIFY (force_env%meta_env)
1111 8 : force_env%qs_env%sim_time = 0.0_dp
1112 8 : force_env%qs_env%sim_step = 0
1113 8 : IF (.NOT. localise) THEN
1114 8 : CALL section_vals_val_set(dft_section, "LOCALIZE%_SECTION_PARAMETERS_", l_val=.FALSE.)
1115 : END IF
1116 8 : CALL qs_mol_dyn(force_env, globenv, hmc_e_initial=hmc_ekin%initial_ekin, hmc_e_final=hmc_ekin%final_ekin)
1117 8 : IF (.NOT. localise) THEN
1118 8 : CALL section_vals_val_set(dft_section, "LOCALIZE%_SECTION_PARAMETERS_", l_val=.TRUE.)
1119 8 : CALL scf_post_calculation_gpw(force_env%qs_env)
1120 : END IF
1121 :
1122 8 : CALL force_env_get(force_env, potential_energy=new_epx)
1123 :
1124 8 : force_env%meta_env => meta_env_saved
1125 8 : CALL tamc_force(force_env, zpot=new_epz)
1126 8 : new_energy = new_epx + new_epz
1127 8 : IF (output_unit > 0) THEN
1128 : WRITE (output_unit, '(a,4(f16.8,1x))') &
1129 4 : "HMC|old Ep, Ekx, Epz, Epx ", old_epx + old_epz, hmc_ekin%initial_ekin, old_epz, old_epx
1130 4 : WRITE (output_unit, '(a,4(f16.8,1x))') "HMC|new Ep, Ekx, Epz, Epx ", new_energy, hmc_ekin%final_ekin, new_epz, new_epx
1131 : END IF
1132 8 : energy_term = new_energy - old_epx - old_epz + hmc_ekin%final_ekin - hmc_ekin%initial_ekin
1133 :
1134 8 : value = -BETA*(energy_term)
1135 : ! to prevent overflows
1136 8 : IF (value > exp_max_val) THEN
1137 : w = 10.0_dp
1138 8 : ELSE IF (value < exp_min_val) THEN
1139 : w = 0.0_dp
1140 : ELSE
1141 0 : w = EXP(value)
1142 : END IF
1143 :
1144 8 : rand = rng_stream%next()
1145 8 : IF (rand < w) THEN
1146 : ! accept the move
1147 0 : moves%hmc%successes = moves%hmc%successes + 1
1148 0 : gmoves%hmc%successes = gmoves%hmc%successes + 1
1149 : ! update energies
1150 0 : energy_check = energy_check + (new_energy - old_epx - old_epz)
1151 0 : old_epx = new_epx
1152 0 : old_epz = new_epz
1153 : ELSE
1154 : ! reset the cell and particle positions
1155 56 : DO iatom = 1, nAtoms
1156 200 : particles_set%els(iatom)%r(1:3) = r(1:3, iatom)
1157 : END DO
1158 16 : DO j = 1, force_env%meta_env%n_colvar
1159 8 : force_env%meta_env%metavar(j)%ss = zbuff(j)
1160 16 : force_env%meta_env%metavar(j)%ff_s = zbuff(j + force_env%meta_env%n_colvar)
1161 : END DO
1162 :
1163 : END IF
1164 :
1165 8 : DEALLOCATE (hmc_ekin)
1166 :
1167 : ! end the timing
1168 8 : CALL timestop(handle)
1169 :
1170 8 : END SUBROUTINE mc_hmc_move
1171 :
1172 : ! **************************************************************************************************
1173 : !> \brief ...
1174 : !> \param force_env ...
1175 : ! **************************************************************************************************
1176 4 : SUBROUTINE metadyn_write_colvar_header(force_env)
1177 : TYPE(force_env_type), POINTER :: force_env
1178 :
1179 : CHARACTER(len=*), PARAMETER :: routineN = 'metadyn_write_colvar_header'
1180 :
1181 : CHARACTER(len=100) :: aux, fmt
1182 : CHARACTER(len=255) :: label1, label2, label3, label4, label5, &
1183 : label6
1184 : INTEGER :: handle, i, iw, m
1185 : TYPE(cp_logger_type), POINTER :: logger
1186 : TYPE(meta_env_type), POINTER :: meta_env
1187 :
1188 2 : NULLIFY (logger, meta_env)
1189 2 : meta_env => force_env%meta_env
1190 2 : IF (.NOT. ASSOCIATED(meta_env)) RETURN
1191 :
1192 2 : CALL timeset(routineN, handle)
1193 2 : logger => cp_get_default_logger()
1194 :
1195 : iw = cp_print_key_unit_nr(logger, meta_env%metadyn_section, &
1196 2 : "PRINT%COLVAR", extension=".metadynLog")
1197 2 : IF (iw > 0) THEN
1198 1 : label1 = ""
1199 1 : label2 = ""
1200 1 : label3 = ""
1201 1 : label4 = ""
1202 1 : label5 = ""
1203 1 : label6 = ""
1204 2 : DO i = 1, meta_env%n_colvar
1205 1 : WRITE (aux, '(a,i0)') "z_", i
1206 1 : label1 = TRIM(label1)//TRIM(aux)
1207 1 : m = 15*i - LEN_TRIM(label1) - 1
1208 12 : label1 = TRIM(label1)//REPEAT(" ", m)//"|"
1209 1 : WRITE (aux, '(a,i0)') "Theta_", i
1210 1 : label2 = TRIM(label2)//TRIM(aux)
1211 1 : m = 15*i - LEN_TRIM(label2) - 1
1212 8 : label2 = TRIM(label2)//REPEAT(" ", m)//"|"
1213 1 : WRITE (aux, '(a,i0)') "F_z", i
1214 1 : label3 = TRIM(label3)//TRIM(aux)
1215 1 : m = 15*i - LEN_TRIM(label3) - 1
1216 11 : label3 = TRIM(label3)//REPEAT(" ", m)//"|"
1217 1 : WRITE (aux, '(a,i0)') "F_h", i
1218 1 : label4 = TRIM(label4)//TRIM(aux)
1219 1 : m = 15*i - LEN_TRIM(label4) - 1
1220 11 : label4 = TRIM(label4)//REPEAT(" ", m)//"|"
1221 1 : WRITE (aux, '(a,i0)') "F_w", i
1222 1 : label5 = TRIM(label5)//TRIM(aux)
1223 1 : m = 15*i - LEN_TRIM(label5) - 1
1224 11 : label5 = TRIM(label5)//REPEAT(" ", m)//"|"
1225 1 : WRITE (aux, '(a,i0)') "v_z", i
1226 1 : label6 = TRIM(label6)//TRIM(aux)
1227 1 : m = 15*i - LEN_TRIM(label6) - 1
1228 12 : label6 = TRIM(label6)//REPEAT(" ", m)//"|"
1229 : END DO
1230 1 : WRITE (fmt, '("(a17,6a",i0 ,",4a15)")') meta_env%n_colvar*15
1231 1 : WRITE (iw, TRIM(fmt)) "#Time[fs] |", &
1232 1 : TRIM(label1), &
1233 1 : TRIM(label2), &
1234 1 : TRIM(label3), &
1235 1 : TRIM(label4), &
1236 1 : TRIM(label5), &
1237 1 : TRIM(label6), &
1238 1 : "Epot_z |", &
1239 1 : "Ene hills |", &
1240 1 : "Epot walls |", &
1241 2 : "Temperature |"
1242 :
1243 : END IF
1244 : CALL cp_print_key_finished_output(iw, logger, meta_env%metadyn_section, &
1245 2 : "PRINT%COLVAR")
1246 :
1247 2 : CALL timestop(handle)
1248 :
1249 : END SUBROUTINE metadyn_write_colvar_header
1250 :
1251 : ! **************************************************************************************************
1252 : !> \brief ...
1253 : !> \param force_env ...
1254 : ! **************************************************************************************************
1255 8 : SUBROUTINE metadyn_write_colvar(force_env)
1256 : TYPE(force_env_type), POINTER :: force_env
1257 :
1258 : CHARACTER(len=*), PARAMETER :: routineN = 'metadyn_write_colvar'
1259 :
1260 : INTEGER :: handle, i, i_c, iw
1261 : REAL(KIND=dp) :: temp
1262 : TYPE(cp_logger_type), POINTER :: logger
1263 : TYPE(meta_env_type), POINTER :: meta_env
1264 : TYPE(metavar_type), POINTER :: cv
1265 :
1266 4 : NULLIFY (logger, meta_env, cv)
1267 4 : meta_env => force_env%meta_env
1268 4 : IF (.NOT. ASSOCIATED(meta_env)) RETURN
1269 :
1270 4 : CALL timeset(routineN, handle)
1271 4 : logger => cp_get_default_logger()
1272 :
1273 4 : IF (meta_env%langevin) THEN
1274 4 : meta_env%ekin_s = 0.0_dp
1275 : ! meta_env%epot_s = 0.0_dp
1276 8 : DO i_c = 1, meta_env%n_colvar
1277 4 : cv => meta_env%metavar(i_c)
1278 8 : meta_env%ekin_s = meta_env%ekin_s + 0.5_dp*cv%mass*cv%vvp**2
1279 : END DO
1280 : END IF
1281 :
1282 : ! write COLVAR file
1283 : iw = cp_print_key_unit_nr(logger, meta_env%metadyn_section, &
1284 4 : "PRINT%COLVAR", extension=".metadynLog")
1285 4 : IF (iw > 0) THEN
1286 2 : IF (meta_env%extended_lagrange) THEN
1287 2 : WRITE (iw, '(f16.8,70f15.8)') meta_env%time*femtoseconds, &
1288 6 : (meta_env%metavar(i)%ss0, i=1, meta_env%n_colvar), &
1289 6 : (meta_env%metavar(i)%ss, i=1, meta_env%n_colvar), &
1290 6 : (meta_env%metavar(i)%ff_s, i=1, meta_env%n_colvar), &
1291 6 : (meta_env%metavar(i)%ff_hills, i=1, meta_env%n_colvar), &
1292 6 : (meta_env%metavar(i)%ff_walls, i=1, meta_env%n_colvar), &
1293 6 : (meta_env%metavar(i)%vvp, i=1, meta_env%n_colvar), &
1294 2 : meta_env%epot_s, &
1295 2 : meta_env%hills_env%energy, &
1296 2 : meta_env%epot_walls, &
1297 28 : (meta_env%ekin_s)*2.0_dp/(REAL(meta_env%n_colvar, KIND=dp))*kelvin
1298 : ELSE
1299 0 : WRITE (iw, '(f16.8,40f13.5)') meta_env%time*femtoseconds, &
1300 0 : (meta_env%metavar(i)%ss0, i=1, meta_env%n_colvar), &
1301 0 : (meta_env%metavar(i)%ff_hills, i=1, meta_env%n_colvar), &
1302 0 : (meta_env%metavar(i)%ff_walls, i=1, meta_env%n_colvar), &
1303 0 : meta_env%hills_env%energy, &
1304 0 : meta_env%epot_walls
1305 : END IF
1306 : END IF
1307 : CALL cp_print_key_finished_output(iw, logger, meta_env%metadyn_section, &
1308 4 : "PRINT%COLVAR")
1309 : ! Temperature for COLVAR
1310 4 : IF (meta_env%extended_lagrange) THEN
1311 4 : temp = meta_env%ekin_s*2.0_dp/(REAL(meta_env%n_colvar, KIND=dp))*kelvin
1312 : meta_env%avg_temp = (meta_env%avg_temp*REAL(meta_env%n_steps, KIND=dp) + &
1313 4 : temp)/REAL(meta_env%n_steps + 1, KIND=dp)
1314 : iw = cp_print_key_unit_nr(logger, meta_env%metadyn_section, &
1315 4 : "PRINT%TEMPERATURE_COLVAR", extension=".metadynLog")
1316 4 : IF (iw > 0) THEN
1317 2 : WRITE (iw, '(T2,79("-"))')
1318 2 : WRITE (iw, '( A,T51,f10.2,T71,f10.2)') ' COLVARS INSTANTANEOUS/AVERAGE TEMPERATURE ', &
1319 4 : temp, meta_env%avg_temp
1320 2 : WRITE (iw, '(T2,79("-"))')
1321 : END IF
1322 : CALL cp_print_key_finished_output(iw, logger, meta_env%metadyn_section, &
1323 4 : "PRINT%TEMPERATURE_COLVAR")
1324 : END IF
1325 4 : CALL timestop(handle)
1326 :
1327 : END SUBROUTINE metadyn_write_colvar
1328 :
1329 : ! **************************************************************************************************
1330 : !> \brief ...
1331 : !> \param force_env ...
1332 : ! **************************************************************************************************
1333 2 : SUBROUTINE setup_velocities_z(force_env)
1334 : TYPE(force_env_type), POINTER :: force_env
1335 :
1336 : INTEGER :: i_c
1337 : REAL(kind=dp) :: ekin_w, fac_t
1338 : TYPE(meta_env_type), POINTER :: meta_env
1339 : TYPE(metavar_type), POINTER :: cv
1340 :
1341 2 : NULLIFY (meta_env)
1342 2 : meta_env => force_env%meta_env
1343 2 : meta_env%ekin_s = 0.0_dp
1344 4 : DO i_c = 1, meta_env%n_colvar
1345 2 : cv => meta_env%metavar(i_c)
1346 2 : cv%vvp = force_env%globenv%gaussian_rng_stream%next()
1347 4 : meta_env%ekin_s = meta_env%ekin_s + 0.5_dp*cv%mass*cv%vvp**2
1348 : END DO
1349 2 : ekin_w = 0.5_dp*meta_env%temp_wanted*REAL(meta_env%n_colvar, KIND=dp)
1350 2 : fac_t = SQRT(ekin_w/MAX(meta_env%ekin_s, 1.0E-8_dp))
1351 4 : DO i_c = 1, meta_env%n_colvar
1352 2 : cv => meta_env%metavar(i_c)
1353 4 : cv%vvp = cv%vvp*fac_t
1354 : END DO
1355 2 : END SUBROUTINE setup_velocities_z
1356 : END MODULE tamc_run
|
