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