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_utils 13 USE cp_files, ONLY: close_file,& 14 open_file 15 USE cp_log_handling, ONLY: cp_get_default_logger,& 16 cp_logger_type,& 17 cp_to_string 18 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 19 cp_print_key_unit_nr 20 USE cp_para_types, ONLY: cp_para_env_type 21 USE cp_subsys_types, ONLY: cp_subsys_type 22 USE force_env_types, ONLY: force_env_get,& 23 force_env_type 24 USE input_constants, ONLY: do_fe_meta,& 25 do_wall_gaussian,& 26 do_wall_m,& 27 do_wall_none,& 28 do_wall_p,& 29 do_wall_quadratic,& 30 do_wall_quartic,& 31 do_wall_reflective 32 USE input_cp2k_free_energy, ONLY: create_metavar_section 33 USE input_enumeration_types, ONLY: enum_i2c,& 34 enumeration_type 35 USE input_keyword_types, ONLY: keyword_get,& 36 keyword_type 37 USE input_section_types, ONLY: section_get_keyword,& 38 section_get_subsection,& 39 section_release,& 40 section_type,& 41 section_vals_get,& 42 section_vals_get_subs_vals,& 43 section_vals_type,& 44 section_vals_val_get 45 USE kinds, ONLY: default_path_length,& 46 dp 47 USE machine, ONLY: m_mov 48 USE message_passing, ONLY: mp_bcast 49 USE metadynamics_types, ONLY: hills_env_type,& 50 meta_env_type,& 51 metadyn_create,& 52 metavar_type,& 53 multiple_walkers_type 54 USE physcon, ONLY: kelvin 55#include "./base/base_uses.f90" 56 57 IMPLICIT NONE 58 PRIVATE 59 60 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'metadynamics_utils' 61 62 PUBLIC :: metadyn_read, & 63 synchronize_multiple_walkers, & 64 add_hill_single, & 65 restart_hills, & 66 get_meta_iter_level, & 67 meta_walls 68 69CONTAINS 70 71! ************************************************************************************************** 72!> \brief reads metadynamics section 73!> \param meta_env ... 74!> \param force_env ... 75!> \param root_section ... 76!> \param para_env ... 77!> \param fe_section ... 78!> \par History 79!> 04.2004 created 80!> \author Teodoro Laino [tlaino] - University of Zurich. 11.2007 81! ************************************************************************************************** 82 SUBROUTINE metadyn_read(meta_env, force_env, root_section, para_env, fe_section) 83 TYPE(meta_env_type), POINTER :: meta_env 84 TYPE(force_env_type), POINTER :: force_env 85 TYPE(section_vals_type), POINTER :: root_section 86 TYPE(cp_para_env_type), POINTER :: para_env 87 TYPE(section_vals_type), OPTIONAL, POINTER :: fe_section 88 89 CHARACTER(len=*), PARAMETER :: routineN = 'metadyn_read', routineP = moduleN//':'//routineN 90 91 CHARACTER(LEN=default_path_length) :: walkers_file_name 92 INTEGER :: handle, i, id_method, n_colvar, n_rep, & 93 number_allocated_colvars 94 INTEGER, DIMENSION(:), POINTER :: walkers_status 95 LOGICAL :: check, explicit 96 REAL(kind=dp) :: dt 97 TYPE(cp_subsys_type), POINTER :: subsys 98 TYPE(section_vals_type), POINTER :: md_section, metadyn_section, & 99 metavar_section, walkers_section 100 101 NULLIFY (subsys) 102 CALL timeset(routineN, handle) 103 104 CALL section_vals_get(fe_section, explicit=explicit) 105 IF (explicit) THEN 106 number_allocated_colvars = 0 107 CALL force_env_get(force_env, subsys=subsys) 108 IF (ASSOCIATED(subsys%colvar_p)) THEN 109 number_allocated_colvars = SIZE(subsys%colvar_p) 110 END IF 111 CALL section_vals_val_get(fe_section, "METHOD", i_val=id_method) 112 IF (id_method /= do_fe_meta) THEN 113 CALL timestop(handle) 114 RETURN 115 ENDIF 116 metadyn_section => section_vals_get_subs_vals(fe_section, "METADYN") 117 CPASSERT(.NOT. ASSOCIATED(meta_env)) 118 119 md_section => section_vals_get_subs_vals(root_section, "MOTION%MD") 120 CALL section_vals_val_get(md_section, "TIMESTEP", r_val=dt) 121 122 metavar_section => section_vals_get_subs_vals(metadyn_section, "METAVAR") 123 CALL section_vals_get(metavar_section, n_repetition=n_colvar) 124 CALL metadyn_create(meta_env, n_colvar=n_colvar, & 125 dt=dt, para_env=para_env, metadyn_section=metadyn_section) 126 127 !Check if using plumed. If so, only get the file name and read nothing else 128 CALL section_vals_val_get(metadyn_section, "USE_PLUMED", l_val=meta_env%use_plumed) 129 IF (meta_env%use_plumed .EQV. .TRUE.) THEN 130 CALL section_vals_val_get(metadyn_section, "PLUMED_INPUT_FILE", c_val=meta_env%plumed_input_file) 131 meta_env%plumed_input_file = TRIM(meta_env%plumed_input_file)//CHAR(0) 132 meta_env%langevin = .FALSE. 133 CALL timestop(handle) 134 RETURN 135 END IF 136 137 CALL section_vals_val_get(metadyn_section, "DO_HILLS", l_val=meta_env%do_hills) 138 CALL section_vals_val_get(metadyn_section, "LAGRANGE", l_val=meta_env%extended_lagrange) 139 CALL section_vals_val_get(metadyn_section, "TAMCSteps", i_val=meta_env%TAMCSteps) 140 IF (meta_env%TAMCSteps < 0) THEN 141 CPABORT("TAMCSteps must be positive!") 142 ENDIF 143 CALL section_vals_val_get(metadyn_section, "Timestep", r_val=meta_env%zdt) 144 IF (meta_env%zdt <= 0.0_dp) THEN 145 CPABORT("Timestep must be positive!") 146 ENDIF 147 CALL section_vals_val_get(metadyn_section, "WW", r_val=meta_env%hills_env%ww) 148 CALL section_vals_val_get(metadyn_section, "NT_HILLS", i_val=meta_env%hills_env%nt_hills) 149 CALL section_vals_val_get(metadyn_section, "MIN_NT_HILLS", i_val=meta_env%hills_env%min_nt_hills) 150 IF (meta_env%hills_env%nt_hills <= 0) THEN 151 meta_env%hills_env%min_nt_hills = meta_env%hills_env%nt_hills 152 CALL cp_warn(__LOCATION__, & 153 "NT_HILLS has a value <= 0; "// & 154 "Setting MIN_NT_HILLS to the same value! "// & 155 "Overriding input specification!") 156 END IF 157 check = meta_env%hills_env%nt_hills >= meta_env%hills_env%min_nt_hills 158 IF (.NOT. check) & 159 CALL cp_abort(__LOCATION__, "MIN_NT_HILLS must have a value smaller or equal to NT_HILLS! "// & 160 "Cross check with the input reference!") 161 !RG Adaptive hills 162 CALL section_vals_val_get(metadyn_section, "MIN_DISP", r_val=meta_env%hills_env%min_disp) 163 CALL section_vals_val_get(metadyn_section, "OLD_HILL_NUMBER", i_val=meta_env%hills_env%old_hill_number) 164 CALL section_vals_val_get(metadyn_section, "OLD_HILL_STEP", i_val=meta_env%hills_env%old_hill_step) 165 166 !Hills tail damping 167 CALL section_vals_val_get(metadyn_section, "HILL_TAIL_CUTOFF", r_val=meta_env%hills_env%tail_cutoff) 168 CALL section_vals_val_get(metadyn_section, "P_EXPONENT", i_val=meta_env%hills_env%p_exp) 169 CALL section_vals_val_get(metadyn_section, "Q_EXPONENT", i_val=meta_env%hills_env%q_exp) 170 171 CALL section_vals_val_get(metadyn_section, "SLOW_GROWTH", l_val=meta_env%hills_env%slow_growth) 172 173 !RG Adaptive hills 174 CALL section_vals_val_get(metadyn_section, "STEP_START_VAL", i_val=meta_env%n_steps) 175 CPASSERT(meta_env%n_steps >= 0) 176 CALL section_vals_val_get(metadyn_section, "NHILLS_START_VAL", & 177 i_val=meta_env%hills_env%n_hills) 178 CALL section_vals_val_get(metadyn_section, "TEMPERATURE", r_val=meta_env%temp_wanted) 179 CALL section_vals_val_get(metadyn_section, "LANGEVIN", l_val=meta_env%langevin) 180 CALL section_vals_val_get(metadyn_section, "TEMP_TOL", explicit=meta_env%tempcontrol, & 181 r_val=meta_env%toll_temp) 182 CALL section_vals_val_get(metadyn_section, "WELL_TEMPERED", l_val=meta_env%well_tempered) 183 CALL section_vals_val_get(metadyn_section, "DELTA_T", explicit=meta_env%hills_env%wtcontrol, & 184 r_val=meta_env%delta_t) 185 CALL section_vals_val_get(metadyn_section, "WTGAMMA", explicit=check, & 186 r_val=meta_env%wtgamma) 187 IF (meta_env%well_tempered) THEN 188 meta_env%hills_env%wtcontrol = meta_env%hills_env%wtcontrol .OR. check 189 check = meta_env%hills_env%wtcontrol 190 IF (.NOT. check) & 191 CALL cp_abort(__LOCATION__, "When using Well-Tempered metadynamics, "// & 192 "DELTA_T (or WTGAMMA) should be explicitly specified.") 193 IF (meta_env%extended_lagrange) & 194 CALL cp_abort(__LOCATION__, & 195 "Well-Tempered metadynamics not possible with extended-lagrangian formulation.") 196 IF (meta_env%hills_env%min_disp > 0.0_dp) & 197 CALL cp_abort(__LOCATION__, & 198 "Well-Tempered metadynamics not possible with Adaptive hills.") 199 ENDIF 200 201 CALL section_vals_val_get(metadyn_section, "COLVAR_AVG_TEMPERATURE_RESTART", & 202 r_val=meta_env%avg_temp) 203 ! Parsing Metavar Section 204 DO i = 1, n_colvar 205 CALL metavar_read(meta_env%metavar(i), meta_env%extended_lagrange, & 206 meta_env%langevin, i, metavar_section) 207 check = (meta_env%metavar(i)%icolvar <= number_allocated_colvars) 208 IF (.NOT. check) & 209 CALL cp_abort(__LOCATION__, & 210 "An error occured in the specification of COLVAR for METAVAR. "// & 211 "Specified COLVAR #("//TRIM(ADJUSTL(cp_to_string(meta_env%metavar(i)%icolvar)))//") "// & 212 "is larger than the maximum number of COLVARS defined in the SUBSYS ("// & 213 TRIM(ADJUSTL(cp_to_string(number_allocated_colvars)))//") !") 214 END DO 215 216 ! Parsing the Multiple Walkers Info 217 IF (meta_env%do_multiple_walkers) THEN 218 NULLIFY (walkers_status) 219 walkers_section => section_vals_get_subs_vals(metadyn_section, "MULTIPLE_WALKERS") 220 221 ! General setup for walkers 222 CALL section_vals_val_get(walkers_section, "WALKER_ID", & 223 i_val=meta_env%multiple_walkers%walker_id) 224 CALL section_vals_val_get(walkers_section, "NUMBER_OF_WALKERS", & 225 i_val=meta_env%multiple_walkers%walkers_tot_nr) 226 CALL section_vals_val_get(walkers_section, "WALKER_COMM_FREQUENCY", & 227 i_val=meta_env%multiple_walkers%walkers_freq_comm) 228 229 ! Handle status and file names 230 ALLOCATE (meta_env%multiple_walkers%walkers_status(meta_env%multiple_walkers%walkers_tot_nr)) 231 ALLOCATE (meta_env%multiple_walkers%walkers_file_name(meta_env%multiple_walkers%walkers_tot_nr)) 232 CALL section_vals_val_get(walkers_section, "WALKERS_STATUS", explicit=explicit) 233 IF (explicit) THEN 234 CALL section_vals_val_get(walkers_section, "WALKERS_STATUS", i_vals=walkers_status) 235 check = (SIZE(walkers_status) == meta_env%multiple_walkers%walkers_tot_nr) 236 IF (.NOT. check) & 237 CALL cp_abort(__LOCATION__, & 238 "Number of Walkers specified in the input does not match with the "// & 239 "size of the WALKERS_STATUS. Please check your input and in case "// & 240 "this is a restart run consider the possibility to switch off the "// & 241 "RESTART_WALKERS in the EXT_RESTART section! ") 242 meta_env%multiple_walkers%walkers_status = walkers_status 243 ELSE 244 meta_env%multiple_walkers%walkers_status = 0 245 END IF 246 meta_env%multiple_walkers%n_hills_local = & 247 meta_env%multiple_walkers%walkers_status(meta_env%multiple_walkers%walker_id) 248 249 CALL section_vals_val_get(walkers_section, "WALKERS_FILE_NAME%_DEFAULT_KEYWORD_", & 250 n_rep_val=n_rep) 251 check = (n_rep == meta_env%multiple_walkers%walkers_tot_nr) 252 IF (.NOT. check) & 253 CALL cp_abort(__LOCATION__, & 254 "Number of Walkers specified in the input does not match with the "// & 255 "number of Walkers File names provided. Please check your input and in case "// & 256 "this is a restart run consider the possibility to switch off the "// & 257 "RESTART_WALKERS in the EXT_RESTART section! ") 258 DO i = 1, n_rep 259 CALL section_vals_val_get(walkers_section, "WALKERS_FILE_NAME%_DEFAULT_KEYWORD_", & 260 i_rep_val=i, c_val=walkers_file_name) 261 meta_env%multiple_walkers%walkers_file_name(i) = walkers_file_name 262 END DO 263 END IF 264 265 ! Print Metadynamics Info 266 CALL print_metadyn_info(meta_env, n_colvar, metadyn_section) 267 END IF 268 269 CALL timestop(handle) 270 271 END SUBROUTINE metadyn_read 272 273! ************************************************************************************************** 274!> \brief prints information on the metadynamics run 275!> \param meta_env ... 276!> \param n_colvar ... 277!> \param metadyn_section ... 278!> \author Teodoro Laino [tlaino] - University of Zurich. 10.2008 279! ************************************************************************************************** 280 SUBROUTINE print_metadyn_info(meta_env, n_colvar, metadyn_section) 281 TYPE(meta_env_type), POINTER :: meta_env 282 INTEGER, INTENT(IN) :: n_colvar 283 TYPE(section_vals_type), POINTER :: metadyn_section 284 285 CHARACTER(len=*), PARAMETER :: routineN = 'print_metadyn_info', & 286 routineP = moduleN//':'//routineN 287 288 CHARACTER(LEN=10) :: my_id, my_tag 289 INTEGER :: handle, i, iw, j 290 TYPE(cp_logger_type), POINTER :: logger 291 TYPE(enumeration_type), POINTER :: enum 292 TYPE(keyword_type), POINTER :: keyword 293 TYPE(section_type), POINTER :: section, wall_section, work_section 294 295 CALL timeset(routineN, handle) 296 297 logger => cp_get_default_logger() 298 iw = cp_print_key_unit_nr(logger, metadyn_section, & 299 "PRINT%PROGRAM_RUN_INFO", extension=".metadynLog") 300 NULLIFY (section, enum, keyword) 301 CALL create_metavar_section(section) 302 wall_section => section_get_subsection(section, "WALL") 303 IF (iw > 0) THEN 304 WRITE (iw, '( /A )') ' METADYN| Meta Dynamics Protocol ' 305 WRITE (iw, '( A,T71,I10)') ' METADYN| Number of interval time steps to spawn hills', & 306 meta_env%hills_env%nt_hills 307 WRITE (iw, '( A,T71,I10)') ' METADYN| Number of previously spawned hills', & 308 meta_env%hills_env%n_hills 309 IF (meta_env%extended_lagrange) THEN 310 WRITE (iw, '( A )') ' METADYN| Extended Lagrangian Scheme ' 311 IF (meta_env%tempcontrol) WRITE (iw, '( A,T71,F10.2)') & 312 ' METADYN| Collective Variables Temperature control', meta_env%toll_temp 313 IF (meta_env%langevin) THEN 314 WRITE (iw, '(A,T71)') ' METADYN| Langevin Thermostat in use for COLVAR ' 315 WRITE (iw, '(A,T71,F10.4)') ' METADYN| Langevin Thermostat. Target Temperature = ', & 316 meta_env%temp_wanted*kelvin 317 ENDIF 318 WRITE (iw, '(A,T71,F10.4)') ' METADYN| COLVARS restarted average temperature ', & 319 meta_env%avg_temp 320 END IF 321 IF (meta_env%do_hills) THEN 322 WRITE (iw, '( A )') ' METADYN| Spawning the Hills ' 323 WRITE (iw, '( A,T71,F10.3)') ' METADYN| Height of the Spawned Gaussian', meta_env%hills_env%ww 324 !RG Adaptive hills 325 IF (meta_env%hills_env%min_disp .GT. 0.0_dp) THEN 326 WRITE (iw, '(A)') ' METADYN| Adapative meta time step is activated' 327 WRITE (iw, '(A,T71,F10.4)') ' METADYN| Minimum displacement for next hill', & 328 meta_env%hills_env%min_disp 329 END IF 330 !RG Adaptive hills 331 END IF 332 333 IF (meta_env%well_tempered) THEN 334 WRITE (iw, '( A )') ' METADYN| Well-Tempered metadynamics ' 335 IF (meta_env%delta_t > EPSILON(1._dp)) THEN 336 WRITE (iw, '( A,T71,F10.3)') ' METADYN| Temperature parameter (Delta T) [K]', meta_env%delta_t*kelvin 337 ELSE 338 WRITE (iw, '( A,T71,F10.3)') ' METADYN| Temperature parameter (gamma)', meta_env%wtgamma 339 END IF 340 END IF 341 342 IF (meta_env%do_multiple_walkers) THEN 343 WRITE (iw, '( A,T71,A10)') ' METADYN| Multiple Walkers', ' ENABLED' 344 WRITE (iw, '( A,T71,I10)') ' METADYN| Number of Multiple Walkers', & 345 meta_env%multiple_walkers%walkers_tot_nr 346 WRITE (iw, '( A,T71,I10)') ' METADYN| Local Walker ID', & 347 meta_env%multiple_walkers%walker_id 348 WRITE (iw, '( A,T71,I10)') ' METADYN| Walker Communication Frequency', & 349 meta_env%multiple_walkers%walkers_freq_comm 350 DO i = 1, meta_env%multiple_walkers%walkers_tot_nr 351 my_tag = "" 352 IF (i == meta_env%multiple_walkers%walker_id) my_tag = " ( Local )" 353 my_id = '( '//TRIM(ADJUSTL(cp_to_string(i)))//' )' 354 WRITE (iw, '(/,A,T71,A10)') ' WALKERS| Walker ID'//TRIM(my_tag), ADJUSTR(my_id) 355 WRITE (iw, '( A,T71,I10)') ' WALKERS| Number of Hills communicated', & 356 meta_env%multiple_walkers%walkers_status(i) 357 WRITE (iw, '( A,T24,A57)') ' WALKERS| Base Filename', & 358 ADJUSTR(meta_env%multiple_walkers%walkers_file_name(i) (1:57)) 359 END DO 360 WRITE (iw, '(/)') 361 END IF 362 363 WRITE (iw, '( A,T71,I10)') ' METADYN| Number of collective variables', meta_env%n_colvar 364 DO i = 1, n_colvar 365 WRITE (iw, '( A )') ' '//'----------------------------------------------------------------------' 366 WRITE (iw, '( A,T71,I10)') ' METAVARS| Collective Variable Number', meta_env%metavar(i)%icolvar 367 IF (meta_env%extended_lagrange) THEN 368 WRITE (iw, '( A,T71,F10.6)') ' METAVARS| Lambda Parameter', meta_env%metavar(i)%lambda 369 WRITE (iw, '( A,T66,F15.6)') ' METAVARS| Collective Variable Mass', meta_env%metavar(i)%mass 370 END IF 371 WRITE (iw, '( A,T71,F10.6)') ' METAVARS| Scaling factor', meta_env%metavar(i)%delta_s 372 IF (meta_env%langevin) WRITE (iw, '( A,T71,F10.6)') ' METAVARS| Friction for Langevin Thermostat', & 373 meta_env%metavar(i)%gamma 374 IF (meta_env%metavar(i)%do_wall) THEN 375 WRITE (iw, '( A,T71,I10)') ' METAVARS| Number of Walls present', SIZE(meta_env%metavar(i)%walls) 376 DO j = 1, SIZE(meta_env%metavar(i)%walls) 377 keyword => section_get_keyword(wall_section, "TYPE") 378 CALL keyword_get(keyword, enum=enum) 379 WRITE (iw, '(/,A,5X,I10,T50,A,T70,A11)') ' METAVARS| Wall Number:', j, 'Type of Wall:', & 380 ADJUSTR(TRIM(enum_i2c(enum, meta_env%metavar(i)%walls(j)%id_type))) 381 ! Type of wall IO 382 SELECT CASE (meta_env%metavar(i)%walls(j)%id_type) 383 CASE (do_wall_none) 384 ! Do Nothing 385 CYCLE 386 CASE (do_wall_reflective) 387 work_section => section_get_subsection(wall_section, "REFLECTIVE") 388 keyword => section_get_keyword(work_section, "DIRECTION") 389 CALL keyword_get(keyword, enum=enum) 390 WRITE (iw, '(A,T70,A11)') ' METAVARS| Wall direction', & 391 ADJUSTR(TRIM(enum_i2c(enum, meta_env%metavar(i)%walls(j)%id_direction))) 392 CASE (do_wall_quadratic) 393 work_section => section_get_subsection(wall_section, "QUADRATIC") 394 keyword => section_get_keyword(work_section, "DIRECTION") 395 CALL keyword_get(keyword, enum=enum) 396 WRITE (iw, '(A,T70,A11)') ' METAVARS| Wall direction', & 397 ADJUSTR(TRIM(enum_i2c(enum, meta_env%metavar(i)%walls(j)%id_direction))) 398 WRITE (iw, '(A,T70,F11.6)') ' METAVARS| Constant K of the quadratic potential', & 399 meta_env%metavar(i)%walls(j)%k_quadratic 400 CASE (do_wall_gaussian) 401 WRITE (iw, '(A,T70,F11.6)') ' METAVARS| Height of the Wall Gaussian', & 402 meta_env%metavar(i)%walls(j)%ww_gauss 403 WRITE (iw, '(A,T70,F11.6)') ' METAVARS| Scale of the Wall Gaussian', & 404 meta_env%metavar(i)%walls(j)%sigma_gauss 405 END SELECT 406 WRITE (iw, '(A,T70,F11.6)') ' METAVARS| Wall location', & 407 meta_env%metavar(i)%walls(j)%pos 408 END DO 409 END IF 410 WRITE (iw, '( A )') ' '//'----------------------------------------------------------------------' 411 ENDDO 412 END IF 413 CALL section_release(section) 414 CALL cp_print_key_finished_output(iw, logger, metadyn_section, "PRINT%PROGRAM_RUN_INFO") 415 416 CALL timestop(handle) 417 418 END SUBROUTINE print_metadyn_info 419 420! ************************************************************************************************** 421!> \brief reads metavar section 422!> \param metavar ... 423!> \param extended_lagrange ... 424!> \param langevin ... 425!> \param icol ... 426!> \param metavar_section ... 427!> \par History 428!> 04.2004 created 429!> \author alessandro laio and fawzi mohamed 430!> Teodoro Laino [tlaino] - University of Zurich. 11.2007 431! ************************************************************************************************** 432 SUBROUTINE metavar_read(metavar, extended_lagrange, langevin, icol, metavar_section) 433 TYPE(metavar_type), INTENT(INOUT) :: metavar 434 LOGICAL, INTENT(IN) :: extended_lagrange, langevin 435 INTEGER, INTENT(IN) :: icol 436 TYPE(section_vals_type), OPTIONAL, POINTER :: metavar_section 437 438 CHARACTER(len=*), PARAMETER :: routineN = 'metavar_read', routineP = moduleN//':'//routineN 439 440 INTEGER :: handle, i, n_walls 441 TYPE(section_vals_type), POINTER :: wall_section, work_section 442 443 CALL timeset(routineN, handle) 444 445 CALL section_vals_val_get(metavar_section, "COLVAR", i_rep_section=icol, i_val=metavar%icolvar) 446 CALL section_vals_val_get(metavar_section, "SCALE", i_rep_section=icol, r_val=metavar%delta_s) 447 ! Walls 448 wall_section => section_vals_get_subs_vals(metavar_section, "WALL", i_rep_section=icol) 449 CALL section_vals_get(wall_section, n_repetition=n_walls) 450 IF (n_walls /= 0) THEN 451 metavar%do_wall = .TRUE. 452 ALLOCATE (metavar%walls(n_walls)) 453 DO i = 1, n_walls 454 CALL section_vals_val_get(wall_section, "TYPE", i_rep_section=i, i_val=metavar%walls(i)%id_type) 455 CALL section_vals_val_get(wall_section, "POSITION", i_rep_section=i, r_val=metavar%walls(i)%pos) 456 SELECT CASE (metavar%walls(i)%id_type) 457 CASE (do_wall_none) 458 ! Just cycle.. 459 CYCLE 460 CASE (do_wall_reflective) 461 work_section => section_vals_get_subs_vals(wall_section, "REFLECTIVE", i_rep_section=i) 462 CALL section_vals_val_get(work_section, "DIRECTION", i_val=metavar%walls(i)%id_direction) 463 CASE (do_wall_quadratic) 464 work_section => section_vals_get_subs_vals(wall_section, "QUADRATIC", i_rep_section=i) 465 CALL section_vals_val_get(work_section, "DIRECTION", i_val=metavar%walls(i)%id_direction) 466 CALL section_vals_val_get(work_section, "K", r_val=metavar%walls(i)%k_quadratic) 467 CASE (do_wall_quartic) 468 work_section => section_vals_get_subs_vals(wall_section, "QUARTIC", i_rep_section=i) 469 CALL section_vals_val_get(work_section, "DIRECTION", i_val=metavar%walls(i)%id_direction) 470 CALL section_vals_val_get(work_section, "K", r_val=metavar%walls(i)%k_quartic) 471 SELECT CASE (metavar%walls(i)%id_direction) 472 CASE (do_wall_m) 473 metavar%walls(i)%pos0 = metavar%walls(i)%pos + (0.05_dp/metavar%walls(i)%k_quartic**(0.25_dp)) 474 CASE (do_wall_p) 475 metavar%walls(i)%pos0 = metavar%walls(i)%pos - (0.05_dp/metavar%walls(i)%k_quartic**(0.25_dp)) 476 END SELECT 477 CASE (do_wall_gaussian) 478 work_section => section_vals_get_subs_vals(wall_section, "GAUSSIAN", i_rep_section=i) 479 CALL section_vals_val_get(work_section, "WW", r_val=metavar%walls(i)%ww_gauss) 480 CALL section_vals_val_get(work_section, "SIGMA", r_val=metavar%walls(i)%sigma_gauss) 481 END SELECT 482 END DO 483 END IF 484 ! Setup few more parameters for extended lagrangian 485 IF (extended_lagrange) THEN 486 CALL section_vals_val_get(metavar_section, "MASS", i_rep_section=icol, r_val=metavar%mass) 487 CALL section_vals_val_get(metavar_section, "LAMBDA", i_rep_section=icol, r_val=metavar%lambda) 488 IF (langevin) THEN 489 CALL section_vals_val_get(metavar_section, "GAMMA", i_rep_section=icol, r_val=metavar%gamma) 490 END IF 491 ENDIF 492 493 CALL timestop(handle) 494 495 END SUBROUTINE metavar_read 496 497! ************************************************************************************************** 498!> \brief Synchronize with the rest of the walkers 499!> \param multiple_walkers ... 500!> \param hills_env ... 501!> \param colvars ... 502!> \param n_colvar ... 503!> \param metadyn_section ... 504!> \author Teodoro Laino [tlaino] - University of Zurich - 10.2008 505! ************************************************************************************************** 506 SUBROUTINE synchronize_multiple_walkers(multiple_walkers, hills_env, colvars, & 507 n_colvar, metadyn_section) 508 TYPE(multiple_walkers_type), POINTER :: multiple_walkers 509 TYPE(hills_env_type), POINTER :: hills_env 510 TYPE(metavar_type), DIMENSION(:), POINTER :: colvars 511 INTEGER, INTENT(IN) :: n_colvar 512 TYPE(section_vals_type), POINTER :: metadyn_section 513 514 CHARACTER(len=*), PARAMETER :: routineN = 'synchronize_multiple_walkers', & 515 routineP = moduleN//':'//routineN 516 517 CHARACTER(LEN=default_path_length) :: filename, tmpname 518 INTEGER :: delta_hills, handle, i, i_hills, ih, iw, & 519 unit_nr 520 LOGICAL :: exist 521 REAL(KIND=dp) :: invdt, ww 522 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: delta_s_save, ss0_save 523 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: delta_s_ss0_buf 524 TYPE(cp_logger_type), POINTER :: logger 525 TYPE(cp_para_env_type), POINTER :: para_env 526 527 CALL timeset(routineN, handle) 528 529 logger => cp_get_default_logger() 530 para_env => logger%para_env 531 532 ! Locally dump information on file.. 533 IF (para_env%ionode) THEN 534 ! Generate file name for the specific Hill 535 i = multiple_walkers%walker_id 536 filename = TRIM(multiple_walkers%walkers_file_name(i))//"_"// & 537 TRIM(ADJUSTL(cp_to_string(multiple_walkers%n_hills_local))) 538 tmpname = TRIM(filename)//".tmp" 539 CALL open_file(file_name=tmpname, file_status="UNKNOWN", & 540 file_form="FORMATTED", file_action="WRITE", & 541 file_position="APPEND", unit_number=unit_nr) 542 WRITE (unit_nr, *) hills_env%ww_history(hills_env%n_hills) 543 DO ih = 1, n_colvar 544 WRITE (unit_nr, *) hills_env%ss_history(ih, hills_env%n_hills) 545 WRITE (unit_nr, *) hills_env%delta_s_history(ih, hills_env%n_hills) 546 END DO 547 IF (hills_env%wtcontrol) WRITE (unit_nr, *) hills_env%invdt_history(hills_env%n_hills) 548 CALL close_file(unit_nr) 549 CALL m_mov(tmpname, filename) 550 END IF 551 552 IF (MODULO(multiple_walkers%n_hills_local, multiple_walkers%walkers_freq_comm) == 0) THEN 553 ! Store colvars information 554 ALLOCATE (ss0_save(n_colvar)) 555 ALLOCATE (delta_s_save(n_colvar)) 556 ALLOCATE (delta_s_ss0_buf(2, 0:n_colvar)) 557 delta_s_ss0_buf = 0 558 DO i = 1, n_colvar 559 ss0_save(i) = colvars(i)%ss0 560 delta_s_save(i) = colvars(i)%delta_s 561 END DO 562 563 ! Watch for other walkers's file and update 564 DO i = 1, multiple_walkers%walkers_tot_nr 565 IF (i == multiple_walkers%walker_id) THEN 566 ! Update local counter 567 multiple_walkers%walkers_status(i) = multiple_walkers%n_hills_local 568 CYCLE 569 END IF 570 571 i_hills = multiple_walkers%walkers_status(i) + 1 572 filename = TRIM(multiple_walkers%walkers_file_name(i))//"_"// & 573 TRIM(ADJUSTL(cp_to_string(i_hills))) 574 575 IF (para_env%ionode) THEN 576 INQUIRE (FILE=TRIM(filename), EXIST=exist) 577 ENDIF 578 CALL mp_bcast(exist, para_env%source, para_env%group) 579 DO WHILE (exist) 580 ! Read information from the walker's file 581 ! We shouldn't care too much about the concurrency of these I/O instructions.. 582 ! In case, they can be fixed in the future.. 583 IF (para_env%ionode) THEN 584 CALL open_file(file_name=filename, file_status="OLD", & 585 file_form="FORMATTED", file_action="READ", & 586 file_position="REWIND", unit_number=unit_nr) 587 READ (unit_nr, *) delta_s_ss0_buf(1, 0) 588 DO ih = 1, n_colvar 589 READ (unit_nr, *) delta_s_ss0_buf(1, ih) 590 READ (unit_nr, *) delta_s_ss0_buf(2, ih) 591 END DO 592 IF (hills_env%wtcontrol) READ (unit_nr, *) delta_s_ss0_buf(2, 0) 593 CALL close_file(unit_nr) 594 ENDIF 595 CALL mp_bcast(delta_s_ss0_buf, para_env%source, para_env%group) 596 ww = delta_s_ss0_buf(1, 0) 597 IF (hills_env%wtcontrol) invdt = delta_s_ss0_buf(2, 0) 598 DO ih = 1, n_colvar 599 colvars(ih)%ss0 = delta_s_ss0_buf(1, ih) 600 colvars(ih)%delta_s = delta_s_ss0_buf(2, ih) 601 ENDDO 602 603 ! Add this hill to the history dependent terms 604 IF (hills_env%wtcontrol) THEN 605 CALL add_hill_single(hills_env, colvars, ww, hills_env%n_hills, n_colvar, invdt=invdt) 606 ELSE 607 CALL add_hill_single(hills_env, colvars, ww, hills_env%n_hills, n_colvar) 608 END IF 609 610 i_hills = i_hills + 1 611 filename = TRIM(multiple_walkers%walkers_file_name(i))//"_"// & 612 TRIM(ADJUSTL(cp_to_string(i_hills))) 613 IF (para_env%ionode) THEN 614 INQUIRE (FILE=TRIM(filename), EXIST=exist) 615 ENDIF 616 CALL mp_bcast(exist, para_env%source, para_env%group) 617 END DO 618 619 delta_hills = i_hills - 1 - multiple_walkers%walkers_status(i) 620 multiple_walkers%walkers_status(i) = i_hills - 1 621 iw = cp_print_key_unit_nr(logger, metadyn_section, "PRINT%PROGRAM_RUN_INFO", & 622 extension=".metadynLog") 623 IF (iw > 0) THEN 624 WRITE (iw, '(T2,A,I0,A,I0,A,I0,A)') 'WALKERS| Walker #', i, '. Reading [', delta_hills, & 625 '] Hills. Total number of Hills acquired [', multiple_walkers%walkers_status(i), ']' 626 END IF 627 CALL cp_print_key_finished_output(iw, logger, metadyn_section, & 628 "PRINT%PROGRAM_RUN_INFO") 629 END DO 630 631 ! Restore colvars information 632 DO i = 1, n_colvar 633 colvars(i)%ss0 = ss0_save(i) 634 colvars(i)%delta_s = delta_s_save(i) 635 END DO 636 DEALLOCATE (ss0_save) 637 DEALLOCATE (delta_s_save) 638 END IF 639 640 CALL timestop(handle) 641 642 END SUBROUTINE synchronize_multiple_walkers 643 644! ************************************************************************************************** 645!> \brief Add a single Hill 646!> \param hills_env ... 647!> \param colvars ... 648!> \param ww ... 649!> \param n_hills ... 650!> \param n_colvar ... 651!> \param invdt ... 652!> \author Teodoro Laino [tlaino] - University of Zurich - 10.2008 653! ************************************************************************************************** 654 SUBROUTINE add_hill_single(hills_env, colvars, ww, n_hills, n_colvar, invdt) 655 TYPE(hills_env_type), POINTER :: hills_env 656 TYPE(metavar_type), DIMENSION(:), POINTER :: colvars 657 REAL(KIND=dp), INTENT(IN) :: ww 658 INTEGER, INTENT(INOUT) :: n_hills 659 INTEGER, INTENT(IN) :: n_colvar 660 REAL(KIND=dp), INTENT(IN), OPTIONAL :: invdt 661 662 CHARACTER(len=*), PARAMETER :: routineN = 'add_hill_single', & 663 routineP = moduleN//':'//routineN 664 665 INTEGER :: handle, i 666 LOGICAL :: wtcontrol 667 REAL(KIND=dp), DIMENSION(:), POINTER :: tnp 668 REAL(KIND=dp), DIMENSION(:, :), POINTER :: tmp 669 670 CALL timeset(routineN, handle) 671 672 wtcontrol = PRESENT(invdt) 673 NULLIFY (tmp, tnp) 674 IF (SIZE(hills_env%ss_history, 2) < n_hills + 1) THEN 675 ALLOCATE (tmp(n_colvar, n_hills + 100)) 676 tmp(:, :n_hills) = hills_env%ss_history 677 tmp(:, n_hills + 1:) = 0.0_dp 678 DEALLOCATE (hills_env%ss_history) 679 hills_env%ss_history => tmp 680 NULLIFY (tmp) 681 ENDIF 682 IF (SIZE(hills_env%delta_s_history, 2) < n_hills + 1) THEN 683 ALLOCATE (tmp(n_colvar, n_hills + 100)) 684 tmp(:, :n_hills) = hills_env%delta_s_history 685 tmp(:, n_hills + 1:) = 0.0_dp 686 DEALLOCATE (hills_env%delta_s_history) 687 hills_env%delta_s_history => tmp 688 NULLIFY (tmp) 689 ENDIF 690 IF (SIZE(hills_env%ww_history) < n_hills + 1) THEN 691 ALLOCATE (tnp(n_hills + 100)) 692 tnp(1:n_hills) = hills_env%ww_history 693 tnp(n_hills + 1:) = 0.0_dp 694 DEALLOCATE (hills_env%ww_history) 695 hills_env%ww_history => tnp 696 NULLIFY (tnp) 697 ENDIF 698 IF (wtcontrol) THEN 699 IF (SIZE(hills_env%invdt_history) < n_hills + 1) THEN 700 ALLOCATE (tnp(n_hills + 100)) 701 tnp(1:n_hills) = hills_env%invdt_history 702 tnp(n_hills + 1:) = 0.0_dp 703 DEALLOCATE (hills_env%invdt_history) 704 hills_env%invdt_history => tnp 705 NULLIFY (tnp) 706 ENDIF 707 ENDIF 708 n_hills = n_hills + 1 709 ! Now add the hill 710 DO i = 1, n_colvar 711 hills_env%ss_history(i, n_hills) = colvars(i)%ss0 712 hills_env%delta_s_history(i, n_hills) = colvars(i)%delta_s 713 ENDDO 714 hills_env%ww_history(n_hills) = ww 715 IF (wtcontrol) hills_env%invdt_history(n_hills) = invdt 716 717 CALL timestop(handle) 718 719 END SUBROUTINE add_hill_single 720 721! ************************************************************************************************** 722!> \brief Restart Hills Information 723!> \param ss_history ... 724!> \param delta_s_history ... 725!> \param ww_history ... 726!> \param ww ... 727!> \param n_hills ... 728!> \param n_colvar ... 729!> \param colvars ... 730!> \param metadyn_section ... 731!> \param invdt_history ... 732!> \author Teodoro Laino [tlaino] - University of Zurich - 10.2008 733! ************************************************************************************************** 734 SUBROUTINE restart_hills(ss_history, delta_s_history, ww_history, ww, & 735 n_hills, n_colvar, colvars, metadyn_section, invdt_history) 736 REAL(KIND=dp), DIMENSION(:, :), POINTER :: ss_history, delta_s_history 737 REAL(KIND=dp), DIMENSION(:), POINTER :: ww_history 738 REAL(KIND=dp) :: ww 739 INTEGER, INTENT(IN) :: n_hills, n_colvar 740 TYPE(metavar_type), DIMENSION(:), POINTER :: colvars 741 TYPE(section_vals_type), POINTER :: metadyn_section 742 REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: invdt_history 743 744 CHARACTER(len=*), PARAMETER :: routineN = 'restart_hills', routineP = moduleN//':'//routineN 745 746 INTEGER :: handle, i, j, ndum 747 LOGICAL :: explicit, wtcontrol 748 REAL(KIND=dp) :: rval 749 REAL(KIND=dp), DIMENSION(:), POINTER :: rvals 750 TYPE(section_vals_type), POINTER :: hills_history 751 752 CALL timeset(routineN, handle) 753 754 wtcontrol = PRESENT(invdt_history) 755 NULLIFY (rvals) 756 hills_history => section_vals_get_subs_vals(metadyn_section, "SPAWNED_HILLS_POS") 757 CALL section_vals_get(hills_history, explicit=explicit) 758 IF (explicit) THEN 759 CALL section_vals_val_get(hills_history, "_DEFAULT_KEYWORD_", n_rep_val=ndum) 760 ! ss_history, delta_s_history, ww_history, invdt_history : deallocate and reallocate with the proper size 761 DEALLOCATE (ss_history) 762 DEALLOCATE (delta_s_history) 763 DEALLOCATE (ww_history) 764 IF (wtcontrol) THEN 765 DEALLOCATE (invdt_history) 766 ENDIF 767 ! 768 CPASSERT(n_hills == ndum) 769 ALLOCATE (ss_history(n_colvar, n_hills)) 770 ALLOCATE (delta_s_history(n_colvar, n_hills)) 771 ALLOCATE (ww_history(n_hills)) 772 IF (wtcontrol) THEN 773 ALLOCATE (invdt_history(n_hills)) 774 ENDIF 775 ! 776 DO i = 1, n_hills 777 CALL section_vals_val_get(hills_history, "_DEFAULT_KEYWORD_", & 778 i_rep_val=i, r_vals=rvals) 779 CPASSERT(SIZE(rvals) == n_colvar) 780 ss_history(1:n_colvar, i) = rvals 781 END DO 782 ! 783 hills_history => section_vals_get_subs_vals(metadyn_section, "SPAWNED_HILLS_SCALE") 784 CALL section_vals_get(hills_history, explicit=explicit) 785 IF (explicit) THEN 786 ! delta_s_history 787 CALL section_vals_val_get(hills_history, "_DEFAULT_KEYWORD_", n_rep_val=ndum) 788 CPASSERT(n_hills == ndum) 789 DO i = 1, n_hills 790 CALL section_vals_val_get(hills_history, "_DEFAULT_KEYWORD_", & 791 i_rep_val=i, r_vals=rvals) 792 CPASSERT(SIZE(rvals) == n_colvar) 793 delta_s_history(1:n_colvar, i) = rvals 794 END DO 795 ELSE 796 CALL cp_warn(__LOCATION__, & 797 "Section SPAWNED_HILLS_SCALE is not present! Setting the scales of the"// & 798 "restarted hills according the parameters specified in the input file. ") 799 DO i = 1, n_hills 800 DO j = 1, n_colvar 801 delta_s_history(j, i) = colvars(i)%delta_s 802 END DO 803 END DO 804 END IF 805 ! 806 hills_history => section_vals_get_subs_vals(metadyn_section, "SPAWNED_HILLS_HEIGHT") 807 CALL section_vals_get(hills_history, explicit=explicit) 808 IF (explicit) THEN 809 ! ww_history 810 CALL section_vals_val_get(hills_history, "_DEFAULT_KEYWORD_", & 811 n_rep_val=ndum) 812 CPASSERT(n_hills == ndum) 813 DO i = 1, n_hills 814 CALL section_vals_val_get(hills_history, "_DEFAULT_KEYWORD_", & 815 i_rep_val=i, r_val=rval) 816 CPASSERT(SIZE(rvals) == n_colvar) 817 ww_history(i) = rval 818 END DO 819 ELSE 820 CALL cp_warn(__LOCATION__, & 821 "Section SPAWNED_HILLS_HEIGHT is not present! Setting the height of the"// & 822 " restarted hills according the parameters specified in the input file. ") 823 ww_history = ww 824 END IF 825 ! 826 hills_history => section_vals_get_subs_vals(metadyn_section, "SPAWNED_HILLS_INVDT") 827 CALL section_vals_get(hills_history, explicit=explicit) 828 IF (wtcontrol) THEN 829 IF (explicit) THEN 830 ! invdt_history 831 CALL section_vals_val_get(hills_history, "_DEFAULT_KEYWORD_", & 832 n_rep_val=ndum) 833 CPASSERT(n_hills == ndum) 834 DO i = 1, n_hills 835 CALL section_vals_val_get(hills_history, "_DEFAULT_KEYWORD_", & 836 i_rep_val=i, r_val=rval) 837 CPASSERT(SIZE(rvals) == n_colvar) 838 invdt_history(i) = rval 839 END DO 840 ELSE 841 CALL cp_warn(__LOCATION__, & 842 "Section SPAWNED_HILLS_INVDT is not present! Restarting from standard"// & 843 " metadynamics run i.e. setting 1/(Delta T) equal to zero. ") 844 invdt_history = 0._dp 845 END IF 846 ELSE 847 IF (explicit) THEN 848 CALL cp_abort(__LOCATION__, & 849 "Found section SPAWNED_HILLS_INVDT while restarting a standard metadynamics run..."// & 850 " Cannot restart metadynamics from well-tempered MetaD runs. ") 851 END IF 852 END IF 853 ENDIF 854 855 CALL timestop(handle) 856 857 END SUBROUTINE restart_hills 858 859! ************************************************************************************************** 860!> \brief Retrieves the iteration level for the metadynamics loop 861!> \param meta_env ... 862!> \param iter_nr ... 863!> \author Teodoro Laino [tlaino] - University of Zurich - 10.2008 864! ************************************************************************************************** 865 SUBROUTINE get_meta_iter_level(meta_env, iter_nr) 866 TYPE(meta_env_type), POINTER :: meta_env 867 INTEGER, INTENT(OUT) :: iter_nr 868 869 CHARACTER(len=*), PARAMETER :: routineN = 'get_meta_iter_level', & 870 routineP = moduleN//':'//routineN 871 872 IF (meta_env%do_multiple_walkers) THEN 873 iter_nr = meta_env%multiple_walkers%n_hills_local 874 ELSE 875 iter_nr = meta_env%hills_env%n_hills 876 END IF 877 878 END SUBROUTINE get_meta_iter_level 879 880! ************************************************************************************************** 881!> \brief ... 882!> \param meta_env ... 883!> \par History 884!> 11.2007 [created] [tlaino] 885!> \author Teodoro Laino - University of Zurich - 11.2007 886! ************************************************************************************************** 887 SUBROUTINE meta_walls(meta_env) 888 TYPE(meta_env_type), POINTER :: meta_env 889 890 CHARACTER(len=*), PARAMETER :: routineN = 'meta_walls', routineP = moduleN//':'//routineN 891 892 INTEGER :: ih, iwall 893 REAL(dp) :: ddp, delta_s, dfunc, diff_ss, dp2, & 894 efunc, ww 895 TYPE(metavar_type), DIMENSION(:), POINTER :: colvars 896 897 colvars => meta_env%metavar 898 ! Forces from the Walls 899 DO ih = 1, SIZE(colvars) 900 IF (colvars(ih)%do_wall) THEN 901 colvars(ih)%epot_walls = 0.0_dp 902 colvars(ih)%ff_walls = 0.0_dp 903 DO iwall = 1, SIZE(colvars(ih)%walls) 904 SELECT CASE (colvars(ih)%walls(iwall)%id_type) 905 CASE (do_wall_reflective, do_wall_none) 906 ! Do Nothing.. treated in the main metadyn function 907 CYCLE 908 CASE (do_wall_quadratic) 909 diff_ss = colvars(ih)%ss0 - colvars(ih)%walls(iwall)%pos 910 IF (colvars(ih)%periodic) THEN 911 ! The difference of a periodic COLVAR is always within [-pi,pi] 912 diff_ss = SIGN(1.0_dp, ASIN(SIN(diff_ss)))*ACOS(COS(diff_ss)) 913 END IF 914 efunc = colvars(ih)%walls(iwall)%k_quadratic*diff_ss**2 915 dfunc = 2.0_dp*colvars(ih)%walls(iwall)%k_quadratic*diff_ss 916 SELECT CASE (colvars(ih)%walls(iwall)%id_direction) 917 CASE (do_wall_p) 918 IF (diff_ss > 0.0_dp) THEN 919 colvars(ih)%ff_walls = colvars(ih)%ff_walls - dfunc 920 colvars(ih)%epot_walls = colvars(ih)%epot_walls + efunc 921 END IF 922 CASE (do_wall_m) 923 IF (diff_ss < 0.0_dp) THEN 924 colvars(ih)%ff_walls = colvars(ih)%ff_walls - dfunc 925 colvars(ih)%epot_walls = colvars(ih)%epot_walls + efunc 926 END IF 927 END SELECT 928 CASE (do_wall_quartic) 929 diff_ss = colvars(ih)%ss0 - colvars(ih)%walls(iwall)%pos0 930 IF (colvars(ih)%periodic) THEN 931 ! The difference of a periodic COLVAR is always within [-pi,pi] 932 diff_ss = SIGN(1.0_dp, ASIN(SIN(diff_ss)))*ACOS(COS(diff_ss)) 933 END IF 934 efunc = colvars(ih)%walls(iwall)%k_quartic*diff_ss*diff_ss**4 935 dfunc = 4.0_dp*colvars(ih)%walls(iwall)%k_quartic*diff_ss**3 936 SELECT CASE (colvars(ih)%walls(iwall)%id_direction) 937 CASE (do_wall_p) 938 IF (diff_ss > 0.0_dp) THEN 939 colvars(ih)%ff_walls = colvars(ih)%ff_walls - dfunc 940 colvars(ih)%epot_walls = colvars(ih)%epot_walls + efunc 941 END IF 942 CASE (do_wall_m) 943 IF (diff_ss < 0.0_dp) THEN 944 colvars(ih)%ff_walls = colvars(ih)%ff_walls - dfunc 945 colvars(ih)%epot_walls = colvars(ih)%epot_walls + efunc 946 END IF 947 END SELECT 948 CASE (do_wall_gaussian) 949 diff_ss = colvars(ih)%ss0 - colvars(ih)%walls(iwall)%pos 950 IF (colvars(ih)%periodic) THEN 951 ! The difference of a periodic COLVAR is always within [-pi,pi] 952 diff_ss = SIGN(1.0_dp, ASIN(SIN(diff_ss)))*ACOS(COS(diff_ss)) 953 END IF 954 ww = colvars(ih)%walls(iwall)%ww_gauss 955 delta_s = colvars(ih)%walls(iwall)%sigma_gauss 956 ddp = (diff_ss)/delta_s 957 dp2 = ddp**2 958 efunc = ww*EXP(-0.5_dp*dp2) 959 dfunc = -efunc*ddp/delta_s 960 colvars(ih)%ff_walls = colvars(ih)%ff_walls - dfunc 961 colvars(ih)%epot_walls = colvars(ih)%epot_walls + efunc 962 END SELECT 963 END DO 964 END IF 965 ENDDO 966 END SUBROUTINE meta_walls 967 968END MODULE metadynamics_utils 969