(git:ed6f26b)
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-2025 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#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
120CONTAINS
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
1069END 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:120
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:1725
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:853
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:1063
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:705
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:662
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, ipi_env)
returns various attributes about the force environment
Definition force_env_types.F:342
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:373
input_constants::do_wall_m
integer, parameter, public do_wall_m
Definition input_constants.F:373
input_constants::do_wall_reflective
integer, parameter, public do_wall_reflective
Definition input_constants.F:367
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:731
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:704
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:1047
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: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
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
colvar_types::colvar_p_type
Definition colvar_types.F:383
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