1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Perform a temperature accelarated hybrid monte carlo (TAHMC) run using QUICKSTEP 8!> \par History 9!> none 10!> \author Alin M Elena 11! ************************************************************************************************** 12MODULE tamc_run 13 14 USE atomic_kind_list_types, ONLY: atomic_kind_list_type 15 USE atomic_kind_types, ONLY: atomic_kind_type 16 USE averages_types, ONLY: average_quantities_type 17 USE barostat_types, ONLY: barostat_type,& 18 create_barostat_type,& 19 release_barostat_type 20 USE bibliography, ONLY: VandenCic2006 21 USE cell_types, ONLY: cell_type 22 USE colvar_methods, ONLY: colvar_eval_glob_f 23 USE colvar_types, ONLY: HBP_colvar_id,& 24 WC_colvar_id,& 25 colvar_p_type 26 USE constraint_fxd, ONLY: fix_atom_control 27 USE cp_external_control, ONLY: external_control 28 USE cp_log_handling, ONLY: cp_get_default_logger,& 29 cp_logger_get_default_io_unit,& 30 cp_logger_type 31 USE cp_output_handling, ONLY: cp_add_iter_level,& 32 cp_iterate,& 33 cp_p_file,& 34 cp_print_key_finished_output,& 35 cp_print_key_should_output,& 36 cp_print_key_unit_nr,& 37 cp_rm_iter_level 38 USE cp_para_types, ONLY: cp_para_env_type 39 USE cp_subsys_types, ONLY: cp_subsys_get,& 40 cp_subsys_type 41 USE cp_units, ONLY: cp_unit_from_cp2k 42 USE distribution_1d_types, ONLY: distribution_1d_type 43 USE force_env_methods, ONLY: force_env_calc_energy_force 44 USE force_env_types, ONLY: force_env_get,& 45 force_env_type 46 USE free_energy_types, ONLY: fe_env_create,& 47 free_energy_type 48 USE global_types, ONLY: global_environment_type 49 USE input_constants, ONLY: langevin_ensemble,& 50 npe_f_ensemble,& 51 npe_i_ensemble,& 52 nph_uniaxial_damped_ensemble,& 53 nph_uniaxial_ensemble,& 54 npt_f_ensemble,& 55 npt_i_ensemble,& 56 reftraj_ensemble 57 USE input_cp2k_check, ONLY: remove_restart_info 58 USE input_cp2k_restarts, ONLY: write_restart 59 USE input_section_types, ONLY: section_vals_get,& 60 section_vals_get_subs_vals,& 61 section_vals_remove_values,& 62 section_vals_type,& 63 section_vals_val_get,& 64 section_vals_val_set 65 USE kinds, ONLY: dp 66 USE machine, ONLY: m_walltime 67 USE mc_environment_types, ONLY: get_mc_env,& 68 mc_env_create,& 69 mc_env_release,& 70 mc_environment_type,& 71 set_mc_env 72 USE mc_misc, ONLY: mc_averages_create,& 73 mc_averages_release 74 USE mc_move_control, ONLY: init_mc_moves,& 75 mc_moves_release 76 USE mc_types, ONLY: get_mc_par,& 77 mc_averages_type,& 78 mc_ekin_type,& 79 mc_moves_type,& 80 mc_simpar_type,& 81 set_mc_par 82 USE md_ener_types, ONLY: create_md_ener,& 83 md_ener_type,& 84 release_md_ener 85 USE md_energies, ONLY: initialize_md_ener,& 86 md_energy 87 USE md_environment_types, ONLY: get_md_env,& 88 md_env_create,& 89 md_env_release,& 90 md_environment_type,& 91 set_md_env 92 USE md_run, ONLY: qs_mol_dyn 93 USE metadynamics_types, ONLY: meta_env_type,& 94 metavar_type,& 95 set_meta_env 96 USE molecule_kind_list_types, ONLY: molecule_kind_list_type 97 USE molecule_kind_types, ONLY: molecule_kind_type 98 USE molecule_list_types, ONLY: molecule_list_type 99 USE molecule_types, ONLY: global_constraint_type,& 100 molecule_type 101 USE parallel_rng_types, ONLY: UNIFORM,& 102 create_rng_stream,& 103 delete_rng_stream,& 104 next_random_number,& 105 rng_stream_type 106 USE particle_list_types, ONLY: particle_list_type 107 USE particle_types, ONLY: particle_type 108 USE physcon, ONLY: boltzmann,& 109 femtoseconds,& 110 joule,& 111 kelvin 112 USE qmmm_util, ONLY: apply_qmmm_walls_reflective 113 USE qs_environment_types, ONLY: get_qs_env 114 USE qs_scf_post_gpw, ONLY: scf_post_calculation_gpw 115 USE reference_manager, ONLY: cite_reference 116 USE reftraj_types, ONLY: create_reftraj,& 117 reftraj_type,& 118 release_reftraj 119 USE reftraj_util, ONLY: initialize_reftraj 120 USE simpar_methods, ONLY: read_md_section 121 USE simpar_types, ONLY: create_simpar_type,& 122 release_simpar_type,& 123 simpar_type 124 USE string_utilities, ONLY: str_comp 125 USE thermal_region_types, ONLY: release_thermal_regions,& 126 thermal_regions_type 127 USE thermal_region_utils, ONLY: create_thermal_regions 128 USE thermostat_methods, ONLY: create_thermostats 129 USE thermostat_types, ONLY: release_thermostats,& 130 thermostats_type 131 USE virial_methods, ONLY: virial_evaluate 132 USE virial_types, ONLY: virial_type 133 USE wiener_process, ONLY: create_wiener_process,& 134 create_wiener_process_cv 135!!!!! monte carlo part 136#include "../../base/base_uses.f90" 137 138 IMPLICIT NONE 139 140 PRIVATE 141 142 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'tamc_run' 143 144 PUBLIC :: qs_tamc 145 146CONTAINS 147 148! ************************************************************************************************** 149!> \brief Driver routine for TAHMC 150!> \param force_env ... 151!> \param globenv ... 152!> \param averages ... 153!> \author Alin M Elena 154!> \note it computes the forces using QuickStep. 155! ************************************************************************************************** 156 SUBROUTINE qs_tamc(force_env, globenv, averages) 157 158 TYPE(force_env_type), POINTER :: force_env 159 TYPE(global_environment_type), POINTER :: globenv 160 TYPE(average_quantities_type), OPTIONAL, POINTER :: averages 161 162 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_tamc', routineP = moduleN//':'//routineN 163 164 CHARACTER(LEN=20) :: ensemble 165 INTEGER :: handle, i, initialStep, iprint, isos, & 166 istep, j, md_stride, nmccycles, & 167 output_unit, rand2skip, run_type_id 168 INTEGER, POINTER :: itimes 169 LOGICAL :: check, explicit, my_rm_restart_info, & 170 save_mem, should_stop 171 REAL(KIND=dp) :: auxRandom, inittime, rval, temp, & 172 time_iter_start, time_iter_stop 173 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: An, fz, xieta, zbuff 174 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: r 175 REAL(KIND=dp), POINTER :: constant, time, used_time 176 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds 177 TYPE(barostat_type), POINTER :: barostat 178 TYPE(cell_type), POINTER :: cell 179 TYPE(cp_logger_type), POINTER :: logger 180 TYPE(cp_para_env_type), POINTER :: para_env 181 TYPE(cp_subsys_type), POINTER :: subsys, subsys_i 182 TYPE(distribution_1d_type), POINTER :: local_particles 183 TYPE(free_energy_type), POINTER :: fe_env 184 TYPE(mc_averages_type), POINTER :: MCaverages 185 TYPE(mc_environment_type), POINTER :: mc_env 186 TYPE(mc_moves_type), POINTER :: gmoves, moves 187 TYPE(mc_simpar_type), POINTER :: mc_par 188 TYPE(md_ener_type), POINTER :: md_ener 189 TYPE(md_environment_type), POINTER :: md_env 190 TYPE(meta_env_type), POINTER :: meta_env_saved 191 TYPE(particle_list_type), POINTER :: particles 192 TYPE(reftraj_type), POINTER :: reftraj 193 TYPE(rng_stream_type), POINTER :: rng_stream, rng_stream_mc 194 TYPE(section_vals_type), POINTER :: constraint_section, force_env_section, & 195 free_energy_section, fs_section, global_section, mc_section, md_section, motion_section, & 196 reftraj_section, subsys_section, work_section 197 TYPE(simpar_type), POINTER :: simpar 198 TYPE(thermal_regions_type), POINTER :: thermal_regions 199 TYPE(thermostats_type), POINTER :: thermostats 200 TYPE(virial_type), POINTER :: virial 201 202 initialStep = 0 203 inittime = 0.0_dp 204 205 CALL timeset(routineN, handle) 206 my_rm_restart_info = .TRUE. 207 NULLIFY (md_env, para_env, fs_section, virial) 208 para_env => force_env%para_env 209 motion_section => section_vals_get_subs_vals(force_env%root_section, "MOTION") 210 md_section => section_vals_get_subs_vals(motion_section, "MD") 211 212 ! Real call to MD driver - Low Level 213 CALL md_env_create(md_env, md_section, para_env, force_env=force_env) 214 IF (PRESENT(averages)) CALL set_md_env(md_env, averages=averages) 215 216 CPASSERT(ASSOCIATED(globenv)) 217 CPASSERT(ASSOCIATED(force_env)) 218 219 NULLIFY (particles, cell, simpar, itimes, used_time, subsys, & 220 md_ener, thermostats, barostat, reftraj, force_env_section, & 221 reftraj_section, work_section, atomic_kinds, & 222 local_particles, time, fe_env, free_energy_section, & 223 constraint_section, thermal_regions, subsys_i) 224 logger => cp_get_default_logger() 225 para_env => force_env%para_env 226 227 global_section => section_vals_get_subs_vals(force_env%root_section, "GLOBAL") 228 free_energy_section => section_vals_get_subs_vals(motion_section, "FREE_ENERGY") 229 constraint_section => section_vals_get_subs_vals(motion_section, "CONSTRAINT") 230 CALL section_vals_val_get(global_section, "SAVE_MEM", l_val=save_mem) 231 232 CALL section_vals_val_get(global_section, "RUN_TYPE", i_val=run_type_id) 233 234 CALL create_simpar_type(simpar) 235 force_env_section => force_env%force_env_section 236 subsys_section => section_vals_get_subs_vals(force_env_section, "SUBSYS") 237 CALL cp_add_iter_level(logger%iter_info, "MD") 238 CALL cp_iterate(logger%iter_info, iter_nr=initialStep) 239 ! Read MD section 240 CALL read_md_section(simpar, motion_section, md_section) 241 ! Setup print_keys 242 simpar%info_constraint = cp_print_key_unit_nr(logger, constraint_section, & 243 "CONSTRAINT_INFO", extension=".shakeLog", log_filename=.FALSE.) 244 simpar%lagrange_multipliers = cp_print_key_unit_nr(logger, constraint_section, & 245 "LAGRANGE_MULTIPLIERS", extension=".LagrangeMultLog", log_filename=.FALSE.) 246 simpar%dump_lm = BTEST(cp_print_key_should_output(logger%iter_info, constraint_section, & 247 "LAGRANGE_MULTIPLIERS"), cp_p_file) 248 249 ! Create the structure for the md energies 250 CALL create_md_ener(md_ener) 251 CALL set_md_env(md_env, md_ener=md_ener) 252 CALL release_md_ener(md_ener) 253 254 ! If requested setup Thermostats 255 CALL create_thermostats(thermostats, md_section, force_env, simpar, para_env, & 256 globenv, global_section) 257 258 ! If requested setup Barostat 259 CALL create_barostat_type(barostat, md_section, force_env, simpar, globenv) 260 261 ! If requested setup different thermal regions 262 CALL create_thermal_regions(thermal_regions, md_section, simpar, force_env) 263 264 CALL set_md_env(md_env, thermostats=thermostats, barostat=barostat, thermal_regions=thermal_regions) 265 266 IF (simpar%ensemble == reftraj_ensemble) THEN 267 reftraj_section => section_vals_get_subs_vals(md_section, "REFTRAJ") 268 CALL create_reftraj(reftraj, reftraj_section, para_env) 269 CALL set_md_env(md_env, reftraj=reftraj) 270 CALL release_reftraj(reftraj) 271 END IF 272 273 CALL force_env_get(force_env, subsys=subsys, cell=cell, & 274 force_env_section=force_env_section) 275 276 ! Set V0 if needed 277 IF (simpar%ensemble == nph_uniaxial_ensemble .OR. simpar%ensemble == nph_uniaxial_damped_ensemble) THEN 278 IF (simpar%v0 == 0._dp) simpar%v0 = cell%deth 279 ENDIF 280 281 ! Initialize velocities possibly applying constraints at the zeroth MD step 282! ! ! CALL section_vals_val_get(motion_section,"PRINT%RESTART%SPLIT_RESTART_FILE",& 283! ! ! l_val=write_binary_restart_file) 284!! let us see if this created all the trouble 285! CALL setup_velocities(force_env,simpar,globenv,md_env,md_section,constraint_section, & 286! write_binary_restart_file) 287 288 ! Setup Free Energy Calculation (if required) 289 CALL fe_env_create(fe_env, free_energy_section) 290 CALL set_md_env(md_env=md_env, simpar=simpar, fe_env=fe_env, cell=cell, & 291 force_env=force_env) 292 293 ! Possibly initialize Wiener processes 294 IF (simpar%ensemble == langevin_ensemble) CALL create_wiener_process(md_env) 295 time_iter_start = m_walltime() 296 297 CALL get_md_env(md_env, force_env=force_env, itimes=itimes, constant=constant, & 298 md_ener=md_ener, t=time, used_time=used_time) 299 300 ! Attach the time counter of the meta_env to the one of the MD 301 CALL set_meta_env(force_env%meta_env, time=time) 302 ! Initialize the md_ener structure 303 304 force_env%meta_env%dt = force_env%meta_env%zdt 305 CALL initialize_md_ener(md_ener, force_env, simpar) 306! force_env%meta_env%dt=force_env%meta_env%zdt 307 308!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! MC setup up 309!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 310 311 NULLIFY (mc_env, mc_par, rng_stream_mc, MCaverages) 312 313 CALL section_vals_get(force_env_section, n_repetition=isos) 314 CPASSERT(isos == 1) 315! set some values...will use get_globenv if that ever comes around 316 317! initialize the random numbers 318! IF (para_env%ionode) THEN 319 CALL create_rng_stream(rng_stream=rng_stream_mc, & 320 name="Random numbers for monte carlo acc/rej", & 321 distribution_type=UNIFORM) 322! ENDIF 323!!!!! this shoudl go in a routine hmc_read 324 325 NULLIFY (mc_section) 326 ALLOCATE (mc_par) 327 328 mc_section => section_vals_get_subs_vals(force_env%root_section, & 329 "MOTION%MC") 330 CALL section_vals_val_get(mc_section, "ENSEMBLE", & 331 c_val=ensemble) 332 CPASSERT(str_comp(ensemble, "TRADITIONAL")) 333 CALL section_vals_val_get(mc_section, "NSTEP", & 334 i_val=nmccycles) 335 CPASSERT(nmccycles > 0) 336 CALL section_vals_val_get(mc_section, "IPRINT", & 337 i_val=iprint) 338 CALL section_vals_val_get(mc_section, "RANDOMTOSKIP", i_val=rand2skip) 339 CPASSERT(rand2skip >= 0) 340 temp = cp_unit_from_cp2k(simpar%temp_ext, "K") 341! 342 343 CALL set_mc_par(mc_par, ensemble=ensemble, nstep=nmccycles, iprint=iprint, temperature=temp, & 344 beta=1.0_dp/temp/boltzmann*joule, exp_max_val=0.9_dp*LOG(HUGE(0.0_dp)), & 345 exp_min_val=0.9_dp*LOG(TINY(0.0_dp)), max_val=HUGE(0.0_dp), min_val=0.0_dp, & 346 source=para_env%source, group=para_env%group, ionode=para_env%ionode, rand2skip=rand2skip) 347 348 output_unit = cp_logger_get_default_io_unit(logger) 349 IF (output_unit > 0) THEN 350 WRITE (output_unit, '(a,a)') "HMC| Hybrid Monte Carlo Scheme " 351 WRITE (output_unit, '(a,a)') "HMC| Ensemble ", ADJUSTL(ensemble) 352 WRITE (output_unit, '(a,i0)') "HMC| MC Cycles ", nmccycles 353 WRITE (output_unit, '(a,i0,a)') "HMC| Print every ", iprint, " cycles" 354 WRITE (output_unit, '(a,i0)') "HMC| Number of random numbers to skip ", rand2skip 355 WRITE (output_unit, '(a,f16.8,a)') "HMC| Temperature ", temp, "K" 356 ENDIF 357 358 CALL force_env_get(force_env, subsys=subsys) 359 360 CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, & 361 particles=particles, virial=virial) 362 363 DO i = 1, rand2skip 364 auxRandom = next_random_number(rng_stream_mc) 365 DO j = 1, 3*SIZE(particles%els) 366 auxRandom = next_random_number(globenv%gaussian_rng_stream) 367 ENDDO 368 ENDDO 369 370 CALL mc_env_create(mc_env) 371 CALL set_mc_env(mc_env, mc_par=mc_par, force_env=force_env) 372!!!!!!!end mc setup 373 374 ! Check for ensembles requiring the stress tensor - takes into account the possibility for 375 ! multiple force_evals 376 IF ((simpar%ensemble == npt_i_ensemble) .OR. & 377 (simpar%ensemble == npt_f_ensemble) .OR. & 378 (simpar%ensemble == npe_f_ensemble) .OR. & 379 (simpar%ensemble == npe_i_ensemble) .OR. & 380 (simpar%ensemble == nph_uniaxial_ensemble) .OR. & 381 (simpar%ensemble == nph_uniaxial_damped_ensemble)) THEN 382 check = virial%pv_availability 383 IF (.NOT. check) & 384 CALL cp_abort(__LOCATION__, & 385 "Virial evaluation not requested for this run in the input file! "// & 386 "You may consider to switch on the virial evaluation with the keyword: STRESS_TENSOR."// & 387 "Be sure the method you are using can compute the virial!") 388 IF (ASSOCIATED(force_env%sub_force_env)) THEN 389 DO i = 1, SIZE(force_env%sub_force_env) 390 IF (ASSOCIATED(force_env%sub_force_env(i)%force_env)) THEN 391 CALL force_env_get(force_env%sub_force_env(i)%force_env, subsys=subsys_i) 392 CALL cp_subsys_get(subsys_i, virial=virial) 393 check = check .AND. virial%pv_availability 394 END IF 395 END DO 396 END IF 397 IF (.NOT. check) & 398 CALL cp_abort(__LOCATION__, & 399 "Virial evaluation not requested for all the force_eval sections present in"// & 400 " the input file! You have to switch on the virial evaluation with the keyword: STRESS_TENSOR "// & 401 " in each force_eval section. Be sure the method you are using can compute the virial!") 402 END IF 403 404 ! Computing Forces at zero MD step 405 IF (simpar%ensemble /= reftraj_ensemble) THEN 406 CALL section_vals_val_get(md_section, "STEP_START_VAL", i_val=itimes) 407 CALL section_vals_val_get(md_section, "TIME_START_VAL", r_val=time) 408 CALL section_vals_val_get(md_section, "ECONS_START_VAL", r_val=constant) 409 CALL section_vals_val_set(md_section, "STEP_START_VAL", i_val=initialStep) 410 CALL section_vals_val_set(md_section, "TIME_START_VAL", r_val=inittime) 411 initialStep = itimes 412 CALL cp_iterate(logger%iter_info, iter_nr=itimes) 413 IF (save_mem) THEN 414 work_section => section_vals_get_subs_vals(subsys_section, "VELOCITY") 415 CALL section_vals_remove_values(work_section) 416 work_section => section_vals_get_subs_vals(subsys_section, "SHELL_VELOCITY") 417 CALL section_vals_remove_values(work_section) 418 work_section => section_vals_get_subs_vals(subsys_section, "CORE_VELOCITY") 419 CALL section_vals_remove_values(work_section) 420 END IF 421 422! CALL force_env_calc_energy_force (force_env, calc_force=.TRUE.) 423 meta_env_saved => force_env%meta_env 424 NULLIFY (force_env%meta_env) 425 CALL force_env_calc_energy_force(force_env, calc_force=.FALSE.) 426 force_env%meta_env => meta_env_saved 427 428 IF (ASSOCIATED(force_env%qs_env)) THEN 429! force_env%qs_env%sim_time=time 430! force_env%qs_env%sim_step=itimes 431 force_env%qs_env%sim_time = 0.0_dp 432 force_env%qs_env%sim_step = 0 433 END IF 434 ! Warm-up engines for metadynamics 435 IF (ASSOCIATED(force_env%meta_env)) THEN 436 IF (force_env%meta_env%langevin) THEN 437 CALL create_wiener_process_cv(force_env%meta_env) 438 NULLIFY (rng_stream) 439 DO j = 1, (rand2skip - 1)/nmccycles 440 DO i = 1, force_env%meta_env%n_colvar 441 rng_stream => force_env%meta_env%rng(i)%stream 442 auxRandom = next_random_number(rng_stream) 443 auxRandom = next_random_number(rng_stream) 444 ENDDO 445 ENDDO 446 ENDIF 447! IF (force_env%meta_env%well_tempered) THEN 448! force_env%meta_env%wttemperature = simpar%temp_ext 449! IF (force_env%meta_env%wtgamma>EPSILON(1._dp)) THEN 450! dummy=force_env%meta_env%wttemperature*(force_env%meta_env%wtgamma-1._dp) 451! IF (force_env%meta_env%delta_t>EPSILON(1._dp)) THEN 452! check=ABS(force_env%meta_env%delta_t-dummy)<1.E+3_dp*EPSILON(1._dp) 453! IF(.NOT.check) CALL cp_abort(__LOCATION__,& 454! "Inconsistency between DELTA_T and WTGAMMA (both specified):"//& 455! " please, verify that DELTA_T=(WTGAMMA-1)*TEMPERATURE") 456! ELSE 457! force_env%meta_env%delta_t = dummy 458! ENDIF 459! ELSE 460! force_env%meta_env%wtgamma = 1._dp & 461! + force_env%meta_env%delta_t/force_env%meta_env%wttemperature 462! ENDIF 463! force_env%meta_env%invdt = 1._dp/force_env%meta_env%delta_t 464! ENDIF 465 CALL tamc_force(force_env) 466! CALL metadyn_write_colvar(force_env) 467 ENDIF 468 469 IF (simpar%do_respa) THEN 470 CALL force_env_calc_energy_force(force_env%sub_force_env(1)%force_env, & 471 calc_force=.TRUE.) 472 END IF 473 474! CALL force_env_get( force_env, subsys=subsys) 475! 476! CALL cp_subsys_get(subsys,atomic_kinds=atomic_kinds,local_particles=local_particles,& 477! particles=particles) 478 479 CALL virial_evaluate(atomic_kinds%els, particles%els, local_particles, & 480 virial, force_env%para_env%group) 481 482 CALL md_energy(md_env, md_ener) 483! CALL md_write_output(md_env) !inits the print env at itimes == 0 also writes trajectories 484 md_stride = 1 485 ELSE 486 CALL get_md_env(md_env, reftraj=reftraj) 487 CALL initialize_reftraj(reftraj, reftraj_section, md_env) 488 itimes = reftraj%info%first_snapshot - 1 489 md_stride = reftraj%info%stride 490 END IF 491 492 CALL cp_print_key_finished_output(simpar%info_constraint, logger, & 493 constraint_section, "CONSTRAINT_INFO") 494 CALL cp_print_key_finished_output(simpar%lagrange_multipliers, logger, & 495 constraint_section, "LAGRANGE_MULTIPLIERS") 496 CALL init_mc_moves(moves) 497 CALL init_mc_moves(gmoves) 498 ALLOCATE (r(1:3, SIZE(particles%els))) 499! ALLOCATE (r_old(1:3,size(particles%els))) 500 CALL mc_averages_create(MCaverages) 501 !!!!! some more buffers 502 ! Allocate random number for Langevin Thermostat acting on COLVARS 503 ALLOCATE (xieta(2*force_env%meta_env%n_colvar)) 504 xieta(:) = 0.0_dp 505 ALLOCATE (An(force_env%meta_env%n_colvar)) 506 An(:) = 0.0_dp 507 ALLOCATE (fz(force_env%meta_env%n_colvar)) 508 fz(:) = 0.0_dp 509 ALLOCATE (zbuff(2*force_env%meta_env%n_colvar)) 510 zbuff(:) = 0.0_dp 511 512 IF (output_unit > 0) THEN 513 WRITE (output_unit, '(a)') "HMC|==== Initial average forces" 514 ENDIF 515 CALL metadyn_write_colvar_header(force_env) 516 moves%hmc%attempts = 0 517 moves%hmc%successes = 0 518 gmoves%hmc%attempts = 0 519 gmoves%hmc%successes = 0 520 IF (initialStep == 0) THEN 521!!! if we come from a restart we shall properly compute the average force 522!!! read the average force up to now 523 DO i = 1, force_env%meta_env%n_colvar 524 fs_section => section_vals_get_subs_vals(force_env%meta_env%metadyn_section, "EXT_LAGRANGE_FS") 525 CALL section_vals_get(fs_section, explicit=explicit) 526 IF (explicit) THEN 527 CALL section_vals_val_get(fs_section, "_DEFAULT_KEYWORD_", & 528 i_rep_val=i, r_val=rval) 529 fz(i) = rval*rand2skip 530 END IF 531 ENDDO 532 533 CALL HMCsampler(globenv, force_env, MCaverages, r, mc_par, moves, gmoves, rng_stream_mc, output_unit, & 534 fz, zbuff, nskip=rand2skip) 535 CALL cp_iterate(logger%iter_info, last=.FALSE., iter_nr=0) 536 CALL section_vals_val_set(mc_section, "RANDOMTOSKIP", i_val=rand2skip + nmccycles) 537 CALL write_restart(md_env=md_env, root_section=force_env%root_section) 538 ENDIF 539 IF (output_unit > 0) THEN 540 WRITE (output_unit, '(a)') "HMC|==== end initial average forces" 541 ENDIF 542! call set_md_env(md_env, init=.FALSE.) 543 544 CALL metadyn_write_colvar(force_env) 545 546 DO istep = 1, force_env%meta_env%TAMCSteps 547 ! Increase counters 548 itimes = itimes + 1 549 time = time + force_env%meta_env%dt 550 IF (output_unit > 0) THEN 551 WRITE (output_unit, '(a)') "HMC|===================================" 552 WRITE (output_unit, '(a,1x,i0)') "HMC| on z step ", istep 553 ENDIF 554 !needed when electric field fields are applied 555 IF (ASSOCIATED(force_env%qs_env)) THEN 556 force_env%qs_env%sim_time = time 557 force_env%qs_env%sim_step = itimes 558 force_env%meta_env%time = force_env%qs_env%sim_time 559 END IF 560 561 CALL cp_iterate(logger%iter_info, last=(istep == force_env%meta_env%TAMCSteps), iter_nr=itimes) 562 ! Open possible Shake output units 563 simpar%info_constraint = cp_print_key_unit_nr(logger, constraint_section, "CONSTRAINT_INFO", & 564 extension=".shakeLog", log_filename=.FALSE.) 565 simpar%lagrange_multipliers = cp_print_key_unit_nr( & 566 logger, constraint_section, & 567 "LAGRANGE_MULTIPLIERS", extension=".LagrangeMultLog", log_filename=.FALSE.) 568 simpar%dump_lm = BTEST(cp_print_key_should_output(logger%iter_info, constraint_section, & 569 "LAGRANGE_MULTIPLIERS"), cp_p_file) 570 571 ! Velocity Verlet Integrator 572 573 moves%hmc%attempts = 0 574 moves%hmc%successes = 0 575 CALL langevinVEC(md_env, globenv, mc_env, moves, gmoves, r, & 576 rng_stream_mc, xieta, An, fz, MCaverages, zbuff) 577 578 ! Close Shake output if requested... 579 CALL cp_print_key_finished_output(simpar%info_constraint, logger, & 580 constraint_section, "CONSTRAINT_INFO") 581 CALL cp_print_key_finished_output(simpar%lagrange_multipliers, logger, & 582 constraint_section, "LAGRANGE_MULTIPLIERS") 583 CALL cp_iterate(logger%iter_info, iter_nr=initialStep) 584 CALL metadyn_write_colvar(force_env) 585 ! Free Energy calculation 586! CALL free_energy_evaluate(md_env,should_stop,free_energy_section) 587 588 ![AME:UB] IF (should_stop) EXIT 589 590 ! Test for <PROJECT_NAME>.EXIT_MD or for WALL_TIME to exit 591 ! Default: 592 ! IF so we don't overwrite the restart or append to the trajectory 593 ! because the execution could in principle stop inside the SCF where energy 594 ! and forces are not converged. 595 ! But: 596 ! You can force to print the last step (for example if the method used 597 ! to compute energy and forces is not SCF based) activating the print_key 598 ! MOTION%MD%PRINT%FORCE_LAST. 599 CALL external_control(should_stop, "MD", globenv=globenv) 600 IF (should_stop) THEN 601 CALL cp_iterate(logger%iter_info, last=.TRUE., iter_nr=itimes) 602! CALL md_output(md_env,md_section,force_env%root_section,should_stop) 603 EXIT 604 END IF 605 606! IF(simpar%ensemble /= reftraj_ensemble) THEN 607! CALL md_energy(md_env, md_ener) 608! CALL temperature_control(simpar, md_env, md_ener, force_env, logger) 609! CALL comvel_control(md_ener, force_env, md_section, logger) 610! CALL angvel_control(md_ener, force_env, md_section, logger) 611! ELSE 612! CALL md_ener_reftraj(md_env, md_ener) 613! END IF 614 615 time_iter_stop = m_walltime() 616 used_time = time_iter_stop - time_iter_start 617 time_iter_start = time_iter_stop 618 619!!!!! this writes the restart... 620! CALL md_output(md_env,md_section,force_env%root_section,should_stop) 621 622! IF(simpar%ensemble == reftraj_ensemble ) THEN 623! CALL write_output_reftraj(md_env) 624! END IF 625 626 IF (output_unit > 0) THEN 627 WRITE (output_unit, '(a,1x,i0)') "HMC| end z step ", istep 628 WRITE (output_unit, '(a)') "HMC|===================================" 629 ENDIF 630 END DO 631 CALL cp_iterate(logger%iter_info, last=.TRUE., iter_nr=itimes) 632 force_env%qs_env%sim_time = 0.0_dp 633 force_env%qs_env%sim_step = 0 634 rand2skip = rand2skip + nmccycles*force_env%meta_env%TAMCSteps 635 IF (initialStep == 0) rand2skip = rand2skip + nmccycles 636 CALL section_vals_val_set(mc_section, "RANDOMTOSKIP", i_val=rand2skip) 637 638 CALL write_restart(md_env=md_env, root_section=force_env%root_section) 639! if we need the final kinetic energy for Hybrid Monte Carlo 640! hmc_ekin%final_ekin=md_ener%ekin 641 642 ! Remove the iteration level 643 CALL cp_rm_iter_level(logger%iter_info, "MD") 644 645 ! Deallocate Thermostats and Barostats 646 CALL release_thermostats(thermostats) 647 CALL release_barostat_type(barostat) 648 CALL release_simpar_type(simpar) 649 CALL release_thermal_regions(thermal_regions) 650 651 CALL md_env_release(md_env) 652 ! Clean restartable sections.. 653 IF (my_rm_restart_info) CALL remove_restart_info(force_env%root_section) 654! IF (para_env%ionode) THEN 655 CALL delete_rng_stream(rng_stream_mc) 656! ENDIF 657 CALL MC_ENV_RELEASE(mc_env) 658 DEALLOCATE (mc_par) 659 CALL MC_MOVES_RELEASE(moves) 660 CALL MC_MOVES_RELEASE(gmoves) 661 DEALLOCATE (r) 662! DEALLOCATE(r_old) 663 DEALLOCATE (xieta) 664 DEALLOCATE (An) 665 DEALLOCATE (fz) 666 DEALLOCATE (zbuff) 667 CALL mc_averages_release(MCaverages) 668 CALL timestop(handle) 669 670 END SUBROUTINE qs_tamc 671 672! ************************************************************************************************** 673!> \brief Propagates velocities for z half a step 674!> \param force_env ... 675!> \param An ... 676!> \author Alin M Elena 677!> \note Vanden-Eijnden Ciccotti C.Phys.Letter 429 (2006) 310-316 678! ************************************************************************************************** 679 SUBROUTINE tamc_velocities_colvar(force_env, An) 680 TYPE(force_env_type), POINTER :: force_env 681 REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: An 682 683 CHARACTER(len=*), PARAMETER :: routineN = 'tamc_velocities_colvar', & 684 routineP = moduleN//':'//routineN 685 686 INTEGER :: handle, i_c 687 REAL(kind=dp) :: dt, fft, sigma 688 TYPE(cp_logger_type), POINTER :: logger 689 TYPE(meta_env_type), POINTER :: meta_env 690 TYPE(metavar_type), POINTER :: cv 691 692 NULLIFY (logger, meta_env, cv) 693 meta_env => force_env%meta_env 694 CALL timeset(routineN, handle) 695 logger => cp_get_default_logger() 696 ! Add citation 697 IF (meta_env%langevin) CALL cite_reference(VandenCic2006) 698 dt = meta_env%dt 699 700 ! Evolve Velocities 701 meta_env%epot_walls = 0.0_dp 702 DO i_c = 1, meta_env%n_colvar 703 cv => meta_env%metavar(i_c) 704 fft = cv%ff_s + cv%ff_hills 705 sigma = SQRT((meta_env%temp_wanted*kelvin)*2.0_dp*(boltzmann/joule)*cv%gamma/cv%mass) 706 cv%vvp = cv%vvp + 0.5_dp*dt*(fft/cv%mass - cv%gamma*cv%vvp)*(1.0_dp - 0.25_dp*dt*cv%gamma) + An(i_c) 707 meta_env%epot_walls = meta_env%epot_walls + cv%epot_walls 708 ENDDO 709 CALL timestop(handle) 710 END SUBROUTINE tamc_velocities_colvar 711 712! ************************************************************************************************** 713!> \brief propagates z one step 714!> \param force_env ... 715!> \param xieta ... 716!> \author Alin M Elena 717!> \note Vanden-Eijnden Ciccotti C.Phys.Letter 429 (2006) 310-316 718! ************************************************************************************************** 719 SUBROUTINE tamc_position_colvar(force_env, xieta) 720 TYPE(force_env_type), POINTER :: force_env 721 REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: xieta 722 723 CHARACTER(len=*), PARAMETER :: routineN = 'tamc_position_colvar', & 724 routineP = moduleN//':'//routineN 725 726 INTEGER :: handle, i_c 727 REAL(kind=dp) :: dt, sigma 728 TYPE(cp_logger_type), POINTER :: logger 729 TYPE(meta_env_type), POINTER :: meta_env 730 TYPE(metavar_type), POINTER :: cv 731 732 NULLIFY (logger, meta_env, cv) 733 meta_env => force_env%meta_env 734! IF (.NOT.ASSOCIATED(meta_env)) RETURN 735 736 CALL timeset(routineN, handle) 737 logger => cp_get_default_logger() 738 739 ! Add citation 740 IF (meta_env%langevin) CALL cite_reference(VandenCic2006) 741 dt = meta_env%dt 742 743 ! Update of ss0 744 DO i_c = 1, meta_env%n_colvar 745 cv => meta_env%metavar(i_c) 746 sigma = SQRT((meta_env%temp_wanted*kelvin)*2.0_dp*(boltzmann/joule)*cv%gamma/cv%mass) 747! cv%ss0 =cv%ss0 +dt*cv%vvp 748 cv%ss0 = cv%ss0 + dt*cv%vvp + dt*SQRT(dt/12.0_dp)*sigma*xieta(i_c + meta_env%n_colvar) 749 IF (cv%periodic) THEN 750 ! A periodic COLVAR is always within [-pi,pi] 751 cv%ss0 = SIGN(1.0_dp, ASIN(SIN(cv%ss0)))*ACOS(COS(cv%ss0)) 752 END IF 753 ENDDO 754 CALL timestop(handle) 755 756 END SUBROUTINE tamc_position_colvar 757 758! ************************************************************************************************** 759!> \brief Computes forces on z 760!> #details also can be used to get the potenzial evergy of z 761!> \param force_env ... 762!> \param zpot ... 763!> \author Alin M Elena 764! ************************************************************************************************** 765 SUBROUTINE tamc_force(force_env, zpot) 766 TYPE(force_env_type), POINTER :: force_env 767 REAL(KIND=dp), INTENT(inout), OPTIONAL :: zpot 768 769 CHARACTER(len=*), PARAMETER :: routineN = 'tamc_force', routineP = moduleN//':'//routineN 770 771 INTEGER :: handle, i, i_c, icolvar, ii 772 LOGICAL :: explicit 773 REAL(kind=dp) :: diff_ss, dt, rval 774 TYPE(colvar_p_type), DIMENSION(:), POINTER :: colvar_p 775 TYPE(cp_logger_type), POINTER :: logger 776 TYPE(cp_subsys_type), POINTER :: subsys 777 TYPE(meta_env_type), POINTER :: meta_env 778 TYPE(metavar_type), POINTER :: cv 779 TYPE(particle_list_type), POINTER :: particles 780 TYPE(section_vals_type), POINTER :: ss0_section, ss_section, vvp_section 781 782 NULLIFY (logger, meta_env) 783 meta_env => force_env%meta_env 784! IF (.NOT.ASSOCIATED(meta_env)) RETURN 785 786 CALL timeset(routineN, handle) 787 logger => cp_get_default_logger() 788 NULLIFY (colvar_p, subsys, cv, ss0_section, vvp_section, ss_section) 789 CALL force_env_get(force_env, subsys=subsys) 790 791 dt = meta_env%dt 792 IF (.NOT. meta_env%restart) meta_env%n_steps = meta_env%n_steps + 1 793 ! compute ss and the derivative of ss with respect to the atomic positions 794 DO i_c = 1, meta_env%n_colvar 795 cv => meta_env%metavar(i_c) 796 icolvar = cv%icolvar 797 CALL colvar_eval_glob_f(icolvar, force_env) 798 cv%ss = subsys%colvar_p(icolvar)%colvar%ss 799 ! Restart for Extended Lagrangian Metadynamics 800 IF (meta_env%restart) THEN 801 ! Initialize the position of the collective variable in the extended lagrange 802 ss0_section => section_vals_get_subs_vals(meta_env%metadyn_section, "EXT_LAGRANGE_SS0") 803 CALL section_vals_get(ss0_section, explicit=explicit) 804 IF (explicit) THEN 805 CALL section_vals_val_get(ss0_section, "_DEFAULT_KEYWORD_", & 806 i_rep_val=i_c, r_val=rval) 807 cv%ss0 = rval 808 ELSE 809 cv%ss0 = cv%ss 810 END IF 811 vvp_section => section_vals_get_subs_vals(meta_env%metadyn_section, "EXT_LAGRANGE_VVP") 812 CALL section_vals_get(vvp_section, explicit=explicit) 813 IF (explicit) THEN 814 CALL setup_velocities_z(force_env) 815 CALL section_vals_val_get(vvp_section, "_DEFAULT_KEYWORD_", & 816 i_rep_val=i_c, r_val=rval) 817 cv%vvp = rval 818 ELSE 819 CALL setup_velocities_z(force_env) 820 ENDIF 821 ss_section => section_vals_get_subs_vals(meta_env%metadyn_section, "EXT_LAGRANGE_SS") 822 CALL section_vals_get(ss_section, explicit=explicit) 823 IF (explicit) THEN 824 CALL section_vals_val_get(ss_section, "_DEFAULT_KEYWORD_", & 825 i_rep_val=i_c, r_val=rval) 826 cv%ss = rval 827 END IF 828 END IF 829 ! 830 ENDDO 831 ! forces on the atoms 832 NULLIFY (particles) 833 CALL cp_subsys_get(subsys, colvar_p=colvar_p, & 834 particles=particles) 835 836 meta_env%restart = .FALSE. 837 meta_env%epot_s = 0.0_dp 838 meta_env%epot_walls = 0.0_dp 839 DO i_c = 1, meta_env%n_colvar 840 cv => meta_env%metavar(i_c) 841 diff_ss = cv%ss - cv%ss0 842 IF (cv%periodic) THEN 843 ! The difference of a periodic COLVAR is always within [-pi,pi] 844 diff_ss = SIGN(1.0_dp, ASIN(SIN(diff_ss)))*ACOS(COS(diff_ss)) 845 END IF 846 cv%epot_s = 0.5_dp*cv%lambda*diff_ss*diff_ss 847 cv%ff_s = cv%lambda*(diff_ss) 848 meta_env%epot_s = meta_env%epot_s + cv%epot_s 849 icolvar = cv%icolvar 850 851 DO ii = 1, colvar_p(icolvar)%colvar%n_atom_s 852 i = colvar_p(icolvar)%colvar%i_atom(ii) 853 particles%els(i)%f = particles%els(i)%f - cv%ff_s*colvar_p(icolvar)%colvar%dsdr(:, ii) 854 ENDDO 855 856 ENDDO 857 IF (PRESENT(zpot)) zpot = meta_env%epot_s 858 CALL fix_atom_control(force_env) 859 860 CALL timestop(handle) 861 END SUBROUTINE tamc_force 862 863! ************************************************************************************************** 864!> \brief propagates one time step both z systems and samples the x system 865!> \param md_env ... 866!> \param globenv ... 867!> \param mc_env ... 868!> \param moves ... 869!> \param gmoves ... 870!> \param r ... 871!> \param rng_stream_mc ... 872!> \param xieta ... 873!> \param An ... 874!> \param fz ... 875!> \param averages ... 876!> \param zbuff ... 877!> \author Alin M Elena 878! ************************************************************************************************** 879 SUBROUTINE langevinVEC(md_env, globenv, mc_env, moves, gmoves, r, & 880 rng_stream_mc, xieta, An, fz, averages, zbuff) 881 882 TYPE(md_environment_type), POINTER :: md_env 883 TYPE(global_environment_type), POINTER :: globenv 884 TYPE(mc_environment_type), POINTER :: mc_env 885 TYPE(mc_moves_type), POINTER :: moves, gmoves 886 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: r 887 TYPE(rng_stream_type), POINTER :: rng_stream_mc 888 REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: xieta, An, fz 889 TYPE(mc_averages_type), INTENT(INOUT), POINTER :: averages 890 REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: zbuff 891 892 CHARACTER(len=*), PARAMETER :: routineN = 'langevinVEC', routineP = moduleN//':'//routineN 893 894 INTEGER :: iprint, ivar, nparticle, nparticle_kind, & 895 nstep, output_unit 896 REAL(KIND=dp) :: dt, gamma, mass, sigma 897 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds 898 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 899 TYPE(cell_type), POINTER :: cell 900 TYPE(cp_logger_type), POINTER :: logger 901 TYPE(cp_para_env_type), POINTER :: para_env 902 TYPE(cp_subsys_type), POINTER :: subsys 903 TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles 904 TYPE(force_env_type), POINTER :: force_env 905 TYPE(global_constraint_type), POINTER :: gci 906 TYPE(mc_simpar_type), POINTER :: mc_par 907 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds 908 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set 909 TYPE(molecule_list_type), POINTER :: molecules 910 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 911 TYPE(particle_list_type), POINTER :: particles 912 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 913 TYPE(rng_stream_type), POINTER :: rng_stream 914 TYPE(simpar_type), POINTER :: simpar 915 TYPE(virial_type), POINTER :: virial 916 917 NULLIFY (logger, mc_par) 918 logger => cp_get_default_logger() 919 output_unit = cp_logger_get_default_io_unit(logger) 920 921 NULLIFY (rng_stream) 922! quantitites to be nullified for the get_md_env 923 NULLIFY (simpar, force_env, para_env) 924! quantities to be nullified for the force_env_get environment 925 NULLIFY (subsys, cell) 926! quantitites to be nullified for the cp_subsys_get 927 NULLIFY (atomic_kinds, local_particles, particles, local_molecules, molecules, molecule_kinds, gci) 928 929 CALL get_md_env(md_env=md_env, simpar=simpar, force_env=force_env, & 930 para_env=para_env) 931 CALL get_mc_env(mc_env, mc_par=mc_par) 932 CALL get_mc_par(mc_par, nstep=nstep, iprint=iprint) 933 934 dt = simpar%dt 935 CALL force_env_get(force_env=force_env, subsys=subsys, cell=cell) 936 937!!!! this bit should vanish once I understand what the hell is with it 938 939! ! Do some checks on coordinates and box 940 CALL apply_qmmm_walls_reflective(force_env) 941 942 CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, & 943 particles=particles, local_molecules=local_molecules, molecules=molecules, & 944 molecule_kinds=molecule_kinds, gci=gci, virial=virial) 945 946 nparticle_kind = atomic_kinds%n_els 947 atomic_kind_set => atomic_kinds%els 948 molecule_kind_set => molecule_kinds%els 949 950 nparticle = particles%n_els 951 particle_set => particles%els 952 molecule_set => molecules%els 953 CPASSERT(ASSOCIATED(force_env%meta_env)) 954 CPASSERT(force_env%meta_env%langevin) 955 ! *** Velocity Verlet for Langevin *** v(t)--> v(t+1/2) 956 !!!!!! noise xi is in the first half, eta in the second half 957 DO ivar = 1, force_env%meta_env%n_colvar 958 rng_stream => force_env%meta_env%rng(ivar)%stream 959 xieta(ivar) = next_random_number(rng_stream) 960 xieta(ivar + force_env%meta_env%n_colvar) = next_random_number(rng_stream) 961 gamma = force_env%meta_env%metavar(ivar)%gamma 962 mass = force_env%meta_env%metavar(ivar)%mass 963 sigma = SQRT((force_env%meta_env%temp_wanted*kelvin)*2.0_dp*(boltzmann/joule)*gamma/mass) 964 An(ivar) = 0.5_dp*SQRT(dt)*sigma*(xieta(ivar)*(1.0_dp - 0.5_dp*dt*gamma) - & 965 dt*gamma*xieta(ivar + force_env%meta_env%n_colvar)/SQRT(12.0_dp)) 966 ENDDO 967! *** Velocity Verlet for Langeving *** v(t)--> v(t+1/2) 968 CALL tamc_velocities_colvar(force_env, An) 969! *** Velocity Verlet for Langevin S(t)->S(t+1) 970 CALL tamc_position_colvar(force_env, xieta) 971!!!!! start zHMC sampler 972 CALL HMCsampler(globenv, force_env, averages, r, mc_par, moves, gmoves, rng_stream_mc, output_unit, fz, zbuff) 973 974! CALL final_mc_write(mc_par,tmp_moves,& 975! output_unit,energy_check,& 976! initial_energy,final_energy,& 977! averages) 978 979!!!!!!!!!!!!!!!!!!!! end zHMC sampler 980 ! *** Velocity Verlet for Langeving *** v(t+1/2)--> v(t+1) 981 CALL tamc_velocities_colvar(force_env, An) 982! CALL virial_evaluate ( atomic_kind_set, particle_set, & 983! local_particles, virial, para_env%group) 984 985 END SUBROUTINE langevinVEC 986 987! ************************************************************************************************** 988!> \brief Driver routin for the canonical sampler using modified HMC 989!> \param globenv ... 990!> \param force_env ... 991!> \param averages ... 992!> \param r ... 993!> \param mc_par ... 994!> \param moves ... 995!> \param gmoves ... 996!> \param rng_stream_mc ... 997!> \param output_unit ... 998!> \param fz ... 999!> \param zbuff ... 1000!> \param nskip ... 1001!> \author Alin M Elena 1002!> \note at the end of this routine %ff_s will contain mean force 1003! ************************************************************************************************** 1004 1005 SUBROUTINE HMCsampler(globenv, force_env, averages, r, mc_par, moves, gmoves, rng_stream_mc, output_unit, & 1006 fz, zbuff, nskip) 1007 TYPE(global_environment_type), POINTER :: globenv 1008 TYPE(force_env_type), POINTER :: force_env 1009 TYPE(mc_averages_type), POINTER :: averages 1010 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: r 1011 TYPE(mc_simpar_type), POINTER :: mc_par 1012 TYPE(mc_moves_type), POINTER :: moves, gmoves 1013 TYPE(rng_stream_type), POINTER :: rng_stream_mc 1014 INTEGER, INTENT(IN) :: output_unit 1015 REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: fz, zbuff 1016 INTEGER, INTENT(IN), OPTIONAL :: nskip 1017 1018 INTEGER :: i, iprint, ishift, it1, j, nsamples, & 1019 nstep 1020 REAL(KIND=dp) :: energy_check, old_epx, old_epz, t1 1021 TYPE(meta_env_type), POINTER :: meta_env_saved 1022 1023 IF (PRESENT(nskip)) THEN 1024 nsamples = nskip 1025 ishift = nskip 1026 ELSE 1027 nsamples = 0 1028 fz = 0.0_dp 1029 ishift = 0 1030 ENDIF 1031! lrestart = .false. 1032! if (present(logger) .and. present(iter)) THEN 1033! lrestart=.true. 1034! ENDIF 1035 CALL get_mc_par(mc_par, nstep=nstep, iprint=iprint) 1036 meta_env_saved => force_env%meta_env 1037 NULLIFY (force_env%meta_env) 1038 CALL force_env_get(force_env, potential_energy=old_epx) 1039 force_env%meta_env => meta_env_saved 1040 1041 old_epz = force_env%meta_env%epot_s 1042!!! average energy will be wrong on restarts 1043 averages%ave_energy = 0.0_dp 1044 t1 = force_env%qs_env%sim_time 1045 it1 = force_env%qs_env%sim_step 1046 IF (output_unit > 0) THEN 1047 WRITE (output_unit, '(a,l4)') "HMC|restart? ", force_env%meta_env%restart 1048 WRITE (output_unit, '(a,3(f16.8,1x))') & 1049 "HMC|Ep, Epx, Epz ", old_epx + force_env%meta_env%epot_s, old_epx, force_env%meta_env%epot_s 1050 WRITE (output_unit, '(a)') "#HMC| No | z.. | theta.. | ff_z... | ff_z/n |" 1051 ENDIF 1052 DO i = 1, nstep 1053 IF (MOD(i, iprint) == 0 .AND. (output_unit > 0)) THEN 1054 WRITE (output_unit, '(a,1x,i0)') "HMC|========== On Monte Carlo cycle ", i + ishift 1055 WRITE (output_unit, '(a)') "HMC| Attempting a minitrajectory move" 1056 WRITE (output_unit, '(a,1x,i0)') "HMC| start mini-trajectory", i + ishift 1057 WRITE (output_unit, '(a,1x,i0,1x)', advance="no") "#HMC|0 ", i + ishift 1058 DO j = 1, force_env%meta_env%n_colvar 1059 WRITE (output_unit, '(f16.8,1x,f16.8,1x,f16.8)', advance="no") force_env%meta_env%metavar(j)%ss0, & 1060 force_env%meta_env%metavar(j)%ss, & 1061 force_env%meta_env%metavar(j)%ff_s !,fz(j)/real(i+ishift,dp) 1062 ENDDO 1063 WRITE (output_unit, *) 1064 ENDIF 1065 1066 CALL mc_hmc_move(mc_par, force_env, globenv, moves, gmoves, old_epx, old_epz, energy_check, & 1067 r, output_unit, rng_stream_mc, zbuff) 1068 ! check averages... 1069 ! force average for z needed too... 1070 averages%ave_energy = averages%ave_energy*REAL(i - 1, dp)/REAL(i, dp) + & 1071 old_epx/REAL(i, dp) 1072 DO j = 1, force_env%meta_env%n_colvar 1073 fz(j) = fz(j) + force_env%meta_env%metavar(j)%ff_s 1074 ENDDO 1075 IF (output_unit > 0) THEN 1076 WRITE (output_unit, '(a,1x,i0)') "HMC|end mini-trajectory", i + ishift 1077!!!!!!!! this prints z and theta(x) --ss0,ss-- needed to determine an acceptable k then 1078 ! the instanteneous force and some instanteneous average for force 1079 WRITE (output_unit, '(a,1x,i0,1x)', advance="no") "#HMC|1 ", i + ishift 1080 DO j = 1, force_env%meta_env%n_colvar 1081 WRITE (output_unit, '(f16.8,1x,f16.8,1x,f16.8,1x,f16.8)', advance="no") force_env%meta_env%metavar(j)%ss0, & 1082 force_env%meta_env%metavar(j)%ss, & 1083 force_env%meta_env%metavar(j)%ff_s, fz(j)/REAL(i + ishift, dp) 1084 ENDDO 1085 WRITE (output_unit, *) 1086 ENDIF 1087 nsamples = nsamples + 1 1088 IF (MOD(i, iprint) == 0 .AND. (output_unit > 0)) THEN 1089 WRITE (output_unit, '(a,f16.8)') "HMC| Running average for potential energy ", averages%ave_energy 1090 WRITE (output_unit, '(a,1x,i0)') "HMC|======== End Monte Carlo cycle ", i + ishift 1091 ENDIF 1092! IF (lrestart) THEN 1093! k=nstep/5 1094! IF(MOD(i,k) == 0) THEN 1095! force_env%qs_env%sim_time=t1 1096! force_env%qs_env%sim_step=it1 1097! DO j=1,force_env%meta_env%n_colvar 1098! force_env%meta_env%metavar(j)%ff_s=fz(j)/real(i+ishift,dp) 1099! ENDDO 1100! ! CALL cp_iterate(logger%iter_info,last=.FALSE.,iter_nr=-1) 1101! CALL section_vals_val_set(mcsec,"RANDOMTOSKIP",i_val=i+ishift) 1102! CALL write_restart(md_env=mdenv,root_section=force_env%root_section) 1103! ! CALL cp_iterate(logger%iter_info,last=.FALSE.,iter_nr=iter) 1104! ENDIF 1105! ENDIF 1106 ENDDO 1107 force_env%qs_env%sim_time = t1 1108 force_env%qs_env%sim_step = it1 1109 IF (output_unit > 0) THEN 1110 WRITE (output_unit, '(a,i0,a,i0,a,f16.8)') "HMC| local acceptance ratio: ", moves%hmc%successes, "/", & 1111 moves%hmc%attempts, "=", REAL(moves%hmc%successes, dp)/REAL(moves%hmc%attempts, dp) 1112 WRITE (output_unit, '(a,i0,a,i0,a,f16.8)') "HMC| global acceptance ratio: ", gmoves%hmc%successes, "/", & 1113 gmoves%hmc%attempts, "=", REAL(gmoves%hmc%successes, dp)/REAL(gmoves%hmc%attempts, dp) 1114 ENDIF 1115 !average force 1116 DO j = 1, force_env%meta_env%n_colvar 1117 force_env%meta_env%metavar(j)%ff_s = fz(j)/nsamples 1118 ENDDO 1119 END SUBROUTINE HMCsampler 1120 1121! ************************************************************************************************** 1122!> \brief performs a hybrid Monte Carlo move 1123!> \param mc_par ... 1124!> \param force_env ... 1125!> \param globenv ... 1126!> \param moves ... 1127!> \param gmoves ... 1128!> \param old_epx ... 1129!> \param old_epz ... 1130!> \param energy_check ... 1131!> \param r ... 1132!> \param output_unit ... 1133!> \param rng_stream ... 1134!> \param zbuff ... 1135!> \author Alin M Elena 1136!> \note It runs a NVE trajectory, followed by localisation and accepts rejects 1137!> using the biased Hamiltonian, rather than the traditional guiding Hamiltonian 1138! ************************************************************************************************** 1139 SUBROUTINE mc_hmc_move(mc_par, force_env, globenv, moves, gmoves, old_epx, old_epz, & 1140 energy_check, r, output_unit, rng_stream, zbuff) 1141 1142 TYPE(mc_simpar_type), POINTER :: mc_par 1143 TYPE(force_env_type), POINTER :: force_env 1144 TYPE(global_environment_type), POINTER :: globenv 1145 TYPE(mc_moves_type), POINTER :: moves, gmoves 1146 REAL(KIND=dp), INTENT(INOUT) :: old_epx, old_epz, energy_check 1147 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: r 1148 INTEGER, INTENT(IN) :: output_unit 1149 TYPE(rng_stream_type), POINTER :: rng_stream 1150 REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: zbuff 1151 1152 CHARACTER(LEN=*), PARAMETER :: routineN = 'mc_hmc_move', routineP = moduleN//':'//routineN 1153 1154 INTEGER :: group, handle, iatom, j, nAtoms, source 1155 LOGICAL :: ionode, localise 1156 REAL(KIND=dp) :: BETA, energy_term, exp_max_val, & 1157 exp_min_val, new_energy, new_epx, & 1158 new_epz, rand, value, w 1159 TYPE(cp_subsys_type), POINTER :: oldsys 1160 TYPE(mc_ekin_type), POINTER :: hmc_ekin 1161 TYPE(meta_env_type), POINTER :: meta_env_saved 1162 TYPE(particle_list_type), POINTER :: particles_set 1163 TYPE(section_vals_type), POINTER :: dft_section, input 1164 1165! begin the timing of the subroutine 1166 1167 CALL timeset(routineN, handle) 1168 1169 CALL get_qs_env(force_env%qs_env, input=input) 1170 dft_section => section_vals_get_subs_vals(input, "DFT") 1171 1172! get a bunch of stuff from mc_par 1173 CALL get_mc_par(mc_par, ionode=ionode, & 1174 BETA=BETA, exp_max_val=exp_max_val, & 1175 exp_min_val=exp_min_val, source=source, group=group) 1176 1177! nullify some pointers 1178! NULLIFY(particles_set,oldsys,hmc_ekin) 1179 NULLIFY (particles_set, oldsys, meta_env_saved, hmc_ekin) 1180 ! now let's grab the particle positions 1181 CALL force_env_get(force_env, subsys=oldsys) 1182 CALL cp_subsys_get(oldsys, particles=particles_set) 1183 nAtoms = SIZE(particles_set%els) 1184! do some allocation 1185 1186 ALLOCATE (hmc_ekin) 1187 1188! record the attempt 1189 moves%hmc%attempts = moves%hmc%attempts + 1 1190 gmoves%hmc%attempts = gmoves%hmc%attempts + 1 1191 1192! save the old coordinates just in case we need to go back 1193 DO iatom = 1, nAtoms 1194 r(1:3, iatom) = particles_set%els(iatom)%r(1:3) 1195 ENDDO 1196 localise = .TRUE. 1197! the same for collective variables data should be made,ss first half and ff_s the last half 1198 DO j = 1, force_env%meta_env%n_colvar 1199 zbuff(j) = force_env%meta_env%metavar(j)%ss 1200 zbuff(j + force_env%meta_env%n_colvar) = force_env%meta_env%metavar(j)%ff_s 1201 IF ((oldsys%colvar_p(force_env%meta_env%metavar(j)%icolvar)%colvar%type_id == HBP_colvar_id) .OR. & 1202 (oldsys%colvar_p(force_env%meta_env%metavar(j)%icolvar)%colvar%type_id == WC_colvar_id)) THEN 1203 localise = .FALSE. 1204 ENDIF 1205 ENDDO 1206 1207! now run the MD simulation 1208 meta_env_saved => force_env%meta_env 1209 NULLIFY (force_env%meta_env) 1210 force_env%qs_env%sim_time = 0.0_dp 1211 force_env%qs_env%sim_step = 0 1212 IF (.NOT. localise) THEN 1213 CALL section_vals_val_set(dft_section, "LOCALIZE%_SECTION_PARAMETERS_", l_val=.FALSE.) 1214 ENDIF 1215 CALL qs_mol_dyn(force_env, globenv, hmc_e_initial=hmc_ekin%initial_ekin, hmc_e_final=hmc_ekin%final_ekin) 1216 IF (.NOT. localise) THEN 1217 CALL section_vals_val_set(dft_section, "LOCALIZE%_SECTION_PARAMETERS_", l_val=.TRUE.) 1218 CALL scf_post_calculation_gpw(force_env%qs_env) 1219 ENDIF 1220 1221 CALL force_env_get(force_env, potential_energy=new_epx) 1222 1223 force_env%meta_env => meta_env_saved 1224 CALL tamc_force(force_env, zpot=new_epz) 1225 new_energy = new_epx + new_epz 1226 IF (output_unit > 0) THEN 1227 WRITE (output_unit, '(a,4(f16.8,1x))') & 1228 "HMC|old Ep, Ekx, Epz, Epx ", old_epx + old_epz, hmc_ekin%initial_ekin, old_epz, old_epx 1229 WRITE (output_unit, '(a,4(f16.8,1x))') "HMC|new Ep, Ekx, Epz, Epx ", new_energy, hmc_ekin%final_ekin, new_epz, new_epx 1230 ENDIF 1231 energy_term = new_energy - old_epx - old_epz + hmc_ekin%final_ekin - hmc_ekin%initial_ekin 1232 1233 value = -BETA*(energy_term) 1234! to prevent overflows 1235 IF (value > exp_max_val) THEN 1236 w = 10.0_dp 1237 ELSEIF (value < exp_min_val) THEN 1238 w = 0.0_dp 1239 ELSE 1240 w = EXP(value) 1241 ENDIF 1242 1243 rand = next_random_number(rng_stream) 1244 IF (rand < w) THEN 1245! accept the move 1246 moves%hmc%successes = moves%hmc%successes + 1 1247 gmoves%hmc%successes = gmoves%hmc%successes + 1 1248! update energies 1249 energy_check = energy_check + (new_energy - old_epx - old_epz) 1250 old_epx = new_epx 1251 old_epz = new_epz 1252 ELSE 1253! reset the cell and particle positions 1254 DO iatom = 1, nAtoms 1255 particles_set%els(iatom)%r(1:3) = r(1:3, iatom) 1256 ENDDO 1257 DO j = 1, force_env%meta_env%n_colvar 1258 force_env%meta_env%metavar(j)%ss = zbuff(j) 1259 force_env%meta_env%metavar(j)%ff_s = zbuff(j + force_env%meta_env%n_colvar) 1260 ENDDO 1261 1262 ENDIF 1263 1264 DEALLOCATE (hmc_ekin) 1265 1266! end the timing 1267 CALL timestop(handle) 1268 1269 END SUBROUTINE mc_hmc_move 1270 1271! ************************************************************************************************** 1272!> \brief ... 1273!> \param force_env ... 1274! ************************************************************************************************** 1275 SUBROUTINE metadyn_write_colvar_header(force_env) 1276 TYPE(force_env_type), POINTER :: force_env 1277 1278 CHARACTER(len=*), PARAMETER :: routineN = 'metadyn_write_colvar_header', & 1279 routineP = moduleN//':'//routineN 1280 1281 CHARACTER(len=100) :: aux, fmt 1282 CHARACTER(len=255) :: label1, label2, label3, label4, label5, & 1283 label6 1284 INTEGER :: handle, i, iw, m 1285 TYPE(cp_logger_type), POINTER :: logger 1286 TYPE(meta_env_type), POINTER :: meta_env 1287 1288 NULLIFY (logger, meta_env) 1289 meta_env => force_env%meta_env 1290 IF (.NOT. ASSOCIATED(meta_env)) RETURN 1291 1292 CALL timeset(routineN, handle) 1293 logger => cp_get_default_logger() 1294 1295 iw = cp_print_key_unit_nr(logger, meta_env%metadyn_section, & 1296 "PRINT%COLVAR", extension=".metadynLog") 1297 IF (iw > 0) THEN 1298 label1 = "" 1299 label2 = "" 1300 label3 = "" 1301 label4 = "" 1302 label5 = "" 1303 label6 = "" 1304 DO i = 1, meta_env%n_colvar 1305 WRITE (aux, '(a,i0)') "z_", i 1306 label1 = TRIM(label1)//TRIM(aux) 1307 m = 15*i - LEN_TRIM(label1) - 1 1308 label1 = TRIM(label1)//REPEAT(" ", m)//"|" 1309 WRITE (aux, '(a,i0)') "Theta_", i 1310 label2 = TRIM(label2)//TRIM(aux) 1311 m = 15*i - LEN_TRIM(label2) - 1 1312 label2 = TRIM(label2)//REPEAT(" ", m)//"|" 1313 WRITE (aux, '(a,i0)') "F_z", i 1314 label3 = TRIM(label3)//TRIM(aux) 1315 m = 15*i - LEN_TRIM(label3) - 1 1316 label3 = TRIM(label3)//REPEAT(" ", m)//"|" 1317 WRITE (aux, '(a,i0)') "F_h", i 1318 label4 = TRIM(label4)//TRIM(aux) 1319 m = 15*i - LEN_TRIM(label4) - 1 1320 label4 = TRIM(label4)//REPEAT(" ", m)//"|" 1321 WRITE (aux, '(a,i0)') "F_w", i 1322 label5 = TRIM(label5)//TRIM(aux) 1323 m = 15*i - LEN_TRIM(label5) - 1 1324 label5 = TRIM(label5)//REPEAT(" ", m)//"|" 1325 WRITE (aux, '(a,i0)') "v_z", i 1326 label6 = TRIM(label6)//TRIM(aux) 1327 m = 15*i - LEN_TRIM(label6) - 1 1328 label6 = TRIM(label6)//REPEAT(" ", m)//"|" 1329 ENDDO 1330 WRITE (fmt, '("(a17,6a",i0 ,",4a15)")') meta_env%n_colvar*15 1331 WRITE (iw, TRIM(fmt)) "#Time[fs] |", & 1332 TRIM(label1), & 1333 TRIM(label2), & 1334 TRIM(label3), & 1335 TRIM(label4), & 1336 TRIM(label5), & 1337 TRIM(label6), & 1338 "Epot_z |", & 1339 "Ene hills |", & 1340 "Epot walls |", & 1341 "Temperature |" 1342 1343 END IF 1344 CALL cp_print_key_finished_output(iw, logger, meta_env%metadyn_section, & 1345 "PRINT%COLVAR") 1346 1347 CALL timestop(handle) 1348 1349 END SUBROUTINE metadyn_write_colvar_header 1350 1351! ************************************************************************************************** 1352!> \brief ... 1353!> \param force_env ... 1354! ************************************************************************************************** 1355 SUBROUTINE metadyn_write_colvar(force_env) 1356 TYPE(force_env_type), POINTER :: force_env 1357 1358 CHARACTER(len=*), PARAMETER :: routineN = 'metadyn_write_colvar', & 1359 routineP = moduleN//':'//routineN 1360 1361 INTEGER :: handle, i, i_c, iw 1362 REAL(KIND=dp) :: temp 1363 TYPE(cp_logger_type), POINTER :: logger 1364 TYPE(meta_env_type), POINTER :: meta_env 1365 TYPE(metavar_type), POINTER :: cv 1366 1367 NULLIFY (logger, meta_env, cv) 1368 meta_env => force_env%meta_env 1369 IF (.NOT. ASSOCIATED(meta_env)) RETURN 1370 1371 CALL timeset(routineN, handle) 1372 logger => cp_get_default_logger() 1373 1374 IF (meta_env%langevin) THEN 1375 meta_env%ekin_s = 0.0_dp 1376! meta_env%epot_s = 0.0_dp 1377 DO i_c = 1, meta_env%n_colvar 1378 cv => meta_env%metavar(i_c) 1379 meta_env%ekin_s = meta_env%ekin_s + 0.5_dp*cv%mass*cv%vvp**2 1380 ENDDO 1381 END IF 1382 1383 ! write COLVAR file 1384 iw = cp_print_key_unit_nr(logger, meta_env%metadyn_section, & 1385 "PRINT%COLVAR", extension=".metadynLog") 1386 IF (iw > 0) THEN 1387 IF (meta_env%extended_lagrange) THEN 1388 WRITE (iw, '(f16.8,70f15.8)') meta_env%time*femtoseconds, & 1389 (meta_env%metavar(i)%ss0, i=1, meta_env%n_colvar), & 1390 (meta_env%metavar(i)%ss, i=1, meta_env%n_colvar), & 1391 (meta_env%metavar(i)%ff_s, i=1, meta_env%n_colvar), & 1392 (meta_env%metavar(i)%ff_hills, i=1, meta_env%n_colvar), & 1393 (meta_env%metavar(i)%ff_walls, i=1, meta_env%n_colvar), & 1394 (meta_env%metavar(i)%vvp, i=1, meta_env%n_colvar), & 1395 meta_env%epot_s, & 1396 meta_env%hills_env%energy, & 1397 meta_env%epot_walls, & 1398 (meta_env%ekin_s)*2.0_dp/(REAL(meta_env%n_colvar, KIND=dp))*kelvin 1399 ELSE 1400 WRITE (iw, '(f16.8,40f13.5)') meta_env%time*femtoseconds, & 1401 (meta_env%metavar(i)%ss0, i=1, meta_env%n_colvar), & 1402 (meta_env%metavar(i)%ff_hills, i=1, meta_env%n_colvar), & 1403 (meta_env%metavar(i)%ff_walls, i=1, meta_env%n_colvar), & 1404 meta_env%hills_env%energy, & 1405 meta_env%epot_walls 1406 END IF 1407 END IF 1408 CALL cp_print_key_finished_output(iw, logger, meta_env%metadyn_section, & 1409 "PRINT%COLVAR") 1410 ! Temperature for COLVAR 1411 IF (meta_env%extended_lagrange) THEN 1412 temp = meta_env%ekin_s*2.0_dp/(REAL(meta_env%n_colvar, KIND=dp))*kelvin 1413 meta_env%avg_temp = (meta_env%avg_temp*REAL(meta_env%n_steps, KIND=dp) + & 1414 temp)/REAL(meta_env%n_steps + 1, KIND=dp) 1415 iw = cp_print_key_unit_nr(logger, meta_env%metadyn_section, & 1416 "PRINT%TEMPERATURE_COLVAR", extension=".metadynLog") 1417 IF (iw > 0) THEN 1418 WRITE (iw, '(T2,79("-"))') 1419 WRITE (iw, '( A,T51,f10.2,T71,f10.2)') ' COLVARS INSTANTANEOUS/AVERAGE TEMPERATURE ', & 1420 temp, meta_env%avg_temp 1421 WRITE (iw, '(T2,79("-"))') 1422 ENDIF 1423 CALL cp_print_key_finished_output(iw, logger, meta_env%metadyn_section, & 1424 "PRINT%TEMPERATURE_COLVAR") 1425 END IF 1426 CALL timestop(handle) 1427 1428 END SUBROUTINE metadyn_write_colvar 1429 1430! ************************************************************************************************** 1431!> \brief ... 1432!> \param force_env ... 1433! ************************************************************************************************** 1434 SUBROUTINE setup_velocities_z(force_env) 1435 TYPE(force_env_type), POINTER :: force_env 1436 1437 INTEGER :: i_c 1438 REAL(kind=dp) :: ekin_w, fac_t 1439 TYPE(meta_env_type), POINTER :: meta_env 1440 TYPE(metavar_type), POINTER :: cv 1441 1442 NULLIFY (meta_env) 1443 meta_env => force_env%meta_env 1444 meta_env%ekin_s = 0.0_dp 1445 DO i_c = 1, meta_env%n_colvar 1446 cv => meta_env%metavar(i_c) 1447 cv%vvp = next_random_number(force_env%globenv%gaussian_rng_stream) 1448 meta_env%ekin_s = meta_env%ekin_s + 0.5_dp*cv%mass*cv%vvp**2 1449 END DO 1450 ekin_w = 0.5_dp*meta_env%temp_wanted*REAL(meta_env%n_colvar, KIND=dp) 1451 fac_t = SQRT(ekin_w/MAX(meta_env%ekin_s, 1.0E-8_dp)) 1452 DO i_c = 1, meta_env%n_colvar 1453 cv => meta_env%metavar(i_c) 1454 cv%vvp = cv%vvp*fac_t 1455 ENDDO 1456 END SUBROUTINE setup_velocities_z 1457END MODULE tamc_run 1458