Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 : ! Initialize velocities possibly applying constraints at the zeroth MD step
273 : ! ! ! CALL section_vals_val_get(motion_section,"PRINT%RESTART%SPLIT_RESTART_FILE",&
274 : ! ! ! l_val=write_binary_restart_file)
275 : !! let us see if this created all the trouble
276 : ! CALL setup_velocities(force_env,simpar,globenv,md_env,md_section,constraint_section, &
277 : ! write_binary_restart_file)
278 :
279 : ! Setup Free Energy Calculation (if required)
280 2 : CALL fe_env_create(fe_env, free_energy_section)
281 : CALL set_md_env(md_env=md_env, simpar=simpar, fe_env=fe_env, cell=cell, &
282 2 : force_env=force_env)
283 :
284 : ! Possibly initialize Wiener processes
285 2 : IF (simpar%ensemble == langevin_ensemble) CALL create_wiener_process(md_env)
286 2 : time_iter_start = m_walltime()
287 :
288 : CALL get_md_env(md_env, force_env=force_env, itimes=itimes, constant=constant, &
289 2 : md_ener=md_ener, t=time, used_time=used_time)
290 :
291 : ! Attach the time counter of the meta_env to the one of the MD
292 2 : CALL set_meta_env(force_env%meta_env, time=time)
293 : ! Initialize the md_ener structure
294 :
295 2 : force_env%meta_env%dt = force_env%meta_env%zdt
296 2 : CALL initialize_md_ener(md_ener, force_env, simpar)
297 : ! force_env%meta_env%dt=force_env%meta_env%zdt
298 :
299 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! MC setup up
300 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
301 :
302 2 : NULLIFY (mc_env, mc_par, MCaverages)
303 :
304 2 : CALL section_vals_get(force_env_section, n_repetition=isos)
305 2 : CPASSERT(isos == 1)
306 : ! set some values...will use get_globenv if that ever comes around
307 :
308 : ! initialize the random numbers
309 : ! IF (para_env%is_source()) THEN
310 : rng_stream_mc = rng_stream_type(name="Random numbers for monte carlo acc/rej", &
311 2 : distribution_type=UNIFORM)
312 : ! ENDIF
313 : !!!!! this should go in a routine hmc_read
314 :
315 2 : NULLIFY (mc_section)
316 2 : ALLOCATE (mc_par)
317 :
318 : mc_section => section_vals_get_subs_vals(force_env%root_section, &
319 2 : "MOTION%MC")
320 : CALL section_vals_val_get(mc_section, "ENSEMBLE", &
321 2 : c_val=ensemble)
322 2 : CPASSERT(str_comp(ensemble, "TRADITIONAL"))
323 : CALL section_vals_val_get(mc_section, "NSTEP", &
324 2 : i_val=nmccycles)
325 2 : CPASSERT(nmccycles > 0)
326 : CALL section_vals_val_get(mc_section, "IPRINT", &
327 2 : i_val=iprint)
328 2 : CALL section_vals_val_get(mc_section, "RANDOMTOSKIP", i_val=rand2skip)
329 2 : CPASSERT(rand2skip >= 0)
330 2 : temp = cp_unit_from_cp2k(simpar%temp_ext, "K")
331 : !
332 :
333 : CALL set_mc_par(mc_par, ensemble=ensemble, nstep=nmccycles, iprint=iprint, temperature=temp, &
334 : beta=1.0_dp/temp/boltzmann*joule, exp_max_val=0.9_dp*LOG(HUGE(0.0_dp)), &
335 : exp_min_val=0.9_dp*LOG(TINY(0.0_dp)), max_val=HUGE(0.0_dp), min_val=0.0_dp, &
336 2 : source=para_env%source, group=para_env, ionode=para_env%is_source(), rand2skip=rand2skip)
337 :
338 2 : output_unit = cp_logger_get_default_io_unit(logger)
339 2 : IF (output_unit > 0) THEN
340 1 : WRITE (output_unit, '(a,a)') "HMC| Hybrid Monte Carlo Scheme "
341 1 : WRITE (output_unit, '(a,a)') "HMC| Ensemble ", ADJUSTL(ensemble)
342 1 : WRITE (output_unit, '(a,i0)') "HMC| MC Cycles ", nmccycles
343 1 : WRITE (output_unit, '(a,i0,a)') "HMC| Print every ", iprint, " cycles"
344 1 : WRITE (output_unit, '(a,i0)') "HMC| Number of random numbers to skip ", rand2skip
345 1 : WRITE (output_unit, '(a,f16.8,a)') "HMC| Temperature ", temp, "K"
346 : END IF
347 :
348 2 : CALL force_env_get(force_env, subsys=subsys)
349 :
350 : CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, &
351 2 : particles=particles, virial=virial)
352 :
353 2 : DO i = 1, rand2skip
354 0 : auxRandom = rng_stream_mc%next()
355 2 : DO j = 1, 3*SIZE(particles%els)
356 0 : auxRandom = globenv%gaussian_rng_stream%next()
357 : END DO
358 : END DO
359 :
360 2 : ALLOCATE (mc_env)
361 2 : CALL mc_env_create(mc_env)
362 2 : CALL set_mc_env(mc_env, mc_par=mc_par, force_env=force_env)
363 : !!!!!!!end mc setup
364 :
365 : ! Check for ensembles requiring the stress tensor - takes into account the possibility for
366 : ! multiple force_evals
367 : IF ((simpar%ensemble == npt_i_ensemble) .OR. &
368 : (simpar%ensemble == npt_ia_ensemble) .OR. &
369 : (simpar%ensemble == npt_f_ensemble) .OR. &
370 : (simpar%ensemble == npe_f_ensemble) .OR. &
371 : (simpar%ensemble == npe_i_ensemble) .OR. &
372 2 : (simpar%ensemble == nph_uniaxial_ensemble) .OR. &
373 : (simpar%ensemble == nph_uniaxial_damped_ensemble)) THEN
374 0 : check = virial%pv_availability
375 0 : IF (.NOT. check) &
376 : CALL cp_abort(__LOCATION__, &
377 : "Virial evaluation not requested for this run in the input file! "// &
378 : "You may consider to switch on the virial evaluation with the keyword: STRESS_TENSOR. "// &
379 0 : "Be sure the method you are using can compute the virial!")
380 0 : IF (ASSOCIATED(force_env%sub_force_env)) THEN
381 0 : DO i = 1, SIZE(force_env%sub_force_env)
382 0 : IF (ASSOCIATED(force_env%sub_force_env(i)%force_env)) THEN
383 0 : CALL force_env_get(force_env%sub_force_env(i)%force_env, subsys=subsys_i)
384 0 : CALL cp_subsys_get(subsys_i, virial=virial)
385 0 : check = check .AND. virial%pv_availability
386 : END IF
387 : END DO
388 : END IF
389 0 : IF (.NOT. check) &
390 : CALL cp_abort(__LOCATION__, &
391 : "Virial evaluation not requested for all the force_eval sections present in"// &
392 : " the input file! You have to switch on the virial evaluation with the keyword: STRESS_TENSOR"// &
393 0 : " in each force_eval section. Be sure the method you are using can compute the virial!")
394 : END IF
395 :
396 : ! Computing Forces at zero MD step
397 2 : IF (simpar%ensemble /= reftraj_ensemble) THEN
398 2 : CALL section_vals_val_get(md_section, "STEP_START_VAL", i_val=itimes)
399 2 : CALL section_vals_val_get(md_section, "TIME_START_VAL", r_val=time)
400 2 : CALL section_vals_val_get(md_section, "ECONS_START_VAL", r_val=constant)
401 2 : CALL section_vals_val_set(md_section, "STEP_START_VAL", i_val=initialStep)
402 2 : CALL section_vals_val_set(md_section, "TIME_START_VAL", r_val=inittime)
403 2 : initialStep = itimes
404 2 : CALL cp_iterate(logger%iter_info, iter_nr=itimes)
405 2 : IF (save_mem) THEN
406 0 : work_section => section_vals_get_subs_vals(subsys_section, "VELOCITY")
407 0 : CALL section_vals_remove_values(work_section)
408 0 : work_section => section_vals_get_subs_vals(subsys_section, "SHELL_VELOCITY")
409 0 : CALL section_vals_remove_values(work_section)
410 0 : work_section => section_vals_get_subs_vals(subsys_section, "CORE_VELOCITY")
411 0 : CALL section_vals_remove_values(work_section)
412 : END IF
413 :
414 : ! CALL force_env_calc_energy_force (force_env, calc_force=.TRUE.)
415 2 : meta_env_saved => force_env%meta_env
416 2 : NULLIFY (force_env%meta_env)
417 2 : CALL force_env_calc_energy_force(force_env, calc_force=.FALSE.)
418 2 : force_env%meta_env => meta_env_saved
419 :
420 2 : IF (ASSOCIATED(force_env%qs_env)) THEN
421 : ! force_env%qs_env%sim_time=time
422 : ! force_env%qs_env%sim_step=itimes
423 2 : force_env%qs_env%sim_time = 0.0_dp
424 2 : force_env%qs_env%sim_step = 0
425 : END IF
426 : ! Warm-up engines for metadynamics
427 2 : IF (ASSOCIATED(force_env%meta_env)) THEN
428 2 : IF (force_env%meta_env%langevin) THEN
429 2 : CALL create_wiener_process_cv(force_env%meta_env)
430 2 : DO j = 1, (rand2skip - 1)/nmccycles
431 2 : DO i = 1, force_env%meta_env%n_colvar
432 0 : auxRandom = force_env%meta_env%rng(i)%next()
433 0 : auxRandom = force_env%meta_env%rng(i)%next()
434 : END DO
435 : END DO
436 : END IF
437 : ! IF (force_env%meta_env%well_tempered) THEN
438 : ! force_env%meta_env%wttemperature = simpar%temp_ext
439 : ! IF (force_env%meta_env%wtgamma>EPSILON(1._dp)) THEN
440 : ! dummy=force_env%meta_env%wttemperature*(force_env%meta_env%wtgamma-1._dp)
441 : ! IF (force_env%meta_env%delta_t>EPSILON(1._dp)) THEN
442 : ! check=ABS(force_env%meta_env%delta_t-dummy)<1.E+3_dp*EPSILON(1._dp)
443 : ! IF(.NOT.check) CALL cp_abort(__LOCATION__,&
444 : ! "Inconsistency between DELTA_T and WTGAMMA (both specified):"//&
445 : ! " please, verify that DELTA_T=(WTGAMMA-1)*TEMPERATURE")
446 : ! ELSE
447 : ! force_env%meta_env%delta_t = dummy
448 : ! ENDIF
449 : ! ELSE
450 : ! force_env%meta_env%wtgamma = 1._dp &
451 : ! + force_env%meta_env%delta_t/force_env%meta_env%wttemperature
452 : ! ENDIF
453 : ! force_env%meta_env%invdt = 1._dp/force_env%meta_env%delta_t
454 : ! ENDIF
455 2 : CALL tamc_force(force_env)
456 : ! CALL metadyn_write_colvar(force_env)
457 : END IF
458 :
459 2 : IF (simpar%do_respa) THEN
460 : CALL force_env_calc_energy_force(force_env%sub_force_env(1)%force_env, &
461 0 : calc_force=.TRUE.)
462 : END IF
463 :
464 : ! CALL force_env_get( force_env, subsys=subsys)
465 : !
466 : ! CALL cp_subsys_get(subsys,atomic_kinds=atomic_kinds,local_particles=local_particles,&
467 : ! particles=particles)
468 :
469 : CALL virial_evaluate(atomic_kinds%els, particles%els, local_particles, &
470 2 : virial, force_env%para_env)
471 :
472 2 : CALL md_energy(md_env, md_ener)
473 : ! CALL md_write_output(md_env) !inits the print env at itimes == 0 also writes trajectories
474 2 : md_stride = 1
475 : ELSE
476 0 : CALL get_md_env(md_env, reftraj=reftraj)
477 0 : CALL initialize_reftraj(reftraj, reftraj_section, md_env)
478 0 : itimes = reftraj%info%first_snapshot - 1
479 0 : md_stride = reftraj%info%stride
480 : END IF
481 :
482 : CALL cp_print_key_finished_output(simpar%info_constraint, logger, &
483 2 : constraint_section, "CONSTRAINT_INFO")
484 : CALL cp_print_key_finished_output(simpar%lagrange_multipliers, logger, &
485 2 : constraint_section, "LAGRANGE_MULTIPLIERS")
486 2 : CALL init_mc_moves(moves)
487 2 : CALL init_mc_moves(gmoves)
488 6 : ALLOCATE (r(1:3, SIZE(particles%els)))
489 : ! ALLOCATE (r_old(1:3,size(particles%els)))
490 2 : CALL mc_averages_create(MCaverages)
491 : !!!!! some more buffers
492 : ! Allocate random number for Langevin Thermostat acting on COLVARS
493 6 : ALLOCATE (xieta(2*force_env%meta_env%n_colvar))
494 6 : xieta(:) = 0.0_dp
495 6 : ALLOCATE (An(force_env%meta_env%n_colvar))
496 4 : An(:) = 0.0_dp
497 4 : ALLOCATE (fz(force_env%meta_env%n_colvar))
498 4 : fz(:) = 0.0_dp
499 4 : ALLOCATE (zbuff(2*force_env%meta_env%n_colvar))
500 6 : zbuff(:) = 0.0_dp
501 :
502 2 : IF (output_unit > 0) THEN
503 1 : WRITE (output_unit, '(a)') "HMC|==== Initial average forces"
504 : END IF
505 2 : CALL metadyn_write_colvar_header(force_env)
506 2 : moves%hmc%attempts = 0
507 2 : moves%hmc%successes = 0
508 2 : gmoves%hmc%attempts = 0
509 2 : gmoves%hmc%successes = 0
510 2 : IF (initialStep == 0) THEN
511 : !!! if we come from a restart we shall properly compute the average force
512 : !!! read the average force up to now
513 4 : DO i = 1, force_env%meta_env%n_colvar
514 2 : fs_section => section_vals_get_subs_vals(force_env%meta_env%metadyn_section, "EXT_LAGRANGE_FS")
515 2 : CALL section_vals_get(fs_section, explicit=explicit)
516 4 : IF (explicit) THEN
517 : CALL section_vals_val_get(fs_section, "_DEFAULT_KEYWORD_", &
518 0 : i_rep_val=i, r_val=rval)
519 0 : fz(i) = rval*rand2skip
520 : END IF
521 : END DO
522 :
523 : CALL HMCsampler(globenv, force_env, MCaverages, r, mc_par, moves, gmoves, rng_stream_mc, output_unit, &
524 2 : fz, zbuff, nskip=rand2skip)
525 2 : CALL cp_iterate(logger%iter_info, last=.FALSE., iter_nr=0)
526 2 : CALL section_vals_val_set(mc_section, "RANDOMTOSKIP", i_val=rand2skip + nmccycles)
527 2 : CALL write_restart(md_env=md_env, root_section=force_env%root_section)
528 : END IF
529 2 : IF (output_unit > 0) THEN
530 1 : WRITE (output_unit, '(a)') "HMC|==== end initial average forces"
531 : END IF
532 : ! call set_md_env(md_env, init=.FALSE.)
533 :
534 2 : CALL metadyn_write_colvar(force_env)
535 :
536 4 : DO istep = 1, force_env%meta_env%TAMCSteps
537 : ! Increase counters
538 2 : itimes = itimes + 1
539 2 : time = time + force_env%meta_env%dt
540 2 : IF (output_unit > 0) THEN
541 1 : WRITE (output_unit, '(a)') "HMC|==================================="
542 1 : WRITE (output_unit, '(a,1x,i0)') "HMC| on z step ", istep
543 : END IF
544 : !needed when electric field fields are applied
545 2 : IF (ASSOCIATED(force_env%qs_env)) THEN
546 2 : force_env%qs_env%sim_time = time
547 2 : force_env%qs_env%sim_step = itimes
548 2 : force_env%meta_env%time = force_env%qs_env%sim_time
549 : END IF
550 :
551 2 : CALL cp_iterate(logger%iter_info, last=(istep == force_env%meta_env%TAMCSteps), iter_nr=itimes)
552 : ! Open possible Shake output units
553 : simpar%info_constraint = cp_print_key_unit_nr(logger, constraint_section, "CONSTRAINT_INFO", &
554 2 : extension=".shakeLog", log_filename=.FALSE.)
555 : simpar%lagrange_multipliers = cp_print_key_unit_nr( &
556 : logger, constraint_section, &
557 2 : "LAGRANGE_MULTIPLIERS", extension=".LagrangeMultLog", log_filename=.FALSE.)
558 : simpar%dump_lm = BTEST(cp_print_key_should_output(logger%iter_info, constraint_section, &
559 2 : "LAGRANGE_MULTIPLIERS"), cp_p_file)
560 :
561 : ! Velocity Verlet Integrator
562 :
563 2 : moves%hmc%attempts = 0
564 2 : moves%hmc%successes = 0
565 : CALL langevinVEC(md_env, globenv, mc_env, moves, gmoves, r, &
566 2 : rng_stream_mc, xieta, An, fz, MCaverages, zbuff)
567 :
568 : ! Close Shake output if requested...
569 : CALL cp_print_key_finished_output(simpar%info_constraint, logger, &
570 2 : constraint_section, "CONSTRAINT_INFO")
571 : CALL cp_print_key_finished_output(simpar%lagrange_multipliers, logger, &
572 2 : constraint_section, "LAGRANGE_MULTIPLIERS")
573 2 : CALL cp_iterate(logger%iter_info, iter_nr=initialStep)
574 2 : CALL metadyn_write_colvar(force_env)
575 : ! Free Energy calculation
576 : ! CALL free_energy_evaluate(md_env,should_stop,free_energy_section)
577 :
578 : ![AME:UB] IF (should_stop) EXIT
579 :
580 : ! Test for <PROJECT_NAME>.EXIT_MD or for WALL_TIME to exit
581 : ! Default:
582 : ! IF so we don't overwrite the restart or append to the trajectory
583 : ! because the execution could in principle stop inside the SCF where energy
584 : ! and forces are not converged.
585 : ! But:
586 : ! You can force to print the last step (for example if the method used
587 : ! to compute energy and forces is not SCF based) activating the print_key
588 : ! MOTION%MD%PRINT%FORCE_LAST.
589 2 : CALL external_control(should_stop, "MD", globenv=globenv)
590 2 : IF (should_stop) THEN
591 0 : CALL cp_iterate(logger%iter_info, last=.TRUE., iter_nr=itimes)
592 : ! CALL md_output(md_env,md_section,force_env%root_section,should_stop)
593 0 : EXIT
594 : END IF
595 :
596 : ! IF(simpar%ensemble /= reftraj_ensemble) THEN
597 : ! CALL md_energy(md_env, md_ener)
598 : ! CALL temperature_control(simpar, md_env, md_ener, force_env, logger)
599 : ! CALL comvel_control(md_ener, force_env, md_section, logger)
600 : ! CALL angvel_control(md_ener, force_env, md_section, logger)
601 : ! ELSE
602 : ! CALL md_ener_reftraj(md_env, md_ener)
603 : ! END IF
604 :
605 2 : time_iter_stop = m_walltime()
606 2 : used_time = time_iter_stop - time_iter_start
607 2 : time_iter_start = time_iter_stop
608 :
609 : !!!!! this writes the restart...
610 : ! CALL md_output(md_env,md_section,force_env%root_section,should_stop)
611 :
612 : ! IF(simpar%ensemble == reftraj_ensemble ) THEN
613 : ! CALL write_output_reftraj(md_env)
614 : ! END IF
615 :
616 6 : IF (output_unit > 0) THEN
617 1 : WRITE (output_unit, '(a,1x,i0)') "HMC| end z step ", istep
618 1 : WRITE (output_unit, '(a)') "HMC|==================================="
619 : END IF
620 : END DO
621 2 : CALL cp_iterate(logger%iter_info, last=.TRUE., iter_nr=itimes)
622 2 : force_env%qs_env%sim_time = 0.0_dp
623 2 : force_env%qs_env%sim_step = 0
624 2 : rand2skip = rand2skip + nmccycles*force_env%meta_env%TAMCSteps
625 2 : IF (initialStep == 0) rand2skip = rand2skip + nmccycles
626 2 : CALL section_vals_val_set(mc_section, "RANDOMTOSKIP", i_val=rand2skip)
627 :
628 2 : CALL write_restart(md_env=md_env, root_section=force_env%root_section)
629 : ! if we need the final kinetic energy for Hybrid Monte Carlo
630 : ! hmc_ekin%final_ekin=md_ener%ekin
631 :
632 : ! Remove the iteration level
633 2 : CALL cp_rm_iter_level(logger%iter_info, "MD")
634 :
635 : ! Deallocate Thermostats and Barostats
636 2 : CALL release_simpar_type(simpar)
637 :
638 2 : CALL md_env_release(md_env)
639 2 : DEALLOCATE (md_env)
640 : ! Clean restartable sections..
641 2 : IF (my_rm_restart_info) CALL remove_restart_info(force_env%root_section)
642 : ! IF (para_env%is_source()) THEN
643 : ! ENDIF
644 2 : CALL MC_ENV_RELEASE(mc_env)
645 2 : DEALLOCATE (mc_env)
646 2 : DEALLOCATE (mc_par)
647 2 : CALL MC_MOVES_RELEASE(moves)
648 2 : CALL MC_MOVES_RELEASE(gmoves)
649 2 : DEALLOCATE (r)
650 : ! DEALLOCATE(r_old)
651 2 : DEALLOCATE (xieta)
652 2 : DEALLOCATE (An)
653 2 : DEALLOCATE (fz)
654 2 : DEALLOCATE (zbuff)
655 2 : CALL mc_averages_release(MCaverages)
656 2 : CALL timestop(handle)
657 :
658 64 : END SUBROUTINE qs_tamc
659 :
660 : ! **************************************************************************************************
661 : !> \brief Propagates velocities for z half a step
662 : !> \param force_env ...
663 : !> \param An ...
664 : !> \author Alin M Elena
665 : !> \note Vanden-Eijnden Ciccotti C.Phys.Letter 429 (2006) 310-316
666 : ! **************************************************************************************************
667 4 : SUBROUTINE tamc_velocities_colvar(force_env, An)
668 : TYPE(force_env_type), POINTER :: force_env
669 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: An
670 :
671 : CHARACTER(len=*), PARAMETER :: routineN = 'tamc_velocities_colvar'
672 :
673 : INTEGER :: handle, i_c
674 : REAL(kind=dp) :: dt, fft, sigma
675 : TYPE(cp_logger_type), POINTER :: logger
676 : TYPE(meta_env_type), POINTER :: meta_env
677 : TYPE(metavar_type), POINTER :: cv
678 :
679 4 : NULLIFY (logger, meta_env, cv)
680 4 : meta_env => force_env%meta_env
681 4 : CALL timeset(routineN, handle)
682 4 : logger => cp_get_default_logger()
683 : ! Add citation
684 4 : IF (meta_env%langevin) CALL cite_reference(VandenCic2006)
685 4 : dt = meta_env%dt
686 :
687 : ! Evolve Velocities
688 4 : meta_env%epot_walls = 0.0_dp
689 8 : DO i_c = 1, meta_env%n_colvar
690 4 : cv => meta_env%metavar(i_c)
691 4 : fft = cv%ff_s + cv%ff_hills
692 4 : sigma = SQRT((meta_env%temp_wanted*kelvin)*2.0_dp*(boltzmann/joule)*cv%gamma/cv%mass)
693 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)
694 8 : meta_env%epot_walls = meta_env%epot_walls + cv%epot_walls
695 : END DO
696 4 : CALL timestop(handle)
697 4 : END SUBROUTINE tamc_velocities_colvar
698 :
699 : ! **************************************************************************************************
700 : !> \brief propagates z one step
701 : !> \param force_env ...
702 : !> \param xieta ...
703 : !> \author Alin M Elena
704 : !> \note Vanden-Eijnden Ciccotti C.Phys.Letter 429 (2006) 310-316
705 : ! **************************************************************************************************
706 2 : SUBROUTINE tamc_position_colvar(force_env, xieta)
707 : TYPE(force_env_type), POINTER :: force_env
708 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: xieta
709 :
710 : CHARACTER(len=*), PARAMETER :: routineN = 'tamc_position_colvar'
711 :
712 : INTEGER :: handle, i_c
713 : REAL(kind=dp) :: dt, sigma
714 : TYPE(cp_logger_type), POINTER :: logger
715 : TYPE(meta_env_type), POINTER :: meta_env
716 : TYPE(metavar_type), POINTER :: cv
717 :
718 2 : NULLIFY (logger, meta_env, cv)
719 2 : meta_env => force_env%meta_env
720 : ! IF (.NOT.ASSOCIATED(meta_env)) RETURN
721 :
722 2 : CALL timeset(routineN, handle)
723 2 : logger => cp_get_default_logger()
724 :
725 : ! Add citation
726 2 : IF (meta_env%langevin) CALL cite_reference(VandenCic2006)
727 2 : dt = meta_env%dt
728 :
729 : ! Update of ss0
730 4 : DO i_c = 1, meta_env%n_colvar
731 2 : cv => meta_env%metavar(i_c)
732 2 : sigma = SQRT((meta_env%temp_wanted*kelvin)*2.0_dp*(boltzmann/joule)*cv%gamma/cv%mass)
733 : ! cv%ss0 =cv%ss0 +dt*cv%vvp
734 2 : cv%ss0 = cv%ss0 + dt*cv%vvp + dt*SQRT(dt/12.0_dp)*sigma*xieta(i_c + meta_env%n_colvar)
735 4 : IF (cv%periodic) THEN
736 : ! A periodic COLVAR is always within [-pi,pi]
737 0 : cv%ss0 = SIGN(1.0_dp, ASIN(SIN(cv%ss0)))*ACOS(COS(cv%ss0))
738 : END IF
739 : END DO
740 2 : CALL timestop(handle)
741 :
742 2 : END SUBROUTINE tamc_position_colvar
743 :
744 : ! **************************************************************************************************
745 : !> \brief Computes forces on z
746 : !> #details also can be used to get the potenzial evergy of z
747 : !> \param force_env ...
748 : !> \param zpot ...
749 : !> \author Alin M Elena
750 : ! **************************************************************************************************
751 10 : SUBROUTINE tamc_force(force_env, zpot)
752 : TYPE(force_env_type), POINTER :: force_env
753 : REAL(KIND=dp), INTENT(inout), OPTIONAL :: zpot
754 :
755 : CHARACTER(len=*), PARAMETER :: routineN = 'tamc_force'
756 :
757 : INTEGER :: handle, i, i_c, icolvar, ii
758 : LOGICAL :: explicit
759 : REAL(kind=dp) :: diff_ss, dt, rval
760 10 : TYPE(colvar_p_type), DIMENSION(:), POINTER :: colvar_p
761 : TYPE(cp_logger_type), POINTER :: logger
762 : TYPE(cp_subsys_type), POINTER :: subsys
763 : TYPE(meta_env_type), POINTER :: meta_env
764 : TYPE(metavar_type), POINTER :: cv
765 : TYPE(particle_list_type), POINTER :: particles
766 : TYPE(section_vals_type), POINTER :: ss0_section, ss_section, vvp_section
767 :
768 10 : NULLIFY (logger, meta_env)
769 10 : meta_env => force_env%meta_env
770 : ! IF (.NOT.ASSOCIATED(meta_env)) RETURN
771 :
772 10 : CALL timeset(routineN, handle)
773 10 : logger => cp_get_default_logger()
774 10 : NULLIFY (colvar_p, subsys, cv, ss0_section, vvp_section, ss_section)
775 10 : CALL force_env_get(force_env, subsys=subsys)
776 :
777 10 : dt = meta_env%dt
778 10 : IF (.NOT. meta_env%restart) meta_env%n_steps = meta_env%n_steps + 1
779 : ! compute ss and the derivative of ss with respect to the atomic positions
780 20 : DO i_c = 1, meta_env%n_colvar
781 10 : cv => meta_env%metavar(i_c)
782 10 : icolvar = cv%icolvar
783 10 : CALL colvar_eval_glob_f(icolvar, force_env)
784 10 : cv%ss = subsys%colvar_p(icolvar)%colvar%ss
785 : ! Restart for Extended Lagrangian Metadynamics
786 20 : IF (meta_env%restart) THEN
787 : ! Initialize the position of the collective variable in the extended lagrange
788 2 : ss0_section => section_vals_get_subs_vals(meta_env%metadyn_section, "EXT_LAGRANGE_SS0")
789 2 : CALL section_vals_get(ss0_section, explicit=explicit)
790 2 : IF (explicit) THEN
791 : CALL section_vals_val_get(ss0_section, "_DEFAULT_KEYWORD_", &
792 2 : i_rep_val=i_c, r_val=rval)
793 2 : cv%ss0 = rval
794 : ELSE
795 0 : cv%ss0 = cv%ss
796 : END IF
797 2 : vvp_section => section_vals_get_subs_vals(meta_env%metadyn_section, "EXT_LAGRANGE_VVP")
798 2 : CALL section_vals_get(vvp_section, explicit=explicit)
799 2 : IF (explicit) THEN
800 0 : CALL setup_velocities_z(force_env)
801 : CALL section_vals_val_get(vvp_section, "_DEFAULT_KEYWORD_", &
802 0 : i_rep_val=i_c, r_val=rval)
803 0 : cv%vvp = rval
804 : ELSE
805 2 : CALL setup_velocities_z(force_env)
806 : END IF
807 2 : ss_section => section_vals_get_subs_vals(meta_env%metadyn_section, "EXT_LAGRANGE_SS")
808 2 : CALL section_vals_get(ss_section, explicit=explicit)
809 2 : IF (explicit) THEN
810 : CALL section_vals_val_get(ss_section, "_DEFAULT_KEYWORD_", &
811 0 : i_rep_val=i_c, r_val=rval)
812 0 : cv%ss = rval
813 : END IF
814 : END IF
815 : !
816 : END DO
817 : ! forces on the atoms
818 10 : NULLIFY (particles)
819 : CALL cp_subsys_get(subsys, colvar_p=colvar_p, &
820 10 : particles=particles)
821 :
822 10 : meta_env%restart = .FALSE.
823 10 : meta_env%epot_s = 0.0_dp
824 10 : meta_env%epot_walls = 0.0_dp
825 20 : DO i_c = 1, meta_env%n_colvar
826 10 : cv => meta_env%metavar(i_c)
827 10 : diff_ss = cv%ss - cv%ss0
828 10 : IF (cv%periodic) THEN
829 : ! The difference of a periodic COLVAR is always within [-pi,pi]
830 0 : diff_ss = SIGN(1.0_dp, ASIN(SIN(diff_ss)))*ACOS(COS(diff_ss))
831 : END IF
832 10 : cv%epot_s = 0.5_dp*cv%lambda*diff_ss*diff_ss
833 10 : cv%ff_s = cv%lambda*(diff_ss)
834 10 : meta_env%epot_s = meta_env%epot_s + cv%epot_s
835 10 : icolvar = cv%icolvar
836 :
837 50 : DO ii = 1, colvar_p(icolvar)%colvar%n_atom_s
838 30 : i = colvar_p(icolvar)%colvar%i_atom(ii)
839 220 : particles%els(i)%f = particles%els(i)%f - cv%ff_s*colvar_p(icolvar)%colvar%dsdr(:, ii)
840 : END DO
841 :
842 : END DO
843 10 : IF (PRESENT(zpot)) zpot = meta_env%epot_s
844 10 : CALL fix_atom_control(force_env)
845 :
846 10 : CALL timestop(handle)
847 10 : END SUBROUTINE tamc_force
848 :
849 : ! **************************************************************************************************
850 : !> \brief propagates one time step both z systems and samples the x system
851 : !> \param md_env ...
852 : !> \param globenv ...
853 : !> \param mc_env ...
854 : !> \param moves ...
855 : !> \param gmoves ...
856 : !> \param r ...
857 : !> \param rng_stream_mc ...
858 : !> \param xieta ...
859 : !> \param An ...
860 : !> \param fz ...
861 : !> \param averages ...
862 : !> \param zbuff ...
863 : !> \author Alin M Elena
864 : ! **************************************************************************************************
865 8 : SUBROUTINE langevinVEC(md_env, globenv, mc_env, moves, gmoves, r, &
866 2 : rng_stream_mc, xieta, An, fz, averages, zbuff)
867 :
868 : TYPE(md_environment_type), POINTER :: md_env
869 : TYPE(global_environment_type), POINTER :: globenv
870 : TYPE(mc_environment_type), POINTER :: mc_env
871 : TYPE(mc_moves_type), POINTER :: moves, gmoves
872 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: r
873 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream_mc
874 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: xieta, An, fz
875 : TYPE(mc_averages_type), INTENT(INOUT), POINTER :: averages
876 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: zbuff
877 :
878 : INTEGER :: iprint, ivar, nparticle, nparticle_kind, &
879 : nstep, output_unit
880 : REAL(KIND=dp) :: dt, gamma, mass, sigma
881 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
882 2 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
883 : TYPE(cell_type), POINTER :: cell
884 : TYPE(cp_logger_type), POINTER :: logger
885 : TYPE(cp_subsys_type), POINTER :: subsys
886 : TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles
887 : TYPE(force_env_type), POINTER :: force_env
888 : TYPE(global_constraint_type), POINTER :: gci
889 : TYPE(mc_simpar_type), POINTER :: mc_par
890 : TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
891 2 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
892 : TYPE(molecule_list_type), POINTER :: molecules
893 2 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
894 : TYPE(mp_para_env_type), POINTER :: para_env
895 : TYPE(particle_list_type), POINTER :: particles
896 2 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
897 : TYPE(simpar_type), POINTER :: simpar
898 : TYPE(virial_type), POINTER :: virial
899 :
900 2 : NULLIFY (logger, mc_par)
901 4 : logger => cp_get_default_logger()
902 2 : output_unit = cp_logger_get_default_io_unit(logger)
903 :
904 : ! quantitites to be nullified for the get_md_env
905 2 : NULLIFY (simpar, force_env, para_env)
906 : ! quantities to be nullified for the force_env_get environment
907 2 : NULLIFY (subsys, cell)
908 : ! quantitites to be nullified for the cp_subsys_get
909 2 : NULLIFY (atomic_kinds, local_particles, particles, local_molecules, molecules, molecule_kinds, gci)
910 :
911 : CALL get_md_env(md_env=md_env, simpar=simpar, force_env=force_env, &
912 2 : para_env=para_env)
913 2 : CALL get_mc_env(mc_env, mc_par=mc_par)
914 2 : CALL get_mc_par(mc_par, nstep=nstep, iprint=iprint)
915 :
916 2 : dt = simpar%dt
917 2 : CALL force_env_get(force_env=force_env, subsys=subsys, cell=cell)
918 :
919 : !!!! this bit should vanish once I understand what the hell is with it
920 :
921 : ! ! Do some checks on coordinates and box
922 2 : CALL apply_qmmm_walls_reflective(force_env)
923 :
924 : CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, &
925 : particles=particles, local_molecules=local_molecules, molecules=molecules, &
926 2 : molecule_kinds=molecule_kinds, gci=gci, virial=virial)
927 :
928 2 : nparticle_kind = atomic_kinds%n_els
929 2 : atomic_kind_set => atomic_kinds%els
930 2 : molecule_kind_set => molecule_kinds%els
931 :
932 2 : nparticle = particles%n_els
933 2 : particle_set => particles%els
934 2 : molecule_set => molecules%els
935 2 : CPASSERT(ASSOCIATED(force_env%meta_env))
936 2 : CPASSERT(force_env%meta_env%langevin)
937 : ! *** Velocity Verlet for Langevin *** v(t)--> v(t+1/2)
938 : !!!!!! noise xi is in the first half, eta in the second half
939 4 : DO ivar = 1, force_env%meta_env%n_colvar
940 2 : xieta(ivar) = force_env%meta_env%rng(ivar)%next()
941 2 : xieta(ivar + force_env%meta_env%n_colvar) = force_env%meta_env%rng(ivar)%next()
942 2 : gamma = force_env%meta_env%metavar(ivar)%gamma
943 2 : mass = force_env%meta_env%metavar(ivar)%mass
944 2 : sigma = SQRT((force_env%meta_env%temp_wanted*kelvin)*2.0_dp*(boltzmann/joule)*gamma/mass)
945 : An(ivar) = 0.5_dp*SQRT(dt)*sigma*(xieta(ivar)*(1.0_dp - 0.5_dp*dt*gamma) - &
946 4 : dt*gamma*xieta(ivar + force_env%meta_env%n_colvar)/SQRT(12.0_dp))
947 : END DO
948 : ! *** Velocity Verlet for Langeving *** v(t)--> v(t+1/2)
949 2 : CALL tamc_velocities_colvar(force_env, An)
950 : ! *** Velocity Verlet for Langevin S(t)->S(t+1)
951 2 : CALL tamc_position_colvar(force_env, xieta)
952 : !!!!! start zHMC sampler
953 2 : CALL HMCsampler(globenv, force_env, averages, r, mc_par, moves, gmoves, rng_stream_mc, output_unit, fz, zbuff)
954 :
955 : ! CALL final_mc_write(mc_par,tmp_moves,&
956 : ! output_unit,energy_check,&
957 : ! initial_energy,final_energy,&
958 : ! averages)
959 :
960 : !!!!!!!!!!!!!!!!!!!! end zHMC sampler
961 : ! *** Velocity Verlet for Langeving *** v(t+1/2)--> v(t+1)
962 2 : CALL tamc_velocities_colvar(force_env, An)
963 : ! CALL virial_evaluate ( atomic_kind_set, particle_set, &
964 : ! local_particles, virial, para_env)
965 :
966 2 : END SUBROUTINE langevinVEC
967 :
968 : ! **************************************************************************************************
969 : !> \brief Driver routin for the canonical sampler using modified HMC
970 : !> \param globenv ...
971 : !> \param force_env ...
972 : !> \param averages ...
973 : !> \param r ...
974 : !> \param mc_par ...
975 : !> \param moves ...
976 : !> \param gmoves ...
977 : !> \param rng_stream_mc ...
978 : !> \param output_unit ...
979 : !> \param fz ...
980 : !> \param zbuff ...
981 : !> \param nskip ...
982 : !> \author Alin M Elena
983 : !> \note at the end of this routine %ff_s will contain mean force
984 : ! **************************************************************************************************
985 :
986 20 : SUBROUTINE HMCsampler(globenv, force_env, averages, r, mc_par, moves, gmoves, rng_stream_mc, output_unit, &
987 4 : fz, zbuff, nskip)
988 : TYPE(global_environment_type), POINTER :: globenv
989 : TYPE(force_env_type), POINTER :: force_env
990 : TYPE(mc_averages_type), POINTER :: averages
991 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: r
992 : TYPE(mc_simpar_type), POINTER :: mc_par
993 : TYPE(mc_moves_type), POINTER :: moves, gmoves
994 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream_mc
995 : INTEGER, INTENT(IN) :: output_unit
996 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: fz, zbuff
997 : INTEGER, INTENT(IN), OPTIONAL :: nskip
998 :
999 : INTEGER :: i, iprint, ishift, it1, j, nsamples, &
1000 : nstep
1001 : REAL(KIND=dp) :: energy_check, old_epx, old_epz, t1
1002 : TYPE(meta_env_type), POINTER :: meta_env_saved
1003 :
1004 4 : IF (PRESENT(nskip)) THEN
1005 2 : nsamples = nskip
1006 2 : ishift = nskip
1007 : ELSE
1008 4 : nsamples = 0
1009 4 : fz = 0.0_dp
1010 : ishift = 0
1011 : END IF
1012 : ! lrestart = .false.
1013 : ! if (present(logger) .and. present(iter)) THEN
1014 : ! lrestart=.true.
1015 : ! ENDIF
1016 4 : CALL get_mc_par(mc_par, nstep=nstep, iprint=iprint)
1017 4 : meta_env_saved => force_env%meta_env
1018 4 : NULLIFY (force_env%meta_env)
1019 4 : CALL force_env_get(force_env, potential_energy=old_epx)
1020 4 : force_env%meta_env => meta_env_saved
1021 :
1022 4 : old_epz = force_env%meta_env%epot_s
1023 : !!! average energy will be wrong on restarts
1024 4 : averages%ave_energy = 0.0_dp
1025 4 : t1 = force_env%qs_env%sim_time
1026 4 : it1 = force_env%qs_env%sim_step
1027 4 : IF (output_unit > 0) THEN
1028 2 : WRITE (output_unit, '(a,l4)') "HMC|restart? ", force_env%meta_env%restart
1029 : WRITE (output_unit, '(a,3(f16.8,1x))') &
1030 2 : "HMC|Ep, Epx, Epz ", old_epx + force_env%meta_env%epot_s, old_epx, force_env%meta_env%epot_s
1031 2 : WRITE (output_unit, '(a)') "#HMC| No | z.. | theta.. | ff_z... | ff_z/n |"
1032 : END IF
1033 12 : DO i = 1, nstep
1034 8 : IF (MOD(i, iprint) == 0 .AND. (output_unit > 0)) THEN
1035 4 : WRITE (output_unit, '(a,1x,i0)') "HMC|========== On Monte Carlo cycle ", i + ishift
1036 4 : WRITE (output_unit, '(a)') "HMC| Attempting a minitrajectory move"
1037 4 : WRITE (output_unit, '(a,1x,i0)') "HMC| start mini-trajectory", i + ishift
1038 4 : WRITE (output_unit, '(a,1x,i0,1x)', advance="no") "#HMC|0 ", i + ishift
1039 8 : DO j = 1, force_env%meta_env%n_colvar
1040 4 : WRITE (output_unit, '(f16.8,1x,f16.8,1x,f16.8)', advance="no") force_env%meta_env%metavar(j)%ss0, &
1041 4 : force_env%meta_env%metavar(j)%ss, &
1042 12 : force_env%meta_env%metavar(j)%ff_s !,fz(j)/real(i+ishift,dp)
1043 : END DO
1044 4 : WRITE (output_unit, *)
1045 : END IF
1046 :
1047 : CALL mc_hmc_move(mc_par, force_env, globenv, moves, gmoves, old_epx, old_epz, energy_check, &
1048 8 : r, output_unit, rng_stream_mc, zbuff)
1049 : ! check averages...
1050 : ! force average for z needed too...
1051 : averages%ave_energy = averages%ave_energy*REAL(i - 1, dp)/REAL(i, dp) + &
1052 8 : old_epx/REAL(i, dp)
1053 16 : DO j = 1, force_env%meta_env%n_colvar
1054 16 : fz(j) = fz(j) + force_env%meta_env%metavar(j)%ff_s
1055 : END DO
1056 8 : IF (output_unit > 0) THEN
1057 4 : WRITE (output_unit, '(a,1x,i0)') "HMC|end mini-trajectory", i + ishift
1058 : !!!!!!!! this prints z and theta(x) --ss0,ss-- needed to determine an acceptable k then
1059 : ! the instanteneous force and some instanteneous average for force
1060 4 : WRITE (output_unit, '(a,1x,i0,1x)', advance="no") "#HMC|1 ", i + ishift
1061 8 : DO j = 1, force_env%meta_env%n_colvar
1062 4 : WRITE (output_unit, '(f16.8,1x,f16.8,1x,f16.8,1x,f16.8)', advance="no") force_env%meta_env%metavar(j)%ss0, &
1063 4 : force_env%meta_env%metavar(j)%ss, &
1064 12 : force_env%meta_env%metavar(j)%ff_s, fz(j)/REAL(i + ishift, dp)
1065 : END DO
1066 4 : WRITE (output_unit, *)
1067 : END IF
1068 8 : nsamples = nsamples + 1
1069 12 : IF (MOD(i, iprint) == 0 .AND. (output_unit > 0)) THEN
1070 4 : WRITE (output_unit, '(a,f16.8)') "HMC| Running average for potential energy ", averages%ave_energy
1071 4 : WRITE (output_unit, '(a,1x,i0)') "HMC|======== End Monte Carlo cycle ", i + ishift
1072 : END IF
1073 : ! IF (lrestart) THEN
1074 : ! k=nstep/5
1075 : ! IF(MOD(i,k) == 0) THEN
1076 : ! force_env%qs_env%sim_time=t1
1077 : ! force_env%qs_env%sim_step=it1
1078 : ! DO j=1,force_env%meta_env%n_colvar
1079 : ! force_env%meta_env%metavar(j)%ff_s=fz(j)/real(i+ishift,dp)
1080 : ! ENDDO
1081 : ! ! CALL cp_iterate(logger%iter_info,last=.FALSE.,iter_nr=-1)
1082 : ! CALL section_vals_val_set(mcsec,"RANDOMTOSKIP",i_val=i+ishift)
1083 : ! CALL write_restart(md_env=mdenv,root_section=force_env%root_section)
1084 : ! ! CALL cp_iterate(logger%iter_info,last=.FALSE.,iter_nr=iter)
1085 : ! ENDIF
1086 : ! ENDIF
1087 : END DO
1088 4 : force_env%qs_env%sim_time = t1
1089 4 : force_env%qs_env%sim_step = it1
1090 4 : IF (output_unit > 0) THEN
1091 2 : WRITE (output_unit, '(a,i0,a,i0,a,f16.8)') "HMC| local acceptance ratio: ", moves%hmc%successes, "/", &
1092 4 : moves%hmc%attempts, "=", REAL(moves%hmc%successes, dp)/REAL(moves%hmc%attempts, dp)
1093 2 : WRITE (output_unit, '(a,i0,a,i0,a,f16.8)') "HMC| global acceptance ratio: ", gmoves%hmc%successes, "/", &
1094 4 : gmoves%hmc%attempts, "=", REAL(gmoves%hmc%successes, dp)/REAL(gmoves%hmc%attempts, dp)
1095 : END IF
1096 : !average force
1097 8 : DO j = 1, force_env%meta_env%n_colvar
1098 8 : force_env%meta_env%metavar(j)%ff_s = fz(j)/nsamples
1099 : END DO
1100 4 : END SUBROUTINE HMCsampler
1101 :
1102 : ! **************************************************************************************************
1103 : !> \brief performs a hybrid Monte Carlo move
1104 : !> \param mc_par ...
1105 : !> \param force_env ...
1106 : !> \param globenv ...
1107 : !> \param moves ...
1108 : !> \param gmoves ...
1109 : !> \param old_epx ...
1110 : !> \param old_epz ...
1111 : !> \param energy_check ...
1112 : !> \param r ...
1113 : !> \param output_unit ...
1114 : !> \param rng_stream ...
1115 : !> \param zbuff ...
1116 : !> \author Alin M Elena
1117 : !> \note It runs a NVE trajectory, followed by localisation and accepts rejects
1118 : !> using the biased Hamiltonian, rather than the traditional guiding Hamiltonian
1119 : ! **************************************************************************************************
1120 32 : SUBROUTINE mc_hmc_move(mc_par, force_env, globenv, moves, gmoves, old_epx, old_epz, &
1121 8 : energy_check, r, output_unit, rng_stream, zbuff)
1122 :
1123 : TYPE(mc_simpar_type), POINTER :: mc_par
1124 : TYPE(force_env_type), POINTER :: force_env
1125 : TYPE(global_environment_type), POINTER :: globenv
1126 : TYPE(mc_moves_type), POINTER :: moves, gmoves
1127 : REAL(KIND=dp), INTENT(INOUT) :: old_epx, old_epz, energy_check
1128 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: r
1129 : INTEGER, INTENT(IN) :: output_unit
1130 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
1131 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: zbuff
1132 :
1133 : CHARACTER(LEN=*), PARAMETER :: routineN = 'mc_hmc_move'
1134 :
1135 : INTEGER :: handle, iatom, j, nAtoms, source
1136 : LOGICAL :: ionode, localise
1137 : REAL(KIND=dp) :: BETA, energy_term, exp_max_val, &
1138 : exp_min_val, new_energy, new_epx, &
1139 : new_epz, rand, value, w
1140 : TYPE(cp_subsys_type), POINTER :: oldsys
1141 : TYPE(mc_ekin_type), POINTER :: hmc_ekin
1142 : TYPE(meta_env_type), POINTER :: meta_env_saved
1143 : TYPE(mp_comm_type) :: group
1144 : TYPE(particle_list_type), POINTER :: particles_set
1145 : TYPE(section_vals_type), POINTER :: dft_section, input
1146 :
1147 : ! begin the timing of the subroutine
1148 :
1149 8 : CALL timeset(routineN, handle)
1150 :
1151 8 : CALL get_qs_env(force_env%qs_env, input=input)
1152 8 : dft_section => section_vals_get_subs_vals(input, "DFT")
1153 :
1154 : ! get a bunch of stuff from mc_par
1155 : CALL get_mc_par(mc_par, ionode=ionode, &
1156 : BETA=BETA, exp_max_val=exp_max_val, &
1157 8 : exp_min_val=exp_min_val, source=source, group=group)
1158 :
1159 : ! nullify some pointers
1160 : ! NULLIFY(particles_set,oldsys,hmc_ekin)
1161 8 : NULLIFY (particles_set, oldsys, meta_env_saved, hmc_ekin)
1162 : ! now let's grab the particle positions
1163 8 : CALL force_env_get(force_env, subsys=oldsys)
1164 8 : CALL cp_subsys_get(oldsys, particles=particles_set)
1165 8 : nAtoms = SIZE(particles_set%els)
1166 : ! do some allocation
1167 :
1168 8 : ALLOCATE (hmc_ekin)
1169 :
1170 : ! record the attempt
1171 8 : moves%hmc%attempts = moves%hmc%attempts + 1
1172 8 : gmoves%hmc%attempts = gmoves%hmc%attempts + 1
1173 :
1174 : ! save the old coordinates just in case we need to go back
1175 56 : DO iatom = 1, nAtoms
1176 200 : r(1:3, iatom) = particles_set%els(iatom)%r(1:3)
1177 : END DO
1178 8 : localise = .TRUE.
1179 : ! the same for collective variables data should be made,ss first half and ff_s the last half
1180 16 : DO j = 1, force_env%meta_env%n_colvar
1181 8 : zbuff(j) = force_env%meta_env%metavar(j)%ss
1182 8 : zbuff(j + force_env%meta_env%n_colvar) = force_env%meta_env%metavar(j)%ff_s
1183 8 : IF ((oldsys%colvar_p(force_env%meta_env%metavar(j)%icolvar)%colvar%type_id == HBP_colvar_id) .OR. &
1184 8 : (oldsys%colvar_p(force_env%meta_env%metavar(j)%icolvar)%colvar%type_id == WC_colvar_id)) THEN
1185 8 : localise = .FALSE.
1186 : END IF
1187 : END DO
1188 :
1189 : ! now run the MD simulation
1190 8 : meta_env_saved => force_env%meta_env
1191 8 : NULLIFY (force_env%meta_env)
1192 8 : force_env%qs_env%sim_time = 0.0_dp
1193 8 : force_env%qs_env%sim_step = 0
1194 8 : IF (.NOT. localise) THEN
1195 8 : CALL section_vals_val_set(dft_section, "LOCALIZE%_SECTION_PARAMETERS_", l_val=.FALSE.)
1196 : END IF
1197 8 : CALL qs_mol_dyn(force_env, globenv, hmc_e_initial=hmc_ekin%initial_ekin, hmc_e_final=hmc_ekin%final_ekin)
1198 8 : IF (.NOT. localise) THEN
1199 8 : CALL section_vals_val_set(dft_section, "LOCALIZE%_SECTION_PARAMETERS_", l_val=.TRUE.)
1200 8 : CALL scf_post_calculation_gpw(force_env%qs_env)
1201 : END IF
1202 :
1203 8 : CALL force_env_get(force_env, potential_energy=new_epx)
1204 :
1205 8 : force_env%meta_env => meta_env_saved
1206 8 : CALL tamc_force(force_env, zpot=new_epz)
1207 8 : new_energy = new_epx + new_epz
1208 8 : IF (output_unit > 0) THEN
1209 : WRITE (output_unit, '(a,4(f16.8,1x))') &
1210 4 : "HMC|old Ep, Ekx, Epz, Epx ", old_epx + old_epz, hmc_ekin%initial_ekin, old_epz, old_epx
1211 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
1212 : END IF
1213 8 : energy_term = new_energy - old_epx - old_epz + hmc_ekin%final_ekin - hmc_ekin%initial_ekin
1214 :
1215 8 : value = -BETA*(energy_term)
1216 : ! to prevent overflows
1217 8 : IF (value > exp_max_val) THEN
1218 : w = 10.0_dp
1219 8 : ELSEIF (value < exp_min_val) THEN
1220 : w = 0.0_dp
1221 : ELSE
1222 0 : w = EXP(value)
1223 : END IF
1224 :
1225 8 : rand = rng_stream%next()
1226 8 : IF (rand < w) THEN
1227 : ! accept the move
1228 0 : moves%hmc%successes = moves%hmc%successes + 1
1229 0 : gmoves%hmc%successes = gmoves%hmc%successes + 1
1230 : ! update energies
1231 0 : energy_check = energy_check + (new_energy - old_epx - old_epz)
1232 0 : old_epx = new_epx
1233 0 : old_epz = new_epz
1234 : ELSE
1235 : ! reset the cell and particle positions
1236 56 : DO iatom = 1, nAtoms
1237 200 : particles_set%els(iatom)%r(1:3) = r(1:3, iatom)
1238 : END DO
1239 16 : DO j = 1, force_env%meta_env%n_colvar
1240 8 : force_env%meta_env%metavar(j)%ss = zbuff(j)
1241 16 : force_env%meta_env%metavar(j)%ff_s = zbuff(j + force_env%meta_env%n_colvar)
1242 : END DO
1243 :
1244 : END IF
1245 :
1246 8 : DEALLOCATE (hmc_ekin)
1247 :
1248 : ! end the timing
1249 8 : CALL timestop(handle)
1250 :
1251 8 : END SUBROUTINE mc_hmc_move
1252 :
1253 : ! **************************************************************************************************
1254 : !> \brief ...
1255 : !> \param force_env ...
1256 : ! **************************************************************************************************
1257 4 : SUBROUTINE metadyn_write_colvar_header(force_env)
1258 : TYPE(force_env_type), POINTER :: force_env
1259 :
1260 : CHARACTER(len=*), PARAMETER :: routineN = 'metadyn_write_colvar_header'
1261 :
1262 : CHARACTER(len=100) :: aux, fmt
1263 : CHARACTER(len=255) :: label1, label2, label3, label4, label5, &
1264 : label6
1265 : INTEGER :: handle, i, iw, m
1266 : TYPE(cp_logger_type), POINTER :: logger
1267 : TYPE(meta_env_type), POINTER :: meta_env
1268 :
1269 2 : NULLIFY (logger, meta_env)
1270 2 : meta_env => force_env%meta_env
1271 2 : IF (.NOT. ASSOCIATED(meta_env)) RETURN
1272 :
1273 2 : CALL timeset(routineN, handle)
1274 2 : logger => cp_get_default_logger()
1275 :
1276 : iw = cp_print_key_unit_nr(logger, meta_env%metadyn_section, &
1277 2 : "PRINT%COLVAR", extension=".metadynLog")
1278 2 : IF (iw > 0) THEN
1279 1 : label1 = ""
1280 1 : label2 = ""
1281 1 : label3 = ""
1282 1 : label4 = ""
1283 1 : label5 = ""
1284 1 : label6 = ""
1285 2 : DO i = 1, meta_env%n_colvar
1286 1 : WRITE (aux, '(a,i0)') "z_", i
1287 1 : label1 = TRIM(label1)//TRIM(aux)
1288 1 : m = 15*i - LEN_TRIM(label1) - 1
1289 12 : label1 = TRIM(label1)//REPEAT(" ", m)//"|"
1290 1 : WRITE (aux, '(a,i0)') "Theta_", i
1291 1 : label2 = TRIM(label2)//TRIM(aux)
1292 1 : m = 15*i - LEN_TRIM(label2) - 1
1293 8 : label2 = TRIM(label2)//REPEAT(" ", m)//"|"
1294 1 : WRITE (aux, '(a,i0)') "F_z", i
1295 1 : label3 = TRIM(label3)//TRIM(aux)
1296 1 : m = 15*i - LEN_TRIM(label3) - 1
1297 11 : label3 = TRIM(label3)//REPEAT(" ", m)//"|"
1298 1 : WRITE (aux, '(a,i0)') "F_h", i
1299 1 : label4 = TRIM(label4)//TRIM(aux)
1300 1 : m = 15*i - LEN_TRIM(label4) - 1
1301 11 : label4 = TRIM(label4)//REPEAT(" ", m)//"|"
1302 1 : WRITE (aux, '(a,i0)') "F_w", i
1303 1 : label5 = TRIM(label5)//TRIM(aux)
1304 1 : m = 15*i - LEN_TRIM(label5) - 1
1305 11 : label5 = TRIM(label5)//REPEAT(" ", m)//"|"
1306 1 : WRITE (aux, '(a,i0)') "v_z", i
1307 1 : label6 = TRIM(label6)//TRIM(aux)
1308 1 : m = 15*i - LEN_TRIM(label6) - 1
1309 12 : label6 = TRIM(label6)//REPEAT(" ", m)//"|"
1310 : END DO
1311 1 : WRITE (fmt, '("(a17,6a",i0 ,",4a15)")') meta_env%n_colvar*15
1312 1 : WRITE (iw, TRIM(fmt)) "#Time[fs] |", &
1313 1 : TRIM(label1), &
1314 1 : TRIM(label2), &
1315 1 : TRIM(label3), &
1316 1 : TRIM(label4), &
1317 1 : TRIM(label5), &
1318 1 : TRIM(label6), &
1319 1 : "Epot_z |", &
1320 1 : "Ene hills |", &
1321 1 : "Epot walls |", &
1322 2 : "Temperature |"
1323 :
1324 : END IF
1325 : CALL cp_print_key_finished_output(iw, logger, meta_env%metadyn_section, &
1326 2 : "PRINT%COLVAR")
1327 :
1328 2 : CALL timestop(handle)
1329 :
1330 : END SUBROUTINE metadyn_write_colvar_header
1331 :
1332 : ! **************************************************************************************************
1333 : !> \brief ...
1334 : !> \param force_env ...
1335 : ! **************************************************************************************************
1336 8 : SUBROUTINE metadyn_write_colvar(force_env)
1337 : TYPE(force_env_type), POINTER :: force_env
1338 :
1339 : CHARACTER(len=*), PARAMETER :: routineN = 'metadyn_write_colvar'
1340 :
1341 : INTEGER :: handle, i, i_c, iw
1342 : REAL(KIND=dp) :: temp
1343 : TYPE(cp_logger_type), POINTER :: logger
1344 : TYPE(meta_env_type), POINTER :: meta_env
1345 : TYPE(metavar_type), POINTER :: cv
1346 :
1347 4 : NULLIFY (logger, meta_env, cv)
1348 4 : meta_env => force_env%meta_env
1349 4 : IF (.NOT. ASSOCIATED(meta_env)) RETURN
1350 :
1351 4 : CALL timeset(routineN, handle)
1352 4 : logger => cp_get_default_logger()
1353 :
1354 4 : IF (meta_env%langevin) THEN
1355 4 : meta_env%ekin_s = 0.0_dp
1356 : ! meta_env%epot_s = 0.0_dp
1357 8 : DO i_c = 1, meta_env%n_colvar
1358 4 : cv => meta_env%metavar(i_c)
1359 8 : meta_env%ekin_s = meta_env%ekin_s + 0.5_dp*cv%mass*cv%vvp**2
1360 : END DO
1361 : END IF
1362 :
1363 : ! write COLVAR file
1364 : iw = cp_print_key_unit_nr(logger, meta_env%metadyn_section, &
1365 4 : "PRINT%COLVAR", extension=".metadynLog")
1366 4 : IF (iw > 0) THEN
1367 2 : IF (meta_env%extended_lagrange) THEN
1368 2 : WRITE (iw, '(f16.8,70f15.8)') meta_env%time*femtoseconds, &
1369 6 : (meta_env%metavar(i)%ss0, i=1, meta_env%n_colvar), &
1370 6 : (meta_env%metavar(i)%ss, i=1, meta_env%n_colvar), &
1371 6 : (meta_env%metavar(i)%ff_s, i=1, meta_env%n_colvar), &
1372 6 : (meta_env%metavar(i)%ff_hills, i=1, meta_env%n_colvar), &
1373 6 : (meta_env%metavar(i)%ff_walls, i=1, meta_env%n_colvar), &
1374 6 : (meta_env%metavar(i)%vvp, i=1, meta_env%n_colvar), &
1375 2 : meta_env%epot_s, &
1376 2 : meta_env%hills_env%energy, &
1377 2 : meta_env%epot_walls, &
1378 28 : (meta_env%ekin_s)*2.0_dp/(REAL(meta_env%n_colvar, KIND=dp))*kelvin
1379 : ELSE
1380 0 : WRITE (iw, '(f16.8,40f13.5)') meta_env%time*femtoseconds, &
1381 0 : (meta_env%metavar(i)%ss0, i=1, meta_env%n_colvar), &
1382 0 : (meta_env%metavar(i)%ff_hills, i=1, meta_env%n_colvar), &
1383 0 : (meta_env%metavar(i)%ff_walls, i=1, meta_env%n_colvar), &
1384 0 : meta_env%hills_env%energy, &
1385 0 : meta_env%epot_walls
1386 : END IF
1387 : END IF
1388 : CALL cp_print_key_finished_output(iw, logger, meta_env%metadyn_section, &
1389 4 : "PRINT%COLVAR")
1390 : ! Temperature for COLVAR
1391 4 : IF (meta_env%extended_lagrange) THEN
1392 4 : temp = meta_env%ekin_s*2.0_dp/(REAL(meta_env%n_colvar, KIND=dp))*kelvin
1393 : meta_env%avg_temp = (meta_env%avg_temp*REAL(meta_env%n_steps, KIND=dp) + &
1394 4 : temp)/REAL(meta_env%n_steps + 1, KIND=dp)
1395 : iw = cp_print_key_unit_nr(logger, meta_env%metadyn_section, &
1396 4 : "PRINT%TEMPERATURE_COLVAR", extension=".metadynLog")
1397 4 : IF (iw > 0) THEN
1398 2 : WRITE (iw, '(T2,79("-"))')
1399 2 : WRITE (iw, '( A,T51,f10.2,T71,f10.2)') ' COLVARS INSTANTANEOUS/AVERAGE TEMPERATURE ', &
1400 4 : temp, meta_env%avg_temp
1401 2 : WRITE (iw, '(T2,79("-"))')
1402 : END IF
1403 : CALL cp_print_key_finished_output(iw, logger, meta_env%metadyn_section, &
1404 4 : "PRINT%TEMPERATURE_COLVAR")
1405 : END IF
1406 4 : CALL timestop(handle)
1407 :
1408 : END SUBROUTINE metadyn_write_colvar
1409 :
1410 : ! **************************************************************************************************
1411 : !> \brief ...
1412 : !> \param force_env ...
1413 : ! **************************************************************************************************
1414 2 : SUBROUTINE setup_velocities_z(force_env)
1415 : TYPE(force_env_type), POINTER :: force_env
1416 :
1417 : INTEGER :: i_c
1418 : REAL(kind=dp) :: ekin_w, fac_t
1419 : TYPE(meta_env_type), POINTER :: meta_env
1420 : TYPE(metavar_type), POINTER :: cv
1421 :
1422 2 : NULLIFY (meta_env)
1423 2 : meta_env => force_env%meta_env
1424 2 : meta_env%ekin_s = 0.0_dp
1425 4 : DO i_c = 1, meta_env%n_colvar
1426 2 : cv => meta_env%metavar(i_c)
1427 2 : cv%vvp = force_env%globenv%gaussian_rng_stream%next()
1428 4 : meta_env%ekin_s = meta_env%ekin_s + 0.5_dp*cv%mass*cv%vvp**2
1429 : END DO
1430 2 : ekin_w = 0.5_dp*meta_env%temp_wanted*REAL(meta_env%n_colvar, KIND=dp)
1431 2 : fac_t = SQRT(ekin_w/MAX(meta_env%ekin_s, 1.0E-8_dp))
1432 4 : DO i_c = 1, meta_env%n_colvar
1433 2 : cv => meta_env%metavar(i_c)
1434 4 : cv%vvp = cv%vvp*fac_t
1435 : END DO
1436 2 : END SUBROUTINE setup_velocities_z
1437 : END MODULE tamc_run
|
