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