Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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 2 : 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 2 : CALL timeset(routineN, handle)
150 2 : CPASSERT(ASSOCIATED(force_env))
151 2 : CPASSERT(ASSOCIATED(simpar))
152 :
153 : #if defined (__PLUMED2)
154 2 : NULLIFY (cell, subsys)
155 2 : CALL force_env_get(force_env, subsys=subsys, cell=cell)
156 2 : CALL pbc_cp2k_plumed_getset_cell(cell, set=.TRUE.) !Store the cell pointer for later use.
157 2 : timestep_plumed = simpar%dt
158 2 : 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 2 : plumedavailable = plumed_installed()
167 :
168 2 : IF (plumedavailable <= 0) THEN
169 0 : CPABORT("NO PLUMED IN YOUR LD_LIBRARY_PATH?")
170 : ELSE
171 2 : CALL plumed_gcreate()
172 2 : IF (cp2k_is_parallel) CALL plumed_gcmd_int("setMPIFComm"//CHAR(0), force_env%para_env%get_handle())
173 2 : CALL plumed_gcmd_int("setRealPrecision"//CHAR(0), 8)
174 2 : CALL plumed_gcmd_double("setMDEnergyUnits"//CHAR(0), kjmol) ! Hartree to kjoule/mol
175 2 : CALL plumed_gcmd_double("setMDLengthUnits"//CHAR(0), angstrom*0.1_dp) ! Bohr to nm
176 2 : CALL plumed_gcmd_double("setMDTimeUnits"//CHAR(0), picoseconds) ! atomic units to ps
177 2 : CALL plumed_gcmd_char("setPlumedDat"//CHAR(0), force_env%meta_env%plumed_input_file//CHAR(0))
178 2 : CALL plumed_gcmd_char("setLogFile"//CHAR(0), "PLUMED.OUT"//CHAR(0))
179 2 : CALL plumed_gcmd_int("setNatoms"//CHAR(0), natom_plumed)
180 2 : CALL plumed_gcmd_char("setMDEngine"//CHAR(0), "cp2k"//CHAR(0))
181 2 : CALL plumed_gcmd_double("setTimestep"//CHAR(0), timestep_plumed)
182 2 : 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 2 : CALL timestop(handle)
193 :
194 2 : END SUBROUTINE metadyn_initialise_plumed
195 :
196 : ! **************************************************************************************************
197 : !> \brief makes sure PLUMED is shut down cleanly
198 : ! **************************************************************************************************
199 2 : SUBROUTINE metadyn_finalise_plumed()
200 :
201 : CHARACTER(LEN=*), PARAMETER :: routineN = 'metadyn_finalise_plumed'
202 :
203 : INTEGER :: handle
204 :
205 2 : CALL timeset(routineN, handle)
206 :
207 : #if defined (__PLUMED2)
208 2 : CALL plumed_gfinalize()
209 : #endif
210 :
211 2 : CALL timestop(handle)
212 :
213 2 : END SUBROUTINE metadyn_finalise_plumed
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 41575 : 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 41575 : 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 41575 : DIMENSION(:) :: pos_plumed_x, pos_plumed_y, pos_plumed_z
249 : REAL(KIND=dp), ALLOCATABLE, &
250 41575 : DIMENSION(:) :: force_plumed_x, force_plumed_y, force_plumed_z
251 : #endif
252 :
253 41575 : CALL timeset(routineN, handle)
254 :
255 : ! Apply Metadynamics
256 41575 : IF (ASSOCIATED(force_env%meta_env)) THEN
257 13794 : IF (force_env%meta_env%use_plumed .EQV. .TRUE.) THEN
258 10 : plumed_itimes = itimes
259 : #if defined (__PLUMED2)
260 10 : CALL force_env_get(force_env, meta_env=meta_env, subsys=subsys, cell=cell)
261 10 : natom_plumed = subsys%particles%n_els
262 30 : ALLOCATE (pos_plumed_x(natom_plumed))
263 20 : ALLOCATE (pos_plumed_y(natom_plumed))
264 20 : ALLOCATE (pos_plumed_z(natom_plumed))
265 :
266 20 : ALLOCATE (force_plumed_x(natom_plumed))
267 20 : ALLOCATE (force_plumed_y(natom_plumed))
268 20 : ALLOCATE (force_plumed_z(natom_plumed))
269 :
270 20 : ALLOCATE (mass_plumed(natom_plumed))
271 :
272 100 : force_plumed_x(:) = 0.0d0
273 100 : force_plumed_y(:) = 0.0d0
274 100 : force_plumed_z(:) = 0.0d0
275 :
276 100 : DO i_part = 1, natom_plumed
277 90 : pos_plumed_x(i_part) = subsys%particles%els(i_part)%r(1)
278 90 : pos_plumed_y(i_part) = subsys%particles%els(i_part)%r(2)
279 90 : pos_plumed_z(i_part) = subsys%particles%els(i_part)%r(3)
280 100 : mass_plumed(i_part) = subsys%particles%els(i_part)%atomic_kind%mass
281 : END DO
282 :
283 : ! stupid cell conversion, to be fixed
284 10 : celltbt(1, 1) = cell%hmat(1, 1)
285 10 : celltbt(1, 2) = cell%hmat(1, 2)
286 10 : celltbt(1, 3) = cell%hmat(1, 3)
287 10 : celltbt(2, 1) = cell%hmat(2, 1)
288 10 : celltbt(2, 2) = cell%hmat(2, 2)
289 10 : celltbt(2, 3) = cell%hmat(2, 3)
290 10 : celltbt(3, 1) = cell%hmat(3, 1)
291 10 : celltbt(3, 2) = cell%hmat(3, 2)
292 10 : celltbt(3, 3) = cell%hmat(3, 3)
293 :
294 : ! virial
295 10 : plu_vir(:, :) = 0.0d0
296 :
297 10 : CALL force_env_get(force_env, potential_energy=stpcfg)
298 :
299 10 : CALL plumed_gcmd_int("setStep"//CHAR(0), plumed_itimes)
300 10 : CALL plumed_gcmd_doubles("setPositionsX"//CHAR(0), pos_plumed_x(:))
301 10 : CALL plumed_gcmd_doubles("setPositionsY"//CHAR(0), pos_plumed_y(:))
302 10 : CALL plumed_gcmd_doubles("setPositionsZ"//CHAR(0), pos_plumed_z(:))
303 10 : CALL plumed_gcmd_doubles("setMasses"//CHAR(0), mass_plumed(:))
304 10 : CALL plumed_gcmd_doubles("setBox"//CHAR(0), celltbt)
305 10 : CALL plumed_gcmd_doubles("setVirial"//CHAR(0), plu_vir) ! PluMed Output, NOT ACTUALLY USED !
306 10 : CALL plumed_gcmd_double("setEnergy"//CHAR(0), stpcfg)
307 10 : CALL plumed_gcmd_doubles("setForcesX"//CHAR(0), force_plumed_x(:)) ! PluMed Output !
308 10 : CALL plumed_gcmd_doubles("setForcesY"//CHAR(0), force_plumed_y(:)) ! PluMed Output !
309 10 : CALL plumed_gcmd_doubles("setForcesZ"//CHAR(0), force_plumed_z(:)) ! PluMed Output !
310 :
311 10 : CALL plumed_gcmd_int("prepareCalc"//CHAR(0), 0)
312 10 : CALL plumed_gcmd_int("prepareDependencies"//CHAR(0), 0)
313 10 : CALL plumed_gcmd_int("shareData"//CHAR(0), 0)
314 :
315 10 : CALL plumed_gcmd_int("performCalc"//CHAR(0), 0)
316 :
317 100 : DO i_part = 1, natom_plumed
318 90 : subsys%particles%els(i_part)%f(1) = subsys%particles%els(i_part)%f(1) + force_plumed_x(i_part)
319 90 : subsys%particles%els(i_part)%f(2) = subsys%particles%els(i_part)%f(2) + force_plumed_y(i_part)
320 100 : subsys%particles%els(i_part)%f(3) = subsys%particles%els(i_part)%f(3) + force_plumed_z(i_part)
321 : END DO
322 :
323 10 : DEALLOCATE (pos_plumed_x, pos_plumed_y, pos_plumed_z)
324 10 : DEALLOCATE (force_plumed_x, force_plumed_y, force_plumed_z)
325 10 : DEALLOCATE (mass_plumed)
326 :
327 : ! Constraints ONLY of Fixed Atom type
328 20 : 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 13784 : IF (force_env%meta_env%langevin) THEN
339 240 : IF (.NOT. PRESENT(rand)) THEN
340 0 : CPABORT("Langevin on COLVAR not implemented for this MD ensemble!")
341 : END IF
342 : ! *** Velocity Verlet for Langevin S(t)->S(t+1)
343 240 : CALL metadyn_position_colvar(force_env)
344 : ! *** Forces from Vs and S(X(t+1))
345 240 : CALL metadyn_forces(force_env)
346 : ! *** Velocity Verlet for Langeving *** v(t+1/2)--> v(t)
347 240 : CALL metadyn_velocities_colvar(force_env, rand)
348 : ELSE
349 13668 : CALL metadyn_forces(force_env, vel)
350 : END IF
351 : ! *** Write down COVAR informations
352 13784 : CALL metadyn_write_colvar(force_env)
353 : END IF
354 : END IF
355 :
356 41575 : CALL timestop(handle)
357 :
358 41575 : 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 13928 : 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 13928 : 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 13928 : NULLIFY (logger, meta_env)
389 13928 : meta_env => force_env%meta_env
390 0 : IF (.NOT. ASSOCIATED(meta_env)) RETURN
391 :
392 13928 : CALL timeset(routineN, handle)
393 13928 : logger => cp_get_default_logger()
394 13928 : NULLIFY (colvar_p, subsys, cv, ss0_section, vvp_section)
395 13928 : CALL force_env_get(force_env, subsys=subsys)
396 :
397 13928 : dt = meta_env%dt
398 13928 : IF (.NOT. meta_env%restart) meta_env%n_steps = meta_env%n_steps + 1
399 :
400 : ! Initialize velocity
401 13928 : IF (meta_env%restart .AND. meta_env%extended_lagrange) THEN
402 54 : meta_env%ekin_s = 0.0_dp
403 130 : DO i_c = 1, meta_env%n_colvar
404 76 : cv => meta_env%metavar(i_c)
405 76 : cv%vvp = force_env%globenv%gaussian_rng_stream%next()
406 130 : meta_env%ekin_s = meta_env%ekin_s + 0.5_dp*cv%mass*cv%vvp**2
407 : END DO
408 54 : ekin_w = 0.5_dp*meta_env%temp_wanted*REAL(meta_env%n_colvar, KIND=dp)
409 54 : fac_t = SQRT(ekin_w/MAX(meta_env%ekin_s, 1.0E-8_dp))
410 130 : DO i_c = 1, meta_env%n_colvar
411 76 : cv => meta_env%metavar(i_c)
412 130 : cv%vvp = cv%vvp*fac_t
413 : END DO
414 54 : 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 28600 : DO i_c = 1, meta_env%n_colvar
420 14672 : cv => meta_env%metavar(i_c)
421 14672 : icolvar = cv%icolvar
422 14672 : CALL colvar_eval_glob_f(icolvar, force_env)
423 14672 : cv%ss = subsys%colvar_p(icolvar)%colvar%ss
424 :
425 : ! Setup the periodic flag if the COLVAR is (-pi,pi] periodic
426 14672 : cv%periodic = (subsys%colvar_p(icolvar)%colvar%type_id == torsion_colvar_id)
427 :
428 : ! Restart for Extended Lagrangian Metadynamics
429 14672 : IF (meta_env%restart) THEN
430 : ! Initialize the position of the collective variable in the extended lagrange
431 190 : ss0_section => section_vals_get_subs_vals(meta_env%metadyn_section, "EXT_LAGRANGE_SS0")
432 190 : CALL section_vals_get(ss0_section, explicit=explicit)
433 190 : IF (explicit) THEN
434 : CALL section_vals_val_get(ss0_section, "_DEFAULT_KEYWORD_", &
435 4 : i_rep_val=i_c, r_val=rval)
436 4 : cv%ss0 = rval
437 : ELSE
438 186 : cv%ss0 = cv%ss
439 : END IF
440 190 : vvp_section => section_vals_get_subs_vals(meta_env%metadyn_section, "EXT_LAGRANGE_VVP")
441 190 : CALL section_vals_get(vvp_section, explicit=explicit)
442 190 : IF (explicit) THEN
443 : CALL section_vals_val_get(vvp_section, "_DEFAULT_KEYWORD_", &
444 4 : i_rep_val=i_c, r_val=rval)
445 4 : cv%vvp = rval
446 : END IF
447 : END IF
448 : !
449 28600 : IF (.NOT. meta_env%extended_lagrange) THEN
450 12506 : cv%ss0 = cv%ss
451 12506 : cv%vvp = 0.0_dp
452 : END IF
453 : END DO
454 : ! History dependent forces (evaluated at s0)
455 13928 : IF (meta_env%do_hills) CALL hills(meta_env)
456 :
457 : ! Apply walls to the colvars
458 13928 : CALL meta_walls(meta_env)
459 :
460 13928 : meta_env%restart = .FALSE.
461 13928 : IF (.NOT. meta_env%extended_lagrange) THEN
462 12312 : meta_env%ekin_s = 0.0_dp
463 12312 : meta_env%epot_s = 0.0_dp
464 12312 : meta_env%epot_walls = 0.0_dp
465 24818 : DO i_c = 1, meta_env%n_colvar
466 12506 : cv => meta_env%metavar(i_c)
467 12506 : cv%epot_s = 0.0_dp
468 12506 : cv%ff_s = 0.0_dp
469 12506 : meta_env%epot_walls = meta_env%epot_walls + cv%epot_walls
470 12506 : icolvar = cv%icolvar
471 12506 : NULLIFY (particles)
472 : CALL cp_subsys_get(subsys, colvar_p=colvar_p, &
473 12506 : particles=particles)
474 63118 : DO ii = 1, colvar_p(icolvar)%colvar%n_atom_s
475 38300 : i = colvar_p(icolvar)%colvar%i_atom(ii)
476 38300 : fft = cv%ff_hills + cv%ff_walls
477 280606 : particles%els(i)%f = particles%els(i)%f + fft*colvar_p(icolvar)%colvar%dsdr(:, ii)
478 : END DO
479 : END DO
480 : ELSE
481 1616 : meta_env%ekin_s = 0.0_dp
482 1616 : meta_env%epot_s = 0.0_dp
483 1616 : meta_env%epot_walls = 0.0_dp
484 3782 : DO i_c = 1, meta_env%n_colvar
485 2166 : cv => meta_env%metavar(i_c)
486 2166 : diff_ss = cv%ss - cv%ss0
487 2166 : IF (cv%periodic) THEN
488 : ! The difference of a periodic COLVAR is always within [-pi,pi]
489 0 : diff_ss = SIGN(1.0_dp, ASIN(SIN(diff_ss)))*ACOS(COS(diff_ss))
490 : END IF
491 2166 : cv%epot_s = 0.5_dp*cv%lambda*(diff_ss)**2.0_dp
492 2166 : cv%ff_s = cv%lambda*(diff_ss)
493 2166 : icolvar = cv%icolvar
494 : ! forces on the atoms
495 2166 : NULLIFY (particles)
496 : CALL cp_subsys_get(subsys, colvar_p=colvar_p, &
497 2166 : particles=particles)
498 10944 : DO ii = 1, colvar_p(icolvar)%colvar%n_atom_s
499 8778 : i = colvar_p(icolvar)%colvar%i_atom(ii)
500 63612 : 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 3782 : IF (.NOT. meta_env%langevin) THEN
504 1678 : fft = cv%ff_s + cv%ff_hills + cv%ff_walls
505 1678 : cv%vvp = cv%vvp + dt*fft/cv%mass
506 1678 : meta_env%ekin_s = meta_env%ekin_s + 0.5_dp*cv%mass*cv%vvp**2
507 1678 : meta_env%epot_s = meta_env%epot_s + cv%epot_s
508 1678 : meta_env%epot_walls = meta_env%epot_walls + cv%epot_walls
509 : END IF
510 : END DO
511 : ! velocity rescaling on the s0
512 1616 : IF (meta_env%tempcontrol .AND. (.NOT. meta_env%langevin)) THEN
513 604 : ekin_w = 0.5_dp*meta_env%temp_wanted*REAL(meta_env%n_colvar, KIND=dp)
514 604 : tol_ekin = 0.5_dp*meta_env%toll_temp*REAL(meta_env%n_colvar, KIND=dp)
515 604 : IF (ABS(ekin_w - meta_env%ekin_s) > tol_ekin) THEN
516 354 : fac_t = SQRT(ekin_w/MAX(meta_env%ekin_s, 1.0E-8_dp))
517 772 : DO i_c = 1, meta_env%n_colvar
518 418 : cv => meta_env%metavar(i_c)
519 772 : cv%vvp = cv%vvp*fac_t
520 : END DO
521 354 : meta_env%ekin_s = ekin_w
522 : END IF
523 : END IF
524 : ! Reflective Wall only for s0
525 3782 : DO i_c = 1, meta_env%n_colvar
526 2166 : cv => meta_env%metavar(i_c)
527 3782 : IF (cv%do_wall) THEN
528 700 : DO iwall = 1, SIZE(cv%walls)
529 248 : SELECT CASE (cv%walls(iwall)%id_type)
530 : CASE (do_wall_reflective)
531 204 : ss0_test = cv%ss0 + dt*cv%vvp
532 204 : IF (cv%periodic) THEN
533 : ! A periodic COLVAR is always within [-pi,pi]
534 0 : ss0_test = SIGN(1.0_dp, ASIN(SIN(ss0_test)))*ACOS(COS(ss0_test))
535 : END IF
536 452 : SELECT CASE (cv%walls(iwall)%id_direction)
537 : CASE (do_wall_p)
538 102 : IF ((ss0_test > cv%walls(iwall)%pos) .AND. (cv%vvp > 0)) cv%vvp = -cv%vvp
539 : CASE (do_wall_m)
540 204 : 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 1616 : IF (.NOT. meta_env%langevin) THEN
548 3050 : DO i_c = 1, meta_env%n_colvar
549 1678 : cv => meta_env%metavar(i_c)
550 1678 : cv%ss0 = cv%ss0 + dt*cv%vvp
551 3050 : IF (cv%periodic) THEN
552 : ! A periodic COLVAR is always within [-pi,pi]
553 0 : 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 13928 : CALL fix_atom_control(force_env)
560 :
561 : ! Reflective Wall only for ss
562 28600 : DO i_c = 1, meta_env%n_colvar
563 14672 : cv => meta_env%metavar(i_c)
564 28600 : IF (cv%do_wall) THEN
565 23162 : DO iwall = 1, SIZE(cv%walls)
566 11276 : SELECT CASE (cv%walls(iwall)%id_type)
567 : CASE (do_wall_reflective)
568 612 : SELECT CASE (cv%walls(iwall)%id_direction)
569 : CASE (do_wall_p)
570 204 : IF (cv%ss < cv%walls(iwall)%pos) CYCLE
571 204 : check_val = -1.0_dp
572 : CASE (do_wall_m)
573 204 : IF (cv%ss > cv%walls(iwall)%pos) CYCLE
574 408 : check_val = 1.0_dp
575 : END SELECT
576 2 : NULLIFY (particles)
577 2 : icolvar = cv%icolvar
578 2 : CALL cp_subsys_get(subsys, colvar_p=colvar_p, particles=particles)
579 2 : scal = 0.0_dp
580 2 : scalf = 0.0_dp
581 2 : norm = 0.0_dp
582 10 : DO ii = 1, colvar_p(icolvar)%colvar%n_atom_s
583 8 : i = colvar_p(icolvar)%colvar%i_atom(ii)
584 8 : IF (PRESENT(vel)) THEN
585 8 : scal = scal + vel(1, i)*colvar_p(icolvar)%colvar%dsdr(1, ii)
586 8 : scal = scal + vel(2, i)*colvar_p(icolvar)%colvar%dsdr(2, ii)
587 8 : scal = scal + vel(3, i)*colvar_p(icolvar)%colvar%dsdr(3, ii)
588 : ELSE
589 0 : scal = scal + particles%els(i)%v(1)*colvar_p(icolvar)%colvar%dsdr(1, ii)
590 0 : scal = scal + particles%els(i)%v(2)*colvar_p(icolvar)%colvar%dsdr(2, ii)
591 0 : scal = scal + particles%els(i)%v(3)*colvar_p(icolvar)%colvar%dsdr(3, ii)
592 : END IF
593 8 : scalf = scalf + particles%els(i)%f(1)*colvar_p(icolvar)%colvar%dsdr(1, ii)
594 8 : scalf = scalf + particles%els(i)%f(2)*colvar_p(icolvar)%colvar%dsdr(2, ii)
595 8 : scalf = scalf + particles%els(i)%f(3)*colvar_p(icolvar)%colvar%dsdr(3, ii)
596 8 : norm = norm + colvar_p(icolvar)%colvar%dsdr(1, ii)**2
597 8 : norm = norm + colvar_p(icolvar)%colvar%dsdr(2, ii)**2
598 10 : norm = norm + colvar_p(icolvar)%colvar%dsdr(3, ii)**2
599 : END DO
600 2 : IF (norm /= 0.0_dp) scal = scal/norm
601 2 : IF (norm /= 0.0_dp) scalf = scalf/norm
602 :
603 2 : IF (scal*check_val > 0.0_dp) CYCLE
604 11896 : DO ii = 1, colvar_p(icolvar)%colvar%n_atom_s
605 8 : i = colvar_p(icolvar)%colvar%i_atom(ii)
606 8 : IF (PRESENT(vel)) THEN
607 32 : vel(:, i) = vel(:, i) - 2.0_dp*colvar_p(icolvar)%colvar%dsdr(:, ii)*scal
608 : ELSE
609 0 : 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 58 : 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 13928 : CALL timestop(handle)
620 13928 : 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 480 : 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 480 : NULLIFY (logger, meta_env, cv)
644 480 : meta_env => force_env%meta_env
645 480 : IF (.NOT. ASSOCIATED(meta_env)) RETURN
646 :
647 480 : CALL timeset(routineN, handle)
648 480 : logger => cp_get_default_logger()
649 : ! Add citation
650 480 : IF (meta_env%langevin) CALL cite_reference(VandenCic2006)
651 :
652 480 : dt = meta_env%dt
653 : ! History dependent forces (evaluated at s0)
654 480 : IF (meta_env%do_hills) CALL hills(meta_env)
655 :
656 : ! Evolve Velocities
657 480 : meta_env%ekin_s = 0.0_dp
658 480 : meta_env%epot_walls = 0.0_dp
659 1440 : DO i_c = 1, meta_env%n_colvar
660 960 : cv => meta_env%metavar(i_c)
661 960 : diff_ss = cv%ss - cv%ss0
662 960 : IF (cv%periodic) THEN
663 : ! The difference of a periodic COLVAR is always within [-pi,pi]
664 0 : diff_ss = SIGN(1.0_dp, ASIN(SIN(diff_ss)))*ACOS(COS(diff_ss))
665 : END IF
666 960 : cv%epot_s = 0.5_dp*cv%lambda*(diff_ss)**2.0_dp
667 960 : cv%ff_s = cv%lambda*(diff_ss)
668 :
669 960 : fft = cv%ff_s + cv%ff_hills
670 960 : 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 960 : 0.5_dp*SQRT(dt)*sigma*rand(i_c)
673 960 : meta_env%ekin_s = meta_env%ekin_s + 0.5_dp*cv%mass*cv%vvp**2
674 1440 : meta_env%epot_walls = meta_env%epot_walls + cv%epot_walls
675 : END DO
676 480 : 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 480 : 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 240 : NULLIFY (logger, meta_env, cv)
699 240 : meta_env => force_env%meta_env
700 240 : IF (.NOT. ASSOCIATED(meta_env)) RETURN
701 :
702 240 : CALL timeset(routineN, handle)
703 240 : logger => cp_get_default_logger()
704 :
705 : ! Add citation
706 240 : IF (meta_env%langevin) CALL cite_reference(VandenCic2006)
707 240 : dt = meta_env%dt
708 :
709 : ! Update of ss0
710 720 : DO i_c = 1, meta_env%n_colvar
711 480 : cv => meta_env%metavar(i_c)
712 480 : cv%ss0 = cv%ss0 + dt*cv%vvp
713 720 : IF (cv%periodic) THEN
714 : ! A periodic COLVAR is always within [-pi,pi]
715 0 : cv%ss0 = SIGN(1.0_dp, ASIN(SIN(cv%ss0)))*ACOS(COS(cv%ss0))
716 : END IF
717 : END DO
718 240 : 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 27856 : 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 13928 : NULLIFY (logger, meta_env, cv)
741 13928 : meta_env => force_env%meta_env
742 13928 : IF (.NOT. ASSOCIATED(meta_env)) RETURN
743 :
744 13928 : CALL timeset(routineN, handle)
745 13928 : 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 13928 : IF (meta_env%langevin) THEN
753 244 : meta_env%ekin_s = 0.0_dp
754 244 : meta_env%epot_s = 0.0_dp
755 732 : DO i_c = 1, meta_env%n_colvar
756 488 : cv => meta_env%metavar(i_c)
757 488 : diff_ss = cv%ss - cv%ss0
758 488 : IF (cv%periodic) THEN
759 : ! The difference of a periodic COLVAR is always within [-pi,pi]
760 0 : diff_ss = SIGN(1.0_dp, ASIN(SIN(diff_ss)))*ACOS(COS(diff_ss))
761 : END IF
762 488 : cv%epot_s = 0.5_dp*cv%lambda*(diff_ss)**2.0_dp
763 488 : cv%ff_s = cv%lambda*(diff_ss)
764 :
765 488 : meta_env%epot_s = meta_env%epot_s + cv%epot_s
766 732 : 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 13928 : "PRINT%COLVAR", extension=".metadynLog")
773 13928 : IF (iw > 0) THEN
774 6682 : IF (meta_env%extended_lagrange) THEN
775 781 : WRITE (iw, '(f16.8,70f15.8)') meta_env%time*femtoseconds, &
776 2591 : (meta_env%metavar(i)%ss0, i=1, meta_env%n_colvar), &
777 2591 : (meta_env%metavar(i)%ss, i=1, meta_env%n_colvar), &
778 2591 : (meta_env%metavar(i)%ff_s, i=1, meta_env%n_colvar), &
779 2591 : (meta_env%metavar(i)%ff_hills, i=1, meta_env%n_colvar), &
780 2591 : (meta_env%metavar(i)%ff_walls, i=1, meta_env%n_colvar), &
781 2591 : (meta_env%metavar(i)%vvp, i=1, meta_env%n_colvar), &
782 781 : meta_env%epot_s, &
783 781 : meta_env%hills_env%energy, &
784 781 : meta_env%epot_walls, &
785 12422 : (meta_env%ekin_s)*2.0_dp/(REAL(meta_env%n_colvar, KIND=dp))*kelvin
786 : ELSE
787 5901 : WRITE (iw, '(f16.8,40f13.5)') meta_env%time*femtoseconds, &
788 17800 : (meta_env%metavar(i)%ss0, i=1, meta_env%n_colvar), &
789 17800 : (meta_env%metavar(i)%ff_hills, i=1, meta_env%n_colvar), &
790 17800 : (meta_env%metavar(i)%ff_walls, i=1, meta_env%n_colvar), &
791 5901 : meta_env%hills_env%energy, &
792 47499 : meta_env%epot_walls
793 : END IF
794 : END IF
795 : CALL cp_print_key_finished_output(iw, logger, meta_env%metadyn_section, &
796 13928 : "PRINT%COLVAR")
797 :
798 : ! Temperature for COLVAR
799 13928 : IF (meta_env%extended_lagrange) THEN
800 1616 : 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 1616 : temp)/REAL(meta_env%n_steps + 1, KIND=dp)
803 : iw = cp_print_key_unit_nr(logger, meta_env%metadyn_section, &
804 1616 : "PRINT%TEMPERATURE_COLVAR", extension=".metadynLog")
805 1616 : IF (iw > 0) THEN
806 808 : WRITE (iw, '(T2,79("-"))')
807 808 : WRITE (iw, '( A,T51,f10.2,T71,f10.2)') ' COLVARS INSTANTANEOUS/AVERAGE TEMPERATURE ', &
808 1616 : temp, meta_env%avg_temp
809 808 : WRITE (iw, '(T2,79("-"))')
810 : END IF
811 : CALL cp_print_key_finished_output(iw, logger, meta_env%metadyn_section, &
812 1616 : "PRINT%TEMPERATURE_COLVAR")
813 : END IF
814 13928 : 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 11278 : 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 11278 : REAL(KIND=dp), DIMENSION(:), POINTER :: ddp, denf, diff_ss, local_last_hills, &
840 11278 : numf
841 : TYPE(cp_logger_type), POINTER :: logger
842 : TYPE(hills_env_type), POINTER :: hills_env
843 11278 : TYPE(metavar_type), DIMENSION(:), POINTER :: colvars
844 : TYPE(multiple_walkers_type), POINTER :: multiple_walkers
845 :
846 11278 : CALL timeset(routineN, handle)
847 :
848 11278 : NULLIFY (hills_env, multiple_walkers, logger, colvars, ddp, local_last_hills)
849 11278 : hills_env => meta_env%hills_env
850 11278 : logger => cp_get_default_logger()
851 11278 : colvars => meta_env%metavar
852 11278 : n_colvar = meta_env%n_colvar
853 11278 : n_step = meta_env%n_steps
854 :
855 : ! Create a temporary logger level specific for metadynamics
856 11278 : CALL cp_add_iter_level(logger%iter_info, "METADYNAMICS")
857 11278 : CALL get_meta_iter_level(meta_env, iter_nr)
858 11278 : CALL cp_iterate(logger%iter_info, last=.FALSE., iter_nr=iter_nr)
859 :
860 : ! Set-up restart if any
861 11278 : IF (meta_env%hills_env%restart) THEN
862 118 : meta_env%hills_env%restart = .FALSE.
863 118 : 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 2 : 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 116 : 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 11278 : intermeta_steps = n_step - hills_env%old_hill_step
875 11278 : force_gauss = .FALSE.
876 11278 : 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 132 : ALLOCATE (ddp(meta_env%n_colvar))
879 88 : ALLOCATE (local_last_hills(meta_env%n_colvar))
880 :
881 220 : 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 44 : dp2 = 0.0_dp
885 132 : DO i = 1, n_colvar
886 88 : ddp(i) = colvars(i)%ss0 - local_last_hills(i)
887 88 : IF (colvars(i)%periodic) THEN
888 : ! The difference of a periodic COLVAR is always within [-pi,pi]
889 0 : ddp(i) = SIGN(1.0_dp, ASIN(SIN(ddp(i))))*ACOS(COS(ddp(i)))
890 : END IF
891 132 : dp2 = dp2 + ddp(i)**2
892 : END DO
893 44 : dp2 = SQRT(dp2)
894 44 : IF (dp2 > hills_env%min_disp) THEN
895 16 : force_gauss = .TRUE.
896 : iw = cp_print_key_unit_nr(logger, meta_env%metadyn_section, &
897 16 : "PRINT%PROGRAM_RUN_INFO", extension=".metadynLog")
898 16 : IF (iw > 0) THEN
899 : WRITE (UNIT=iw, FMT="(A,F0.6,A,F0.6)") &
900 8 : " METADYN| Threshold value for COLVAR displacement exceeded: ", &
901 16 : dp2, " > ", hills_env%min_disp
902 : END IF
903 : CALL cp_print_key_finished_output(iw, logger, meta_env%metadyn_section, &
904 16 : "PRINT%PROGRAM_RUN_INFO")
905 : END IF
906 44 : DEALLOCATE (ddp)
907 44 : 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 11278 : .AND. (.NOT. meta_env%restart) .AND. (hills_env%nt_hills > 0)) THEN
913 1116 : IF (meta_env%do_multiple_walkers) multiple_walkers => meta_env%multiple_walkers
914 :
915 1116 : n_hills_start = hills_env%n_hills
916 : ! Add the hill corresponding to this location
917 1116 : IF (meta_env%well_tempered) THEN
918 : ! Well-Tempered scaling of hills height
919 : V_now_here = 0._dp
920 6 : DO ih = 1, hills_env%n_hills
921 2 : dp2 = 0._dp
922 4 : DO i = 1, n_colvar
923 2 : diff = colvars(i)%ss0 - hills_env%ss_history(i, ih)
924 2 : IF (colvars(i)%periodic) THEN
925 : ! The difference of a periodic COLVAR is always within [-pi,pi]
926 0 : diff = SIGN(1.0_dp, ASIN(SIN(diff)))*ACOS(COS(diff))
927 : END IF
928 2 : diff = (diff)/hills_env%delta_s_history(i, ih)
929 4 : dp2 = dp2 + diff**2
930 : END DO
931 2 : V_to_fes = 1.0_dp + meta_env%wttemperature*hills_env%invdt_history(ih)
932 6 : V_now_here = V_now_here + hills_env%ww_history(ih)/V_to_fes*EXP(-0.5_dp*dp2)
933 : END DO
934 4 : wtww = hills_env%ww*EXP(-V_now_here*meta_env%invdt)
935 4 : ww = wtww*(1.0_dp + meta_env%wttemperature*meta_env%invdt)
936 4 : CALL add_hill_single(hills_env, colvars, ww, hills_env%n_hills, n_colvar, meta_env%invdt)
937 : ELSE
938 1112 : 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 1116 : IF (meta_env%do_multiple_walkers) multiple_walkers%n_hills_local = multiple_walkers%n_hills_local + 1
942 :
943 1116 : hills_env%old_hill_number = hills_env%n_hills
944 1116 : hills_env%old_hill_step = n_step
945 :
946 : ! Update iteration level for printing
947 1116 : CALL get_meta_iter_level(meta_env, iter_nr)
948 1116 : 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 1116 : "PRINT%PROGRAM_RUN_INFO", extension=".metadynLog")
953 1116 : IF (iw > 0) THEN
954 558 : IF (meta_env%do_multiple_walkers) THEN
955 : WRITE (iw, '(/,1X,"METADYN|",A,I0,A,I0,A,/)') &
956 66 : ' Global/Local Hills number (', hills_env%n_hills, '/', multiple_walkers%n_hills_local, &
957 132 : ') added.'
958 : ELSE
959 492 : 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 1116 : "PRINT%PROGRAM_RUN_INFO")
964 :
965 : ! Handle Multiple Walkers
966 1116 : 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 132 : "PRINT%HILLS", middle_name="LOCAL", extension=".metadynLog")
970 132 : IF (iw > 0) THEN
971 66 : WRITE (iw, '(f12.1,30f13.5)') meta_env%time*femtoseconds, &
972 198 : (hills_env%ss_history(ih, hills_env%n_hills), ih=1, n_colvar), &
973 132 : (hills_env%delta_s_history(ih, hills_env%n_hills), ih=1, n_colvar), &
974 396 : 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 132 : "PRINT%HILLS")
978 :
979 : ! Check the communication buffer of the other walkers
980 : CALL synchronize_multiple_walkers(multiple_walkers, hills_env, colvars, &
981 132 : 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 1116 : "PRINT%HILLS", extension=".metadynLog")
988 1116 : IF (iw > 0) THEN
989 697 : DO i_hills = n_hills_start + 1, hills_env%n_hills
990 373 : WRITE (iw, '(f12.1,30f13.5)') meta_env%time*femtoseconds, &
991 1182 : (hills_env%ss_history(ih, i_hills), ih=1, n_colvar), &
992 809 : (hills_env%delta_s_history(ih, i_hills), ih=1, n_colvar), &
993 2688 : 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 1116 : "PRINT%HILLS")
998 : END IF
999 :
1000 : ! Computes forces due to the hills: history dependent term
1001 33834 : ALLOCATE (ddp(meta_env%n_colvar))
1002 22556 : ALLOCATE (diff_ss(meta_env%n_colvar))
1003 22556 : ALLOCATE (numf(meta_env%n_colvar))
1004 22556 : ALLOCATE (denf(meta_env%n_colvar))
1005 :
1006 11278 : hills_env%energy = 0.0_dp
1007 23594 : DO ih = 1, n_colvar
1008 23594 : colvars(ih)%ff_hills = 0.0_dp
1009 : END DO
1010 123556 : DO ih = 1, hills_env%n_hills
1011 112278 : slow_growth = 1.0_dp
1012 112278 : IF (hills_env%slow_growth .AND. (ih == hills_env%n_hills)) THEN
1013 3804 : slow_growth = REAL(n_step - hills_env%old_hill_step, dp)/REAL(hills_env%nt_hills, dp)
1014 : END IF
1015 112278 : dp2 = 0._dp
1016 112278 : cut_func = 1.0_dp
1017 236078 : DO i = 1, n_colvar
1018 123800 : diff_ss(i) = colvars(i)%ss0 - hills_env%ss_history(i, ih)
1019 123800 : IF (colvars(i)%periodic) THEN
1020 : ! The difference of a periodic COLVAR is always within [-pi,pi]
1021 0 : diff_ss(i) = SIGN(1.0_dp, ASIN(SIN(diff_ss(i))))*ACOS(COS(diff_ss(i)))
1022 : END IF
1023 123800 : 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 0 : ddp(i) = 0.0_dp
1029 : ELSE
1030 123800 : ddp(i) = (diff_ss(i))/hills_env%delta_s_history(i, ih)
1031 : END IF
1032 123800 : dp2 = dp2 + ddp(i)**2
1033 236078 : IF (hills_env%tail_cutoff > 0.0_dp) THEN
1034 38104 : frac = ABS(ddp(i))/hills_env%tail_cutoff
1035 38104 : numf(i) = frac**hills_env%p_exp
1036 38104 : denf(i) = frac**hills_env%q_exp
1037 38104 : 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 112278 : dfunc = hills_env%ww_history(ih)*EXP(-0.5_dp*dp2)*slow_growth
1042 112278 : IF (meta_env%well_tempered) dfunc = dfunc/(1.0_dp + meta_env%wttemperature*hills_env%invdt_history(ih))
1043 112278 : hills_env%energy = hills_env%energy + dfunc*cut_func
1044 247356 : DO i = 1, n_colvar
1045 236078 : 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 123800 : ddp(i)/hills_env%delta_s_history(i, ih)*dfunc*cut_func
1050 123800 : 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 38006 : dfunc*cut_func/ABS(diff_ss(i))
1054 : END IF
1055 : END IF
1056 : END DO
1057 : END DO
1058 11278 : DEALLOCATE (ddp)
1059 11278 : DEALLOCATE (diff_ss)
1060 11278 : DEALLOCATE (numf)
1061 11278 : DEALLOCATE (denf)
1062 :
1063 11278 : CALL cp_rm_iter_level(logger%iter_info, "METADYNAMICS")
1064 :
1065 11278 : CALL timestop(handle)
1066 :
1067 22556 : END SUBROUTINE hills
1068 :
1069 : END MODULE metadynamics
|
