(git:6344619)
Loading...
Searching...
No Matches
metadynamics.F
Go to the documentation of this file.
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 Performs the metadynamics calculation
10!> \par History
11!> 01.2005 created [fawzi and ale]
12!> 11.2007 Teodoro Laino [tlaino] - University of Zurich
13! **************************************************************************************************
14MODULE metadynamics
15 USE cp_log_handling, ONLY: cp_logger_type, cp_get_default_logger
16 USE bibliography, ONLY: vandencic2006
17#if defined (__PLUMED2)
18 USE iso_c_binding, ONLY: c_int, c_double, c_char
19 USE cell_types, ONLY: cell_type, &
20 pbc_cp2k_plumed_getset_cell
21#else
22 USE cell_types, ONLY: cell_type
23#endif
24 USE colvar_methods, ONLY: colvar_eval_glob_f
25 USE colvar_types, ONLY: colvar_p_type, &
26 torsion_colvar_id
27 USE constraint_fxd, ONLY: fix_atom_control
28 USE cp_output_handling, ONLY: cp_add_iter_level, &
29 cp_iterate, &
30 cp_print_key_finished_output, &
31 cp_print_key_unit_nr, &
32 cp_rm_iter_level
33 USE cp_subsys_types, ONLY: cp_subsys_get, &
34 cp_subsys_type
35 USE force_env_types, ONLY: force_env_get, &
36 force_env_type
37 USE input_constants, ONLY: do_wall_m, &
38 do_wall_p, &
39 do_wall_reflective
40 USE input_section_types, ONLY: section_vals_get, &
41 section_vals_get_subs_vals, &
42 section_vals_type, &
43 section_vals_val_get
44 USE kinds, ONLY: dp
45 USE message_passing, ONLY: cp2k_is_parallel
46 USE metadynamics_types, ONLY: hills_env_type, &
47 meta_env_type, &
48 metavar_type, &
49 multiple_walkers_type
50 USE metadynamics_utils, ONLY: add_hill_single, &
51 get_meta_iter_level, &
52 meta_walls, &
53 restart_hills, &
54 synchronize_multiple_walkers
55 USE particle_list_types, ONLY: particle_list_type
56#if defined (__PLUMED2)
57 USE physcon, ONLY: angstrom, &
58 boltzmann, &
59 femtoseconds, &
60 joule, &
61 kelvin, kjmol, &
62 picoseconds
63#else
64 USE physcon, ONLY: boltzmann, &
65 femtoseconds, &
66 joule, &
67 kelvin
68#endif
69 USE reference_manager, ONLY: cite_reference
70 USE simpar_types, ONLY: simpar_type
71#if defined (__PLUMED2)
72 USE virial_types, ONLY: virial_type
73#endif
74#include "./base/base_uses.f90"
75
76 IMPLICIT NONE
77 PRIVATE
78
79 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'metadynamics'
80
81 PUBLIC :: metadyn_forces, metadyn_integrator
82 PUBLIC :: metadyn_velocities_colvar, metadyn_write_colvar
83 PUBLIC :: metadyn_initialise_plumed, metadyn_finalise_plumed
84
85#if defined (__PLUMED2)
86 INTERFACE
87 FUNCTION plumed_installed() RESULT(res) BIND(C, name="plumed_installed")
88 IMPORT :: c_int
89 INTEGER(KIND=C_INT) :: res
90 END FUNCTION plumed_installed
91
92 SUBROUTINE plumed_gcreate() BIND(C, name="plumed_gcreate")
93 END SUBROUTINE plumed_gcreate
94
95 SUBROUTINE plumed_gfinalize() BIND(C, name="plumed_gfinalize")
96 END SUBROUTINE plumed_gfinalize
97
98 SUBROUTINE plumed_gcmd_int(key, value) BIND(C, name="plumed_gcmd")
99 IMPORT :: c_char, c_int
100 CHARACTER(KIND=C_CHAR), DIMENSION(*) :: key
101 INTEGER(KIND=C_INT) :: value
102 END SUBROUTINE plumed_gcmd_int
103
104 SUBROUTINE plumed_gcmd_double(key, value) BIND(C, name="plumed_gcmd")
105 IMPORT :: c_char, c_double
106 CHARACTER(KIND=C_CHAR), DIMENSION(*) :: key
107 ! PLUMED retains pointers passed by setEnergy (and the array setters).
108 REAL(KIND=c_double), TARGET :: value
109 END SUBROUTINE plumed_gcmd_double
110
111 SUBROUTINE plumed_gcmd_doubles(key, value) BIND(C, name="plumed_gcmd")
112 IMPORT :: c_char, c_double
113 CHARACTER(KIND=C_CHAR), DIMENSION(*) :: key
114 REAL(KIND=c_double), TARGET, DIMENSION(*) :: value
115 END SUBROUTINE plumed_gcmd_doubles
116
117 SUBROUTINE plumed_gcmd_char(key, value) BIND(C, name="plumed_gcmd")
118 IMPORT :: c_char
119 CHARACTER(KIND=C_CHAR), DIMENSION(*) :: key, value
120 END SUBROUTINE plumed_gcmd_char
121 END INTERFACE
122#endif
123
124CONTAINS
125
126! **************************************************************************************************
127!> \brief ...
128!> \param force_env ...
129!> \param simpar ...
130!> \param itimes ...
131! **************************************************************************************************
132 SUBROUTINE metadyn_initialise_plumed(force_env, simpar, itimes)
133
134 TYPE(force_env_type), POINTER :: force_env
135 TYPE(simpar_type), POINTER :: simpar
136 INTEGER, INTENT(IN) :: itimes
137
138 CHARACTER(LEN=*), PARAMETER :: routinen = 'metadyn_initialise_plumed'
139
140 INTEGER :: handle
141
142#if defined (__PLUMED2)
143 INTEGER :: natom_plumed
144 REAL(kind=dp) :: timestep_plumed
145 TYPE(cell_type), POINTER :: cell
146 TYPE(cp_subsys_type), POINTER :: subsys
147#endif
148
149#if defined (__PLUMED2)
150 INTEGER :: plumedavailable
151#endif
152
153 CALL timeset(routinen, handle)
154 cpassert(ASSOCIATED(force_env))
155 cpassert(ASSOCIATED(simpar))
156
157#if defined (__PLUMED2)
158 NULLIFY (cell, subsys)
159 CALL force_env_get(force_env, subsys=subsys, cell=cell)
160 CALL pbc_cp2k_plumed_getset_cell(cell, set=.true.) !Store the cell pointer for later use.
161 timestep_plumed = simpar%dt
162 natom_plumed = subsys%particles%n_els
163#else
164 mark_used(force_env)
165 mark_used(simpar)
166#endif
167
168#if !defined (__PLUMED2)
169 mark_used(itimes)
170#endif
171#if defined (__PLUMED2)
172 plumedavailable = plumed_installed()
173
174 IF (plumedavailable <= 0) THEN
175 cpabort("NO PLUMED IN YOUR LD_LIBRARY_PATH?")
176 ELSE
177 CALL plumed_gcreate()
178 IF (cp2k_is_parallel) CALL plumed_gcmd_int("setMPIFComm"//char(0), force_env%para_env%get_handle())
179 CALL plumed_gcmd_int("setRealPrecision"//char(0), 8)
180 CALL plumed_gcmd_double("setMDEnergyUnits"//char(0), kjmol) ! Hartree to kjoule/mol
181 CALL plumed_gcmd_double("setMDLengthUnits"//char(0), angstrom*0.1_dp) ! Bohr to nm
182 CALL plumed_gcmd_double("setMDTimeUnits"//char(0), picoseconds) ! atomic units to ps
183 CALL plumed_gcmd_char("setPlumedDat"//char(0), force_env%meta_env%plumed_input_file//char(0))
184 CALL plumed_gcmd_char("setLogFile"//char(0), "PLUMED.OUT"//char(0))
185 CALL plumed_gcmd_int("setNatoms"//char(0), natom_plumed)
186 CALL plumed_gcmd_char("setMDEngine"//char(0), "cp2k"//char(0))
187 CALL plumed_gcmd_double("setTimestep"//char(0), timestep_plumed)
188 ! CP2K stores kBT in atomic energy units, not kelvin.
189 IF (simpar%temp_ext > 0.0_dp) CALL plumed_gcmd_double("setKbT"//char(0), simpar%temp_ext)
190 IF (itimes > 0) CALL plumed_gcmd_int("setRestart"//char(0), 1)
191 CALL plumed_gcmd_int("init"//char(0), 0)
192 END IF
193#endif
194
195#if !defined (__PLUMED2)
196 CALL cp_abort(__location__, &
197 "Requested to use plumed for metadynamics, but cp2k was"// &
198 " not compiled with plumed support.")
199#endif
200
201 CALL timestop(handle)
202
203 END SUBROUTINE metadyn_initialise_plumed
204
205! **************************************************************************************************
206!> \brief makes sure PLUMED is shut down cleanly
207! **************************************************************************************************
208 SUBROUTINE metadyn_finalise_plumed()
209
210 CHARACTER(LEN=*), PARAMETER :: routinen = 'metadyn_finalise_plumed'
211
212 INTEGER :: handle
213
214 CALL timeset(routinen, handle)
215
216#if defined (__PLUMED2)
217 CALL plumed_gfinalize()
218#endif
219
220 CALL timestop(handle)
221
222 END SUBROUTINE metadyn_finalise_plumed
223
224! **************************************************************************************************
225!> \brief General driver for applying metadynamics
226!> \param force_env ...
227!> \param itimes ...
228!> \param vel ...
229!> \param rand ...
230!> \date 01.2009
231!> \par History
232!> 01.2009 created
233!> \author Teodoro Laino
234! **************************************************************************************************
235 SUBROUTINE metadyn_integrator(force_env, itimes, vel, rand)
236 TYPE(force_env_type), POINTER :: force_env
237 INTEGER, INTENT(IN) :: itimes
238 REAL(kind=dp), DIMENSION(:, :), &
239 INTENT(INOUT), OPTIONAL :: vel
240 REAL(kind=dp), DIMENSION(:), OPTIONAL, &
241 POINTER :: rand
242
243 CHARACTER(len=*), PARAMETER :: routinen = 'metadyn_integrator'
244
245 INTEGER :: handle, plumed_itimes
246#if defined (__PLUMED2)
247 INTEGER :: i_part, natom_plumed
248 TYPE(virial_type), POINTER :: virial
249 TYPE(cell_type), POINTER :: cell
250 TYPE(cp_subsys_type), POINTER :: subsys
251#endif
252#if defined (__PLUMED2)
253 TYPE(meta_env_type), POINTER :: meta_env
254 REAL(kind=dp), ALLOCATABLE, TARGET, DIMENSION(:) :: mass_plumed
255 REAL(kind=dp), TARGET, DIMENSION(3, 3) :: plu_vir, celltbt
256 REAL(kind=dp), TARGET :: stpcfg
257 REAL(kind=dp), ALLOCATABLE, TARGET, &
258 DIMENSION(:) :: pos_plumed_x, pos_plumed_y, pos_plumed_z
259 REAL(kind=dp), ALLOCATABLE, TARGET, &
260 DIMENSION(:) :: force_plumed_x, force_plumed_y, force_plumed_z
261#endif
262
263 CALL timeset(routinen, handle)
264
265 ! Apply Metadynamics
266 IF (ASSOCIATED(force_env%meta_env)) THEN
267 IF (force_env%meta_env%use_plumed .EQV. .true.) THEN
268 plumed_itimes = itimes
269#if defined (__PLUMED2)
270 CALL force_env_get(force_env, meta_env=meta_env, subsys=subsys, cell=cell)
271 CALL cp_subsys_get(subsys, virial=virial)
272 natom_plumed = subsys%particles%n_els
273 ALLOCATE (pos_plumed_x(natom_plumed))
274 ALLOCATE (pos_plumed_y(natom_plumed))
275 ALLOCATE (pos_plumed_z(natom_plumed))
276
277 ALLOCATE (force_plumed_x(natom_plumed))
278 ALLOCATE (force_plumed_y(natom_plumed))
279 ALLOCATE (force_plumed_z(natom_plumed))
280
281 ALLOCATE (mass_plumed(natom_plumed))
282
283 DO i_part = 1, natom_plumed
284 pos_plumed_x(i_part) = subsys%particles%els(i_part)%r(1)
285 pos_plumed_y(i_part) = subsys%particles%els(i_part)%r(2)
286 pos_plumed_z(i_part) = subsys%particles%els(i_part)%r(3)
287 mass_plumed(i_part) = subsys%particles%els(i_part)%atomic_kind%mass
288 ! ENERGY biases rescale physical forces, so these must be inputs too.
289 force_plumed_x(i_part) = subsys%particles%els(i_part)%f(1)
290 force_plumed_y(i_part) = subsys%particles%els(i_part)%f(2)
291 force_plumed_z(i_part) = subsys%particles%els(i_part)%f(3)
292 END DO
293
294 ! stupid cell conversion, to be fixed
295 celltbt(1, 1) = cell%hmat(1, 1)
296 celltbt(1, 2) = cell%hmat(1, 2)
297 celltbt(1, 3) = cell%hmat(1, 3)
298 celltbt(2, 1) = cell%hmat(2, 1)
299 celltbt(2, 2) = cell%hmat(2, 2)
300 celltbt(2, 3) = cell%hmat(2, 3)
301 celltbt(3, 1) = cell%hmat(3, 1)
302 celltbt(3, 2) = cell%hmat(3, 2)
303 celltbt(3, 3) = cell%hmat(3, 3)
304
305 ! PLUMED uses the opposite sign to CP2K's potential virial.
306 plu_vir(:, :) = 0.0d0
307 IF (virial%pv_availability) plu_vir = -virial%pv_virial
308
309 ! Do not feed the previous step's bias back into the ENERGY CV.
310 meta_env%plumed_bias = 0.0_dp
311 CALL force_env_get(force_env, potential_energy=stpcfg)
312
313 CALL plumed_gcmd_int("setStep"//char(0), plumed_itimes)
314 CALL plumed_gcmd_doubles("setPositionsX"//char(0), pos_plumed_x(:))
315 CALL plumed_gcmd_doubles("setPositionsY"//char(0), pos_plumed_y(:))
316 CALL plumed_gcmd_doubles("setPositionsZ"//char(0), pos_plumed_z(:))
317 CALL plumed_gcmd_doubles("setMasses"//char(0), mass_plumed(:))
318 CALL plumed_gcmd_doubles("setBox"//char(0), celltbt)
319 CALL plumed_gcmd_doubles("setVirial"//char(0), plu_vir)
320 CALL plumed_gcmd_double("setEnergy"//char(0), stpcfg)
321 CALL plumed_gcmd_doubles("setForcesX"//char(0), force_plumed_x(:)) ! PluMed Output !
322 CALL plumed_gcmd_doubles("setForcesY"//char(0), force_plumed_y(:)) ! PluMed Output !
323 CALL plumed_gcmd_doubles("setForcesZ"//char(0), force_plumed_z(:)) ! PluMed Output !
324
325 CALL plumed_gcmd_int("prepareCalc"//char(0), 0)
326
327 CALL plumed_gcmd_int("performCalcNoUpdate"//char(0), 0)
328 CALL plumed_gcmd_double("getBias"//char(0), meta_env%plumed_bias)
329 CALL plumed_gcmd_int("update"//char(0), 0)
330 ! PLUMED updates the box force on rank zero. CP2K stores a replicated
331 ! potential virial, including for the barostat on all MPI ranks.
332 CALL force_env%para_env%bcast(plu_vir, source=0)
333 IF (virial%pv_availability) virial%pv_virial = -plu_vir
334
335 DO i_part = 1, natom_plumed
336 subsys%particles%els(i_part)%f(1) = force_plumed_x(i_part)
337 subsys%particles%els(i_part)%f(2) = force_plumed_y(i_part)
338 subsys%particles%els(i_part)%f(3) = force_plumed_z(i_part)
339 END DO
340
341 DEALLOCATE (pos_plumed_x, pos_plumed_y, pos_plumed_z)
342 DEALLOCATE (force_plumed_x, force_plumed_y, force_plumed_z)
343 DEALLOCATE (mass_plumed)
344
345 ! Constraints ONLY of Fixed Atom type
346 CALL fix_atom_control(force_env)
347#endif
348
349#if !defined (__PLUMED2)
350 CALL cp_abort(__location__, &
351 "Requested to use plumed for metadynamics, but cp2k was"// &
352 " not compiled with plumed support.")
353#endif
354
355 ELSE
356 IF (force_env%meta_env%langevin) THEN
357 IF (.NOT. PRESENT(rand)) THEN
358 cpabort("Langevin on COLVAR not implemented for this MD ensemble!")
359 END IF
360 ! *** Velocity Verlet for Langevin S(t)->S(t+1)
361 CALL metadyn_position_colvar(force_env)
362 ! *** Forces from Vs and S(X(t+1))
363 CALL metadyn_forces(force_env)
364 ! *** Velocity Verlet for Langeving *** v(t+1/2)--> v(t)
365 CALL metadyn_velocities_colvar(force_env, rand)
366 ELSE
367 CALL metadyn_forces(force_env, vel)
368 END IF
369 ! *** Write down COVAR informations
370 CALL metadyn_write_colvar(force_env)
371 END IF
372 END IF
373
374 CALL timestop(handle)
375
376 END SUBROUTINE metadyn_integrator
377
378! **************************************************************************************************
379!> \brief add forces to the subsys due to the metadynamics run
380!> possibly modifies the velocites (if reflective walls are applied)
381!> \param force_env ...
382!> \param vel ...
383!> \par History
384!> 04.2004 created
385! **************************************************************************************************
386 SUBROUTINE metadyn_forces(force_env, vel)
387 TYPE(force_env_type), POINTER :: force_env
388 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT), &
389 OPTIONAL :: vel
390
391 CHARACTER(len=*), PARAMETER :: routinen = 'metadyn_forces'
392
393 INTEGER :: handle, i, i_c, icolvar, ii, iwall
394 LOGICAL :: explicit
395 REAL(kind=dp) :: check_val, diff_ss, dt, ekin_w, fac_t, &
396 fft, norm, rval, scal, scalf, &
397 ss0_test, tol_ekin
398 TYPE(colvar_p_type), DIMENSION(:), POINTER :: colvar_p
399 TYPE(cp_logger_type), POINTER :: logger
400 TYPE(cp_subsys_type), POINTER :: subsys
401 TYPE(meta_env_type), POINTER :: meta_env
402 TYPE(metavar_type), POINTER :: cv
403 TYPE(particle_list_type), POINTER :: particles
404 TYPE(section_vals_type), POINTER :: ss0_section, vvp_section
405
406 NULLIFY (logger, meta_env)
407 meta_env => force_env%meta_env
408 IF (.NOT. ASSOCIATED(meta_env)) RETURN
409
410 CALL timeset(routinen, handle)
411 logger => cp_get_default_logger()
412 NULLIFY (colvar_p, subsys, cv, ss0_section, vvp_section)
413 CALL force_env_get(force_env, subsys=subsys)
414
415 dt = meta_env%dt
416 IF (.NOT. meta_env%restart) meta_env%n_steps = meta_env%n_steps + 1
417
418 ! Initialize velocity
419 IF (meta_env%restart .AND. meta_env%extended_lagrange) THEN
420 meta_env%ekin_s = 0.0_dp
421 DO i_c = 1, meta_env%n_colvar
422 cv => meta_env%metavar(i_c)
423 cv%vvp = force_env%globenv%gaussian_rng_stream%next()
424 meta_env%ekin_s = meta_env%ekin_s + 0.5_dp*cv%mass*cv%vvp**2
425 END DO
426 ekin_w = 0.5_dp*meta_env%temp_wanted*real(meta_env%n_colvar, kind=dp)
427 fac_t = sqrt(ekin_w/max(meta_env%ekin_s, 1.0e-8_dp))
428 DO i_c = 1, meta_env%n_colvar
429 cv => meta_env%metavar(i_c)
430 cv%vvp = cv%vvp*fac_t
431 END DO
432 meta_env%ekin_s = 0.0_dp
433 END IF
434
435 ! *** Velocity Verlet for Langevin S(t)->S(t+1)
436 ! compute ss and the derivative of ss with respect to the atomic positions
437 DO i_c = 1, meta_env%n_colvar
438 cv => meta_env%metavar(i_c)
439 icolvar = cv%icolvar
440 CALL colvar_eval_glob_f(icolvar, force_env)
441 cv%ss = subsys%colvar_p(icolvar)%colvar%ss
442
443 ! Setup the periodic flag if the COLVAR is (-pi,pi] periodic
444 cv%periodic = (subsys%colvar_p(icolvar)%colvar%type_id == torsion_colvar_id)
445
446 ! Restart for Extended Lagrangian Metadynamics
447 IF (meta_env%restart) THEN
448 ! Initialize the position of the collective variable in the extended lagrange
449 ss0_section => section_vals_get_subs_vals(meta_env%metadyn_section, "EXT_LAGRANGE_SS0")
450 CALL section_vals_get(ss0_section, explicit=explicit)
451 IF (explicit) THEN
452 CALL section_vals_val_get(ss0_section, "_DEFAULT_KEYWORD_", &
453 i_rep_val=i_c, r_val=rval)
454 cv%ss0 = rval
455 ELSE
456 cv%ss0 = cv%ss
457 END IF
458 vvp_section => section_vals_get_subs_vals(meta_env%metadyn_section, "EXT_LAGRANGE_VVP")
459 CALL section_vals_get(vvp_section, explicit=explicit)
460 IF (explicit) THEN
461 CALL section_vals_val_get(vvp_section, "_DEFAULT_KEYWORD_", &
462 i_rep_val=i_c, r_val=rval)
463 cv%vvp = rval
464 END IF
465 END IF
466 !
467 IF (.NOT. meta_env%extended_lagrange) THEN
468 cv%ss0 = cv%ss
469 cv%vvp = 0.0_dp
470 END IF
471 END DO
472 ! History dependent forces (evaluated at s0)
473 IF (meta_env%do_hills) CALL hills(meta_env)
474
475 ! Apply walls to the colvars
476 CALL meta_walls(meta_env)
477
478 meta_env%restart = .false.
479 IF (.NOT. meta_env%extended_lagrange) THEN
480 meta_env%ekin_s = 0.0_dp
481 meta_env%epot_s = 0.0_dp
482 meta_env%epot_walls = 0.0_dp
483 DO i_c = 1, meta_env%n_colvar
484 cv => meta_env%metavar(i_c)
485 cv%epot_s = 0.0_dp
486 cv%ff_s = 0.0_dp
487 meta_env%epot_walls = meta_env%epot_walls + cv%epot_walls
488 icolvar = cv%icolvar
489 NULLIFY (particles)
490 CALL cp_subsys_get(subsys, colvar_p=colvar_p, &
491 particles=particles)
492 DO ii = 1, colvar_p(icolvar)%colvar%n_atom_s
493 i = colvar_p(icolvar)%colvar%i_atom(ii)
494 fft = cv%ff_hills + cv%ff_walls
495 particles%els(i)%f = particles%els(i)%f + fft*colvar_p(icolvar)%colvar%dsdr(:, ii)
496 END DO
497 END DO
498 ELSE
499 meta_env%ekin_s = 0.0_dp
500 meta_env%epot_s = 0.0_dp
501 meta_env%epot_walls = 0.0_dp
502 DO i_c = 1, meta_env%n_colvar
503 cv => meta_env%metavar(i_c)
504 diff_ss = cv%ss - cv%ss0
505 IF (cv%periodic) THEN
506 ! The difference of a periodic COLVAR is always within [-pi,pi]
507 diff_ss = sign(1.0_dp, asin(sin(diff_ss)))*acos(cos(diff_ss))
508 END IF
509 cv%epot_s = 0.5_dp*cv%lambda*(diff_ss)**2.0_dp
510 cv%ff_s = cv%lambda*(diff_ss)
511 icolvar = cv%icolvar
512 ! forces on the atoms
513 NULLIFY (particles)
514 CALL cp_subsys_get(subsys, colvar_p=colvar_p, &
515 particles=particles)
516 DO ii = 1, colvar_p(icolvar)%colvar%n_atom_s
517 i = colvar_p(icolvar)%colvar%i_atom(ii)
518 particles%els(i)%f = particles%els(i)%f - cv%ff_s*colvar_p(icolvar)%colvar%dsdr(:, ii)
519 END DO
520 ! velocity verlet on the s0 if NOT langevin
521 IF (.NOT. meta_env%langevin) THEN
522 fft = cv%ff_s + cv%ff_hills + cv%ff_walls
523 cv%vvp = cv%vvp + dt*fft/cv%mass
524 meta_env%ekin_s = meta_env%ekin_s + 0.5_dp*cv%mass*cv%vvp**2
525 meta_env%epot_s = meta_env%epot_s + cv%epot_s
526 meta_env%epot_walls = meta_env%epot_walls + cv%epot_walls
527 END IF
528 END DO
529 ! velocity rescaling on the s0
530 IF (meta_env%tempcontrol .AND. (.NOT. meta_env%langevin)) THEN
531 ekin_w = 0.5_dp*meta_env%temp_wanted*real(meta_env%n_colvar, kind=dp)
532 tol_ekin = 0.5_dp*meta_env%toll_temp*real(meta_env%n_colvar, kind=dp)
533 IF (abs(ekin_w - meta_env%ekin_s) > tol_ekin) THEN
534 fac_t = sqrt(ekin_w/max(meta_env%ekin_s, 1.0e-8_dp))
535 DO i_c = 1, meta_env%n_colvar
536 cv => meta_env%metavar(i_c)
537 cv%vvp = cv%vvp*fac_t
538 END DO
539 meta_env%ekin_s = ekin_w
540 END IF
541 END IF
542 ! Reflective Wall only for s0
543 DO i_c = 1, meta_env%n_colvar
544 cv => meta_env%metavar(i_c)
545 IF (cv%do_wall) THEN
546 DO iwall = 1, SIZE(cv%walls)
547 SELECT CASE (cv%walls(iwall)%id_type)
548 CASE (do_wall_reflective)
549 ss0_test = cv%ss0 + dt*cv%vvp
550 IF (cv%periodic) THEN
551 ! A periodic COLVAR is always within [-pi,pi]
552 ss0_test = sign(1.0_dp, asin(sin(ss0_test)))*acos(cos(ss0_test))
553 END IF
554 SELECT CASE (cv%walls(iwall)%id_direction)
555 CASE (do_wall_p)
556 IF ((ss0_test > cv%walls(iwall)%pos) .AND. (cv%vvp > 0)) cv%vvp = -cv%vvp
557 CASE (do_wall_m)
558 IF ((ss0_test < cv%walls(iwall)%pos) .AND. (cv%vvp < 0)) cv%vvp = -cv%vvp
559 END SELECT
560 END SELECT
561 END DO
562 END IF
563 END DO
564 ! Update of ss0 if NOT langevin
565 IF (.NOT. meta_env%langevin) THEN
566 DO i_c = 1, meta_env%n_colvar
567 cv => meta_env%metavar(i_c)
568 cv%ss0 = cv%ss0 + dt*cv%vvp
569 IF (cv%periodic) THEN
570 ! A periodic COLVAR is always within [-pi,pi]
571 cv%ss0 = sign(1.0_dp, asin(sin(cv%ss0)))*acos(cos(cv%ss0))
572 END IF
573 END DO
574 END IF
575 END IF
576 ! Constraints ONLY of Fixed Atom type
577 CALL fix_atom_control(force_env)
578
579 ! Reflective Wall only for ss
580 DO i_c = 1, meta_env%n_colvar
581 cv => meta_env%metavar(i_c)
582 IF (cv%do_wall) THEN
583 DO iwall = 1, SIZE(cv%walls)
584 SELECT CASE (cv%walls(iwall)%id_type)
585 CASE (do_wall_reflective)
586 SELECT CASE (cv%walls(iwall)%id_direction)
587 CASE (do_wall_p)
588 IF (cv%ss < cv%walls(iwall)%pos) cycle
589 check_val = -1.0_dp
590 CASE (do_wall_m)
591 IF (cv%ss > cv%walls(iwall)%pos) cycle
592 check_val = 1.0_dp
593 END SELECT
594 NULLIFY (particles)
595 icolvar = cv%icolvar
596 CALL cp_subsys_get(subsys, colvar_p=colvar_p, particles=particles)
597 scal = 0.0_dp
598 scalf = 0.0_dp
599 norm = 0.0_dp
600 DO ii = 1, colvar_p(icolvar)%colvar%n_atom_s
601 i = colvar_p(icolvar)%colvar%i_atom(ii)
602 IF (PRESENT(vel)) THEN
603 scal = scal + vel(1, i)*colvar_p(icolvar)%colvar%dsdr(1, ii)
604 scal = scal + vel(2, i)*colvar_p(icolvar)%colvar%dsdr(2, ii)
605 scal = scal + vel(3, i)*colvar_p(icolvar)%colvar%dsdr(3, ii)
606 ELSE
607 scal = scal + particles%els(i)%v(1)*colvar_p(icolvar)%colvar%dsdr(1, ii)
608 scal = scal + particles%els(i)%v(2)*colvar_p(icolvar)%colvar%dsdr(2, ii)
609 scal = scal + particles%els(i)%v(3)*colvar_p(icolvar)%colvar%dsdr(3, ii)
610 END IF
611 scalf = scalf + particles%els(i)%f(1)*colvar_p(icolvar)%colvar%dsdr(1, ii)
612 scalf = scalf + particles%els(i)%f(2)*colvar_p(icolvar)%colvar%dsdr(2, ii)
613 scalf = scalf + particles%els(i)%f(3)*colvar_p(icolvar)%colvar%dsdr(3, ii)
614 norm = norm + colvar_p(icolvar)%colvar%dsdr(1, ii)**2
615 norm = norm + colvar_p(icolvar)%colvar%dsdr(2, ii)**2
616 norm = norm + colvar_p(icolvar)%colvar%dsdr(3, ii)**2
617 END DO
618 IF (norm /= 0.0_dp) scal = scal/norm
619 IF (norm /= 0.0_dp) scalf = scalf/norm
620
621 IF (scal*check_val > 0.0_dp) cycle
622 DO ii = 1, colvar_p(icolvar)%colvar%n_atom_s
623 i = colvar_p(icolvar)%colvar%i_atom(ii)
624 IF (PRESENT(vel)) THEN
625 vel(:, i) = vel(:, i) - 2.0_dp*colvar_p(icolvar)%colvar%dsdr(:, ii)*scal
626 ELSE
627 particles%els(i)%v(:) = particles%els(i)%v(:) - 2.0_dp*colvar_p(icolvar)%colvar%dsdr(:, ii)*scal
628 END IF
629 ! Nullify forces along the colvar (this avoids the weird behaviors of the reflective wall)
630 particles%els(i)%f(:) = particles%els(i)%f(:) - colvar_p(icolvar)%colvar%dsdr(:, ii)*scalf
631 END DO
632 END SELECT
633 END DO
634 END IF
635 END DO
636
637 CALL timestop(handle)
638 END SUBROUTINE metadyn_forces
639
640! **************************************************************************************************
641!> \brief Evolves velocities COLVAR according to
642!> Vanden-Eijnden Ciccotti C.Phys.Letter 429 (2006) 310-316
643!> \param force_env ...
644!> \param rand ...
645!> \date 01.2009
646!> \author Fabio Sterpone and Teodoro Laino
647! **************************************************************************************************
648 SUBROUTINE metadyn_velocities_colvar(force_env, rand)
649 TYPE(force_env_type), POINTER :: force_env
650 REAL(kind=dp), DIMENSION(:), INTENT(INOUT), &
651 OPTIONAL :: rand
652
653 CHARACTER(len=*), PARAMETER :: routinen = 'metadyn_velocities_colvar'
654
655 INTEGER :: handle, i_c
656 REAL(kind=dp) :: diff_ss, dt, fft, sigma
657 TYPE(cp_logger_type), POINTER :: logger
658 TYPE(meta_env_type), POINTER :: meta_env
659 TYPE(metavar_type), POINTER :: cv
660
661 NULLIFY (logger, meta_env, cv)
662 meta_env => force_env%meta_env
663 IF (.NOT. ASSOCIATED(meta_env)) RETURN
664
665 CALL timeset(routinen, handle)
666 logger => cp_get_default_logger()
667 ! Add citation
668 IF (meta_env%langevin) CALL cite_reference(vandencic2006)
669
670 dt = meta_env%dt
671 ! History dependent forces (evaluated at s0)
672 IF (meta_env%do_hills) CALL hills(meta_env)
673
674 ! Evolve Velocities
675 meta_env%ekin_s = 0.0_dp
676 meta_env%epot_walls = 0.0_dp
677 DO i_c = 1, meta_env%n_colvar
678 cv => meta_env%metavar(i_c)
679 diff_ss = cv%ss - cv%ss0
680 IF (cv%periodic) THEN
681 ! The difference of a periodic COLVAR is always within [-pi,pi]
682 diff_ss = sign(1.0_dp, asin(sin(diff_ss)))*acos(cos(diff_ss))
683 END IF
684 cv%epot_s = 0.5_dp*cv%lambda*(diff_ss)**2.0_dp
685 cv%ff_s = cv%lambda*(diff_ss)
686
687 fft = cv%ff_s + cv%ff_hills
688 sigma = sqrt((meta_env%temp_wanted*kelvin)*2.0_dp*(boltzmann/joule)*cv%gamma/cv%mass)
689 cv%vvp = cv%vvp + 0.5_dp*dt*fft/cv%mass - 0.5_dp*dt*cv%gamma*cv%vvp + &
690 0.5_dp*sqrt(dt)*sigma*rand(i_c)
691 meta_env%ekin_s = meta_env%ekin_s + 0.5_dp*cv%mass*cv%vvp**2
692 meta_env%epot_walls = meta_env%epot_walls + cv%epot_walls
693 END DO
694 CALL timestop(handle)
695
696 END SUBROUTINE metadyn_velocities_colvar
697
698! **************************************************************************************************
699!> \brief Evolves COLVAR position
700!> Vanden-Eijnden Ciccotti C.Phys.Letter 429 (2006) 310-316
701!> \param force_env ...
702!> \date 01.2009
703!> \author Fabio Sterpone and Teodoro Laino
704! **************************************************************************************************
705 SUBROUTINE metadyn_position_colvar(force_env)
706 TYPE(force_env_type), POINTER :: force_env
707
708 CHARACTER(len=*), PARAMETER :: routinen = 'metadyn_position_colvar'
709
710 INTEGER :: handle, i_c
711 REAL(kind=dp) :: dt
712 TYPE(cp_logger_type), POINTER :: logger
713 TYPE(meta_env_type), POINTER :: meta_env
714 TYPE(metavar_type), POINTER :: cv
715
716 NULLIFY (logger, meta_env, cv)
717 meta_env => force_env%meta_env
718 IF (.NOT. ASSOCIATED(meta_env)) RETURN
719
720 CALL timeset(routinen, handle)
721 logger => cp_get_default_logger()
722
723 ! Add citation
724 IF (meta_env%langevin) CALL cite_reference(vandencic2006)
725 dt = meta_env%dt
726
727 ! Update of ss0
728 DO i_c = 1, meta_env%n_colvar
729 cv => meta_env%metavar(i_c)
730 cv%ss0 = cv%ss0 + dt*cv%vvp
731 IF (cv%periodic) THEN
732 ! A periodic COLVAR is always within [-pi,pi]
733 cv%ss0 = sign(1.0_dp, asin(sin(cv%ss0)))*acos(cos(cv%ss0))
734 END IF
735 END DO
736 CALL timestop(handle)
737
738 END SUBROUTINE metadyn_position_colvar
739
740! **************************************************************************************************
741!> \brief Write down COLVAR information evolved according to
742!> Vanden-Eijnden Ciccotti C.Phys.Letter 429 (2006) 310-316
743!> \param force_env ...
744!> \date 01.2009
745!> \author Fabio Sterpone and Teodoro Laino
746! **************************************************************************************************
747 SUBROUTINE metadyn_write_colvar(force_env)
748 TYPE(force_env_type), POINTER :: force_env
749
750 CHARACTER(len=*), PARAMETER :: routinen = 'metadyn_write_colvar'
751
752 INTEGER :: handle, i, i_c, iw
753 REAL(kind=dp) :: diff_ss, temp
754 TYPE(cp_logger_type), POINTER :: logger
755 TYPE(meta_env_type), POINTER :: meta_env
756 TYPE(metavar_type), POINTER :: cv
757
758 NULLIFY (logger, meta_env, cv)
759 meta_env => force_env%meta_env
760 IF (.NOT. ASSOCIATED(meta_env)) RETURN
761
762 CALL timeset(routinen, handle)
763 logger => cp_get_default_logger()
764
765 ! If Langevin we need to recompute few quantities
766 ! This does not apply to the standard lagrangian scheme since it is
767 ! implemented with a plain Newton integration scheme.. while Langevin
768 ! follows the correct Verlet integration.. This will have to be made
769 ! uniform in the future (Teodoro Laino - 01.2009)
770 IF (meta_env%langevin) THEN
771 meta_env%ekin_s = 0.0_dp
772 meta_env%epot_s = 0.0_dp
773 DO i_c = 1, meta_env%n_colvar
774 cv => meta_env%metavar(i_c)
775 diff_ss = cv%ss - cv%ss0
776 IF (cv%periodic) THEN
777 ! The difference of a periodic COLVAR is always within [-pi,pi]
778 diff_ss = sign(1.0_dp, asin(sin(diff_ss)))*acos(cos(diff_ss))
779 END IF
780 cv%epot_s = 0.5_dp*cv%lambda*(diff_ss)**2.0_dp
781 cv%ff_s = cv%lambda*(diff_ss)
782
783 meta_env%epot_s = meta_env%epot_s + cv%epot_s
784 meta_env%ekin_s = meta_env%ekin_s + 0.5_dp*cv%mass*cv%vvp**2
785 END DO
786 END IF
787
788 ! write COLVAR file
789 iw = cp_print_key_unit_nr(logger, meta_env%metadyn_section, &
790 "PRINT%COLVAR", extension=".metadynLog")
791 IF (iw > 0) THEN
792 IF (meta_env%extended_lagrange) THEN
793 WRITE (iw, '(f16.8,70f15.8)') meta_env%time*femtoseconds, &
794 (meta_env%metavar(i)%ss0, i=1, meta_env%n_colvar), &
795 (meta_env%metavar(i)%ss, i=1, meta_env%n_colvar), &
796 (meta_env%metavar(i)%ff_s, i=1, meta_env%n_colvar), &
797 (meta_env%metavar(i)%ff_hills, i=1, meta_env%n_colvar), &
798 (meta_env%metavar(i)%ff_walls, i=1, meta_env%n_colvar), &
799 (meta_env%metavar(i)%vvp, i=1, meta_env%n_colvar), &
800 meta_env%epot_s, &
801 meta_env%hills_env%energy, &
802 meta_env%epot_walls, &
803 (meta_env%ekin_s)*2.0_dp/(real(meta_env%n_colvar, kind=dp))*kelvin
804 ELSE
805 WRITE (iw, '(f16.8,40f13.5)') meta_env%time*femtoseconds, &
806 (meta_env%metavar(i)%ss0, i=1, meta_env%n_colvar), &
807 (meta_env%metavar(i)%ff_hills, i=1, meta_env%n_colvar), &
808 (meta_env%metavar(i)%ff_walls, i=1, meta_env%n_colvar), &
809 meta_env%hills_env%energy, &
810 meta_env%epot_walls
811 END IF
812 END IF
813 CALL cp_print_key_finished_output(iw, logger, meta_env%metadyn_section, &
814 "PRINT%COLVAR")
815
816 ! Temperature for COLVAR
817 IF (meta_env%extended_lagrange) THEN
818 temp = meta_env%ekin_s*2.0_dp/(real(meta_env%n_colvar, kind=dp))*kelvin
819 meta_env%avg_temp = (meta_env%avg_temp*real(meta_env%n_steps, kind=dp) + &
820 temp)/real(meta_env%n_steps + 1, kind=dp)
821 iw = cp_print_key_unit_nr(logger, meta_env%metadyn_section, &
822 "PRINT%TEMPERATURE_COLVAR", extension=".metadynLog")
823 IF (iw > 0) THEN
824 WRITE (iw, '(T2,79("-"))')
825 WRITE (iw, '( A,T51,f10.2,T71,f10.2)') ' COLVARS INSTANTANEOUS/AVERAGE TEMPERATURE ', &
826 temp, meta_env%avg_temp
827 WRITE (iw, '(T2,79("-"))')
828 END IF
829 CALL cp_print_key_finished_output(iw, logger, meta_env%metadyn_section, &
830 "PRINT%TEMPERATURE_COLVAR")
831 END IF
832 CALL timestop(handle)
833
834 END SUBROUTINE metadyn_write_colvar
835
836! **************************************************************************************************
837!> \brief Major driver for adding hills and computing forces due to the history
838!> dependent term
839!> \param meta_env ...
840!> \par History
841!> 04.2004 created
842!> 10.2008 Teodoro Laino [tlaino] - University of Zurich
843!> Major rewriting and addition of multiple walkers
844! **************************************************************************************************
845 SUBROUTINE hills(meta_env)
846 TYPE(meta_env_type), POINTER :: meta_env
847
848 CHARACTER(len=*), PARAMETER :: routinen = 'hills'
849
850 INTEGER :: handle, i, i_hills, ih, intermeta_steps, &
851 iter_nr, iw, n_colvar, n_hills_start, &
852 n_step
853 LOGICAL :: force_gauss
854 REAL(kind=dp) :: cut_func, dfunc, diff, dp2, frac, &
855 slow_growth, v_now_here, v_to_fes, &
856 wtww, ww
857 REAL(kind=dp), DIMENSION(:), POINTER :: ddp, denf, diff_ss, local_last_hills, &
858 numf
859 TYPE(cp_logger_type), POINTER :: logger
860 TYPE(hills_env_type), POINTER :: hills_env
861 TYPE(metavar_type), DIMENSION(:), POINTER :: colvars
862 TYPE(multiple_walkers_type), POINTER :: multiple_walkers
863
864 CALL timeset(routinen, handle)
865
866 NULLIFY (hills_env, multiple_walkers, logger, colvars, ddp, local_last_hills)
867 hills_env => meta_env%hills_env
868 logger => cp_get_default_logger()
869 colvars => meta_env%metavar
870 n_colvar = meta_env%n_colvar
871 n_step = meta_env%n_steps
872
873 ! Create a temporary logger level specific for metadynamics
874 CALL cp_add_iter_level(logger%iter_info, "METADYNAMICS")
875 CALL get_meta_iter_level(meta_env, iter_nr)
876 CALL cp_iterate(logger%iter_info, last=.false., iter_nr=iter_nr)
877
878 ! Set-up restart if any
879 IF (meta_env%hills_env%restart) THEN
880 meta_env%hills_env%restart = .false.
881 IF (meta_env%well_tempered) THEN
882 CALL restart_hills(hills_env%ss_history, hills_env%delta_s_history, hills_env%ww_history, &
883 hills_env%ww, hills_env%n_hills, n_colvar, colvars, meta_env%metadyn_section, &
884 invdt_history=hills_env%invdt_history)
885 ELSE
886 CALL restart_hills(hills_env%ss_history, hills_env%delta_s_history, hills_env%ww_history, &
887 hills_env%ww, hills_env%n_hills, n_colvar, colvars, meta_env%metadyn_section)
888 END IF
889 END IF
890
891 ! Proceed with normal calculation
892 intermeta_steps = n_step - hills_env%old_hill_step
893 force_gauss = .false.
894 IF ((hills_env%min_disp > 0.0_dp) .AND. (hills_env%old_hill_number > 0) .AND. &
895 (intermeta_steps >= hills_env%min_nt_hills)) THEN
896 ALLOCATE (ddp(meta_env%n_colvar))
897 ALLOCATE (local_last_hills(meta_env%n_colvar))
898
899 local_last_hills(1:n_colvar) = hills_env%ss_history(1:n_colvar, hills_env%old_hill_number)
900
901 !RG Calculate the displacement
902 dp2 = 0.0_dp
903 DO i = 1, n_colvar
904 ddp(i) = colvars(i)%ss0 - local_last_hills(i)
905 IF (colvars(i)%periodic) THEN
906 ! The difference of a periodic COLVAR is always within [-pi,pi]
907 ddp(i) = sign(1.0_dp, asin(sin(ddp(i))))*acos(cos(ddp(i)))
908 END IF
909 dp2 = dp2 + ddp(i)**2
910 END DO
911 dp2 = sqrt(dp2)
912 IF (dp2 > hills_env%min_disp) THEN
913 force_gauss = .true.
914 iw = cp_print_key_unit_nr(logger, meta_env%metadyn_section, &
915 "PRINT%PROGRAM_RUN_INFO", extension=".metadynLog")
916 IF (iw > 0) THEN
917 WRITE (unit=iw, fmt="(A,F0.6,A,F0.6)") &
918 " METADYN| Threshold value for COLVAR displacement exceeded: ", &
919 dp2, " > ", hills_env%min_disp
920 END IF
921 CALL cp_print_key_finished_output(iw, logger, meta_env%metadyn_section, &
922 "PRINT%PROGRAM_RUN_INFO")
923 END IF
924 DEALLOCATE (ddp)
925 DEALLOCATE (local_last_hills)
926 END IF
927
928 !RG keep into account adaptive hills
929 IF (((modulo(intermeta_steps, hills_env%nt_hills) == 0) .OR. force_gauss) &
930 .AND. (.NOT. meta_env%restart) .AND. (hills_env%nt_hills > 0)) THEN
931 IF (meta_env%do_multiple_walkers) multiple_walkers => meta_env%multiple_walkers
932
933 n_hills_start = hills_env%n_hills
934 ! Add the hill corresponding to this location
935 IF (meta_env%well_tempered) THEN
936 ! Well-Tempered scaling of hills height
937 v_now_here = 0._dp
938 DO ih = 1, hills_env%n_hills
939 dp2 = 0._dp
940 DO i = 1, n_colvar
941 diff = colvars(i)%ss0 - hills_env%ss_history(i, ih)
942 IF (colvars(i)%periodic) THEN
943 ! The difference of a periodic COLVAR is always within [-pi,pi]
944 diff = sign(1.0_dp, asin(sin(diff)))*acos(cos(diff))
945 END IF
946 diff = (diff)/hills_env%delta_s_history(i, ih)
947 dp2 = dp2 + diff**2
948 END DO
949 v_to_fes = 1.0_dp + meta_env%wttemperature*hills_env%invdt_history(ih)
950 v_now_here = v_now_here + hills_env%ww_history(ih)/v_to_fes*exp(-0.5_dp*dp2)
951 END DO
952 wtww = hills_env%ww*exp(-v_now_here*meta_env%invdt)
953 ww = wtww*(1.0_dp + meta_env%wttemperature*meta_env%invdt)
954 CALL add_hill_single(hills_env, colvars, ww, hills_env%n_hills, n_colvar, meta_env%invdt)
955 ELSE
956 CALL add_hill_single(hills_env, colvars, hills_env%ww, hills_env%n_hills, n_colvar)
957 END IF
958 ! Update local n_hills counter
959 IF (meta_env%do_multiple_walkers) multiple_walkers%n_hills_local = multiple_walkers%n_hills_local + 1
960
961 hills_env%old_hill_number = hills_env%n_hills
962 hills_env%old_hill_step = n_step
963
964 ! Update iteration level for printing
965 CALL get_meta_iter_level(meta_env, iter_nr)
966 CALL cp_iterate(logger%iter_info, last=.false., iter_nr=iter_nr)
967
968 ! Print just program_run_info
969 iw = cp_print_key_unit_nr(logger, meta_env%metadyn_section, &
970 "PRINT%PROGRAM_RUN_INFO", extension=".metadynLog")
971 IF (iw > 0) THEN
972 IF (meta_env%do_multiple_walkers) THEN
973 WRITE (iw, '(/,1X,"METADYN|",A,I0,A,I0,A,/)') &
974 ' Global/Local Hills number (', hills_env%n_hills, '/', multiple_walkers%n_hills_local, &
975 ') added.'
976 ELSE
977 WRITE (iw, '(/,1X,"METADYN|",A,I0,A,/)') ' Hills number (', hills_env%n_hills, ') added.'
978 END IF
979 END IF
980 CALL cp_print_key_finished_output(iw, logger, meta_env%metadyn_section, &
981 "PRINT%PROGRAM_RUN_INFO")
982
983 ! Handle Multiple Walkers
984 IF (meta_env%do_multiple_walkers) THEN
985 ! Print Local Hills file if requested
986 iw = cp_print_key_unit_nr(logger, meta_env%metadyn_section, &
987 "PRINT%HILLS", middle_name="LOCAL", extension=".metadynLog")
988 IF (iw > 0) THEN
989 WRITE (iw, '(f12.1,30f13.5)') meta_env%time*femtoseconds, &
990 (hills_env%ss_history(ih, hills_env%n_hills), ih=1, n_colvar), &
991 (hills_env%delta_s_history(ih, hills_env%n_hills), ih=1, n_colvar), &
992 hills_env%ww_history(hills_env%n_hills)
993 END IF
994 CALL cp_print_key_finished_output(iw, logger, meta_env%metadyn_section, &
995 "PRINT%HILLS")
996
997 ! Check the communication buffer of the other walkers
998 CALL synchronize_multiple_walkers(multiple_walkers, hills_env, colvars, &
999 n_colvar, meta_env%metadyn_section)
1000 END IF
1001
1002 ! Print Hills file if requested (for multiple walkers this file includes
1003 ! the Hills coming from all the walkers).
1004 iw = cp_print_key_unit_nr(logger, meta_env%metadyn_section, &
1005 "PRINT%HILLS", extension=".metadynLog")
1006 IF (iw > 0) THEN
1007 DO i_hills = n_hills_start + 1, hills_env%n_hills
1008 WRITE (iw, '(f12.1,30f13.5)') meta_env%time*femtoseconds, &
1009 (hills_env%ss_history(ih, i_hills), ih=1, n_colvar), &
1010 (hills_env%delta_s_history(ih, i_hills), ih=1, n_colvar), &
1011 hills_env%ww_history(i_hills)
1012 END DO
1013 END IF
1014 CALL cp_print_key_finished_output(iw, logger, meta_env%metadyn_section, &
1015 "PRINT%HILLS")
1016 END IF
1017
1018 ! Computes forces due to the hills: history dependent term
1019 ALLOCATE (ddp(meta_env%n_colvar))
1020 ALLOCATE (diff_ss(meta_env%n_colvar))
1021 ALLOCATE (numf(meta_env%n_colvar))
1022 ALLOCATE (denf(meta_env%n_colvar))
1023
1024 hills_env%energy = 0.0_dp
1025 DO ih = 1, n_colvar
1026 colvars(ih)%ff_hills = 0.0_dp
1027 END DO
1028 DO ih = 1, hills_env%n_hills
1029 slow_growth = 1.0_dp
1030 IF (hills_env%slow_growth .AND. (ih == hills_env%n_hills)) THEN
1031 slow_growth = real(n_step - hills_env%old_hill_step, dp)/real(hills_env%nt_hills, dp)
1032 END IF
1033 dp2 = 0._dp
1034 cut_func = 1.0_dp
1035 DO i = 1, n_colvar
1036 diff_ss(i) = colvars(i)%ss0 - hills_env%ss_history(i, ih)
1037 IF (colvars(i)%periodic) THEN
1038 ! The difference of a periodic COLVAR is always within [-pi,pi]
1039 diff_ss(i) = sign(1.0_dp, asin(sin(diff_ss(i))))*acos(cos(diff_ss(i)))
1040 END IF
1041 IF (hills_env%delta_s_history(i, ih) == 0.0_dp) THEN
1042 ! trick: scale = 0 is interpreted as infinitely wide Gaussian hill
1043 ! instead of infinitely narrow. This way one can combine several
1044 ! one-dimensional bias potentials in a multi-dimensional metadyn
1045 ! simulation.
1046 ddp(i) = 0.0_dp
1047 ELSE
1048 ddp(i) = (diff_ss(i))/hills_env%delta_s_history(i, ih)
1049 END IF
1050 dp2 = dp2 + ddp(i)**2
1051 IF (hills_env%tail_cutoff > 0.0_dp) THEN
1052 frac = abs(ddp(i))/hills_env%tail_cutoff
1053 numf(i) = frac**hills_env%p_exp
1054 denf(i) = frac**hills_env%q_exp
1055 cut_func = cut_func*(1.0_dp - numf(i))/(1.0_dp - denf(i))
1056 END IF
1057 END DO
1058 ! ff_hills contains the "force" due to the hills
1059 dfunc = hills_env%ww_history(ih)*exp(-0.5_dp*dp2)*slow_growth
1060 IF (meta_env%well_tempered) dfunc = dfunc/(1.0_dp + meta_env%wttemperature*hills_env%invdt_history(ih))
1061 hills_env%energy = hills_env%energy + dfunc*cut_func
1062 DO i = 1, n_colvar
1063 IF (hills_env%delta_s_history(i, ih) /= 0.0_dp) THEN
1064 ! only apply a force when the Gaussian hill has a finite width in
1065 ! this direction
1066 colvars(i)%ff_hills = colvars(i)%ff_hills + &
1067 ddp(i)/hills_env%delta_s_history(i, ih)*dfunc*cut_func
1068 IF (hills_env%tail_cutoff > 0.0_dp .AND. abs(diff_ss(i)) > 10.e-5_dp) THEN
1069 colvars(i)%ff_hills = colvars(i)%ff_hills + &
1070 (hills_env%p_exp*numf(i)/(1.0_dp - numf(i)) - hills_env%q_exp*denf(i)/(1.0_dp - denf(i)))* &
1071 dfunc*cut_func/abs(diff_ss(i))
1072 END IF
1073 END IF
1074 END DO
1075 END DO
1076 DEALLOCATE (ddp)
1077 DEALLOCATE (diff_ss)
1078 DEALLOCATE (numf)
1079 DEALLOCATE (denf)
1080
1081 CALL cp_rm_iter_level(logger%iter_info, "METADYNAMICS")
1082
1083 CALL timestop(handle)
1084
1085 END SUBROUTINE hills
1086
1087END MODULE metadynamics
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:126
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:21
bibliography::vandencic2006
integer, save, public vandencic2006
Definition bibliography.F:36
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
colvar_methods
defines collective variables s({R}) and the derivative of this variable wrt R these can then be used ...
Definition colvar_methods.F:16
colvar_methods::colvar_eval_glob_f
subroutine, public colvar_eval_glob_f(icolvar, force_env)
evaluates the derivatives (dsdr) given and due to the given colvar
Definition colvar_methods.F:1763
colvar_types
Initialize the collective variables types.
Definition colvar_types.F:15
colvar_types::torsion_colvar_id
integer, parameter, public torsion_colvar_id
Definition colvar_types.F:56
constraint_fxd
Definition constraint_fxd.F:12
constraint_fxd::fix_atom_control
subroutine, public fix_atom_control(force_env, w)
allows for fix atom constraints
Definition constraint_fxd.F:65
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:838
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1058
cp_output_handling::cp_iterate
subroutine, public cp_iterate(iteration_info, last, iter_nr, increment, iter_nr_out)
adds one to the actual iteration
Definition cp_output_handling.F:559
cp_output_handling::cp_rm_iter_level
subroutine, public cp_rm_iter_level(iteration_info, level_name, n_rlevel_att)
Removes an iteration level.
Definition cp_output_handling.F:690
cp_output_handling::cp_add_iter_level
subroutine, public cp_add_iter_level(iteration_info, level_name, n_rlevel_new)
Adds an iteration level.
Definition cp_output_handling.F:647
cp_subsys_types
types that represent a subsys, i.e. a part of the system
Definition cp_subsys_types.F:16
cp_subsys_types::cp_subsys_get
subroutine, public cp_subsys_get(subsys, ref_count, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell)
returns information about various attributes of the given subsys
Definition cp_subsys_types.F:364
force_env_types
Interface for the force calculations.
Definition force_env_types.F:18
force_env_types::force_env_get
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env, ipi_env)
returns various attributes about the force environment
Definition force_env_types.F:344
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::do_wall_p
integer, parameter, public do_wall_p
Definition input_constants.F:393
input_constants::do_wall_m
integer, parameter, public do_wall_m
Definition input_constants.F:393
input_constants::do_wall_reflective
integer, parameter, public do_wall_reflective
Definition input_constants.F:387
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:735
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition input_section_types.F:708
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1052
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
message_passing::cp2k_is_parallel
logical, parameter, public cp2k_is_parallel
Definition message_passing.F:114
metadynamics_types
defines types for metadynamics calculation
Definition metadynamics_types.F:13
metadynamics_utils
Performs the metadynamics calculation.
Definition metadynamics_utils.F:14
metadynamics_utils::synchronize_multiple_walkers
subroutine, public synchronize_multiple_walkers(multiple_walkers, hills_env, colvars, n_colvar, metadyn_section)
Synchronize with the rest of the walkers.
Definition metadynamics_utils.F:516
metadynamics_utils::get_meta_iter_level
subroutine, public get_meta_iter_level(meta_env, iter_nr)
Retrieves the iteration level for the metadynamics loop.
Definition metadynamics_utils.F:872
metadynamics_utils::meta_walls
subroutine, public meta_walls(meta_env)
...
Definition metadynamics_utils.F:891
metadynamics_utils::restart_hills
subroutine, public restart_hills(ss_history, delta_s_history, ww_history, ww, n_hills, n_colvar, colvars, metadyn_section, invdt_history)
Restart Hills Information.
Definition metadynamics_utils.F:742
metadynamics_utils::add_hill_single
subroutine, public add_hill_single(hills_env, colvars, ww, n_hills, n_colvar, invdt)
Add a single Hill.
Definition metadynamics_utils.F:662
metadynamics
Performs the metadynamics calculation.
Definition metadynamics.F:14
metadynamics::metadyn_finalise_plumed
subroutine, public metadyn_finalise_plumed()
makes sure PLUMED is shut down cleanly
Definition metadynamics.F:209
metadynamics::metadyn_initialise_plumed
subroutine, public metadyn_initialise_plumed(force_env, simpar, itimes)
...
Definition metadynamics.F:133
metadynamics::metadyn_forces
subroutine, public metadyn_forces(force_env, vel)
add forces to the subsys due to the metadynamics run possibly modifies the velocites (if reflective w...
Definition metadynamics.F:387
metadynamics::metadyn_velocities_colvar
subroutine, public metadyn_velocities_colvar(force_env, rand)
Evolves velocities COLVAR according to Vanden-Eijnden Ciccotti C.Phys.Letter 429 (2006) 310-316.
Definition metadynamics.F:649
metadynamics::metadyn_integrator
subroutine, public metadyn_integrator(force_env, itimes, vel, rand)
General driver for applying metadynamics.
Definition metadynamics.F:236
metadynamics::metadyn_write_colvar
subroutine, public metadyn_write_colvar(force_env)
Write down COLVAR information evolved according to Vanden-Eijnden Ciccotti C.Phys....
Definition metadynamics.F:748
particle_list_types
represent a simple array based list of the given type
Definition particle_list_types.F:15
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::boltzmann
real(kind=dp), parameter, public boltzmann
Definition physcon.F:129
physcon::femtoseconds
real(kind=dp), parameter, public femtoseconds
Definition physcon.F:153
physcon::joule
real(kind=dp), parameter, public joule
Definition physcon.F:159
physcon::kelvin
real(kind=dp), parameter, public kelvin
Definition physcon.F:165
physcon::angstrom
real(kind=dp), parameter, public angstrom
Definition physcon.F:144
physcon::picoseconds
real(kind=dp), parameter, public picoseconds
Definition physcon.F:156
physcon::kjmol
real(kind=dp), parameter, public kjmol
Definition physcon.F:168
reference_manager
provides a uniform framework to add references to CP2K cite and output these
Definition reference_manager.F:25
reference_manager::cite_reference
subroutine, public cite_reference(key)
marks a given reference as cited.
Definition reference_manager.F:85
simpar_types
Type for storing MD parameters.
Definition simpar_types.F:14
virial_types
Definition virial_types.F:13
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
colvar_types::colvar_p_type
Definition colvar_types.F:393
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
cp_subsys_types::cp_subsys_type
represents a system: atoms, molecules, their pos,vel,...
Definition cp_subsys_types.F:89
force_env_types::force_env_type
wrapper to abstract the force evaluation of the various methods
Definition force_env_types.F:145
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
metadynamics_types::hills_env_type
defines types for HILLS
Definition metadynamics_types.F:40
metadynamics_types::meta_env_type
defines meta_env type
Definition metadynamics_types.F:95
metadynamics_types::metavar_type
defines types for COLVAR used in the metadynamics
Definition metadynamics_types.F:69
metadynamics_types::multiple_walkers_type
defines types for multiple walkers run
Definition metadynamics_types.F:82
particle_list_types::particle_list_type
represent a list of objects
Definition particle_list_types.F:46
simpar_types::simpar_type
Simulation parameter type for molecular dynamics.
Definition simpar_types.F:30
virial_types::virial_type
Definition virial_types.F:26