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