1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Methods to performs a path integral run 8!> \author fawzi 9!> \par History 10!> 02.2005 created [fawzi] 11!> 11.2006 modified so it might actually work [hforbert] 12!> 10.2015 added RPMD propagator 13!> 10.2015 added exact harmonic integrator [Felix Uhl] 14!> \note quick & dirty rewrite of my python program 15! ************************************************************************************************** 16MODULE pint_methods 17 18 USE atomic_kind_list_types, ONLY: atomic_kind_list_type 19 USE atomic_kind_types, ONLY: atomic_kind_type,& 20 get_atomic_kind 21 USE bibliography, ONLY: Brieuc2016,& 22 Ceriotti2010,& 23 Ceriotti2012,& 24 cite_reference 25 USE cell_types, ONLY: cell_type 26 USE constraint, ONLY: rattle_control,& 27 shake_control,& 28 shake_update_targets 29 USE constraint_util, ONLY: getold 30 USE cp_external_control, ONLY: external_control 31 USE cp_log_handling, ONLY: cp_get_default_logger,& 32 cp_logger_get_default_io_unit,& 33 cp_logger_type,& 34 cp_to_string 35 USE cp_output_handling, ONLY: cp_add_iter_level,& 36 cp_iterate,& 37 cp_p_file,& 38 cp_print_key_finished_output,& 39 cp_print_key_should_output,& 40 cp_print_key_unit_nr,& 41 cp_rm_iter_level 42 USE cp_para_types, ONLY: cp_para_env_type 43 USE cp_subsys_types, ONLY: cp_subsys_get,& 44 cp_subsys_type 45 USE cp_units, ONLY: cp_unit_from_cp2k,& 46 cp_unit_to_cp2k 47 USE distribution_1d_types, ONLY: distribution_1d_type 48 USE f77_interface, ONLY: f_env_add_defaults,& 49 f_env_rm_defaults,& 50 f_env_type 51 USE force_env_types, ONLY: force_env_get 52 USE gle_system_dynamics, ONLY: gle_cholesky_stab,& 53 gle_matrix_exp,& 54 restart_gle 55 USE gle_system_types, ONLY: gle_dealloc,& 56 gle_init,& 57 gle_thermo_create 58 USE global_types, ONLY: global_environment_type 59 USE helium_interactions, ONLY: helium_intpot_scan 60 USE helium_io, ONLY: helium_write_cubefile 61 USE helium_methods, ONLY: helium_create,& 62 helium_init,& 63 helium_release 64 USE helium_sampling, ONLY: helium_do_run,& 65 helium_step 66 USE helium_types, ONLY: helium_solvent_p_type 67 USE input_constants, ONLY: integrate_exact,& 68 integrate_numeric,& 69 propagator_rpmd,& 70 transformation_normal,& 71 transformation_stage 72 USE input_cp2k_restarts, ONLY: write_restart 73 USE input_section_types, ONLY: & 74 section_type, section_vals_get, section_vals_get_subs_vals, section_vals_release, & 75 section_vals_retain, section_vals_type, section_vals_val_get, section_vals_val_set, & 76 section_vals_val_unset 77 USE kinds, ONLY: default_path_length,& 78 default_string_length,& 79 dp 80 USE machine, ONLY: m_flush,& 81 m_walltime 82 USE mathconstants, ONLY: twopi 83 USE mathlib, ONLY: gcd 84 USE message_passing, ONLY: mp_bcast,& 85 mp_comm_self 86 USE molecule_kind_list_types, ONLY: molecule_kind_list_type 87 USE molecule_kind_types, ONLY: molecule_kind_type 88 USE molecule_list_types, ONLY: molecule_list_type 89 USE molecule_types, ONLY: global_constraint_type,& 90 molecule_type 91 USE parallel_rng_types, ONLY: GAUSSIAN,& 92 create_rng_stream,& 93 delete_rng_stream,& 94 next_random_number,& 95 next_rng_seed,& 96 rng_stream_type 97 USE particle_list_types, ONLY: particle_list_type 98 USE particle_types, ONLY: particle_type 99 USE pint_gle, ONLY: pint_calc_gle_energy,& 100 pint_gle_init,& 101 pint_gle_step 102 USE pint_io, ONLY: pint_write_action,& 103 pint_write_centroids,& 104 pint_write_com,& 105 pint_write_ener,& 106 pint_write_line,& 107 pint_write_rgyr,& 108 pint_write_step_info,& 109 pint_write_trajectory 110 USE pint_normalmode, ONLY: normalmode_calc_uf_h,& 111 normalmode_env_create,& 112 normalmode_init_masses,& 113 normalmode_release 114 USE pint_piglet, ONLY: pint_calc_piglet_energy,& 115 pint_piglet_create,& 116 pint_piglet_init,& 117 pint_piglet_release,& 118 pint_piglet_step 119 USE pint_pile, ONLY: pint_calc_pile_energy,& 120 pint_pile_init,& 121 pint_pile_release,& 122 pint_pile_step 123 USE pint_public, ONLY: pint_levy_walk 124 USE pint_qtb, ONLY: pint_calc_qtb_energy,& 125 pint_qtb_init,& 126 pint_qtb_release,& 127 pint_qtb_step 128 USE pint_staging, ONLY: staging_calc_uf_h,& 129 staging_env_create,& 130 staging_init_masses,& 131 staging_release 132 USE pint_transformations, ONLY: pint_f2uf,& 133 pint_u2x,& 134 pint_x2u 135 USE pint_types, ONLY: & 136 e_conserved_id, e_kin_thermo_id, e_kin_virial_id, e_potential_id, pint_env_type, & 137 thermostat_gle, thermostat_none, thermostat_nose, thermostat_piglet, thermostat_pile, & 138 thermostat_qtb 139 USE replica_methods, ONLY: rep_env_calc_e_f,& 140 rep_env_create 141 USE replica_types, ONLY: rep_env_release,& 142 replica_env_type 143 USE simpar_types, ONLY: create_simpar_type,& 144 release_simpar_type 145#include "../base/base_uses.f90" 146 147 IMPLICIT NONE 148 PRIVATE 149 150 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .TRUE. 151 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pint_methods' 152 INTEGER, SAVE, PRIVATE :: last_pint_id = 0 153 154 PUBLIC :: do_pint_run 155 156CONTAINS 157 158! *************************************************************************** 159!> \brief Create a path integral environment 160!> \param pint_env ... 161!> \param input ... 162!> \param input_declaration ... 163!> \param para_env ... 164!> \par History 165!> Fixed some bugs [hforbert] 166!> Added normal mode transformation [hforbert] 167!> 10.2015 Added RPMD propagator and harmonic integrator [Felix Uhl] 168!> 10.2018 Added centroid constraints [cschran+rperez] 169!> \author fawzi 170!> \note Might return an unassociated pointer in parallel on the processors 171!> that are not needed. 172! ************************************************************************************************** 173 SUBROUTINE pint_create(pint_env, input, input_declaration, para_env) 174 175 TYPE(pint_env_type), POINTER :: pint_env 176 TYPE(section_vals_type), POINTER :: input 177 TYPE(section_type), POINTER :: input_declaration 178 TYPE(cp_para_env_type), POINTER :: para_env 179 180 CHARACTER(len=*), PARAMETER :: routineN = 'pint_create', routineP = moduleN//':'//routineN 181 182 CHARACTER(len=2*default_string_length) :: msg 183 CHARACTER(len=default_path_length) :: output_file_name, project_name 184 INTEGER :: handle, iat, ibead, icont, idim, idir, & 185 ierr, ig, itmp, nrep, prep, stat 186 LOGICAL :: explicit, ltmp 187 REAL(kind=dp) :: dt, mass, omega 188 REAL(kind=dp), DIMENSION(3, 2) :: seed 189 TYPE(cp_subsys_type), POINTER :: subsys 190 TYPE(f_env_type), POINTER :: f_env 191 TYPE(global_constraint_type), POINTER :: gci 192 TYPE(particle_list_type), POINTER :: particles 193 TYPE(replica_env_type), POINTER :: rep_env 194 TYPE(section_vals_type), POINTER :: constraint_section, gle_section, nose_section, & 195 piglet_section, pile_section, pint_section, qtb_section, transform_section 196 197 CALL timeset(routineN, handle) 198 199 NULLIFY (f_env, subsys, particles, nose_section, gle_section, gci) 200 201 CPASSERT(.NOT. ASSOCIATED(pint_env)) 202 CPASSERT(ASSOCIATED(input)) 203 CPASSERT(input%ref_count > 0) 204 NULLIFY (rep_env) 205 pint_section => section_vals_get_subs_vals(input, "MOTION%PINT") 206 CALL section_vals_val_get(pint_section, "p", i_val=nrep) 207 CALL section_vals_val_get(pint_section, "proc_per_replica", & 208 i_val=prep) 209 ! Maybe let the user have his/her way as long as prep is 210 ! within the bounds of number of CPUs?? 211 IF ((prep < 1) .OR. (prep > para_env%num_pe) .OR. & 212 (MOD(prep*nrep, para_env%num_pe) /= 0)) THEN 213 prep = para_env%num_pe/gcd(para_env%num_pe, nrep) 214 IF (para_env%ionode) THEN 215 WRITE (UNIT=msg, FMT=*) "PINT WARNING: Adjusting number of processors per replica to ", prep 216 CPWARN(msg) 217 END IF 218 END IF 219 220 ! replica_env modifies the global input structure - which is wrong - one 221 ! of the side effects is the inifite adding of the -r-N string to the 222 ! project name and the output file name, which corrupts restart files. 223 ! For now: save the project name and output file name and restore them 224 ! after the rep_env_create has executed - the initialization of the 225 ! replicas will run correctly anyways. 226 ! TODO: modify rep_env so that it behaves better 227 CALL section_vals_val_get(input, "GLOBAL%PROJECT_NAME", c_val=project_name) 228 CALL section_vals_val_get(input, "GLOBAL%OUTPUT_FILE_NAME", c_val=output_file_name) 229 CALL rep_env_create(rep_env, para_env=para_env, input=input, & 230 input_declaration=input_declaration, nrep=nrep, prep=prep, row_force=.TRUE.) 231 CALL section_vals_val_set(input, "GLOBAL%PROJECT_NAME", c_val=TRIM(project_name)) 232 IF (LEN_TRIM(output_file_name) .GT. 0) THEN 233 CALL section_vals_val_set(input, "GLOBAL%OUTPUT_FILE_NAME", c_val=TRIM(output_file_name)) 234 ELSE 235 CALL section_vals_val_unset(input, "GLOBAL%OUTPUT_FILE_NAME") 236 END IF 237 IF (.NOT. ASSOCIATED(rep_env)) RETURN 238 239 ALLOCATE (pint_env) 240 241 NULLIFY (pint_env%logger) 242 pint_env%logger => cp_get_default_logger() 243 CALL cp_add_iter_level(pint_env%logger%iter_info, "PINT") 244 245 last_pint_id = last_pint_id + 1 246 pint_env%id_nr = last_pint_id 247 pint_env%ref_count = 1 248 NULLIFY (pint_env%replicas, pint_env%input, pint_env%staging_env, & 249 pint_env%normalmode_env, pint_env%propagator) 250 pint_env%p = nrep 251 pint_env%replicas => rep_env 252 pint_env%ndim = rep_env%ndim 253 pint_env%input => input 254 255 CALL section_vals_retain(pint_env%input) 256 257 ! get first step, last step, number of steps, etc 258 CALL section_vals_val_get(input, "MOTION%PINT%ITERATION", & 259 i_val=itmp) 260 pint_env%first_step = itmp 261 CALL section_vals_val_get(input, "MOTION%PINT%MAX_STEP", & 262 explicit=explicit) 263 IF (explicit) THEN 264 CALL section_vals_val_get(input, "MOTION%PINT%MAX_STEP", & 265 i_val=itmp) 266 pint_env%last_step = itmp 267 pint_env%num_steps = pint_env%last_step - pint_env%first_step 268 ELSE 269 CALL section_vals_val_get(input, "MOTION%PINT%NUM_STEPS", & 270 i_val=itmp) 271 pint_env%num_steps = itmp 272 pint_env%last_step = pint_env%first_step + pint_env%num_steps 273 END IF 274 275 CALL section_vals_val_get(pint_section, "DT", & 276 r_val=pint_env%dt) 277 pint_env%t = pint_env%first_step*pint_env%dt 278 279 CALL section_vals_val_get(pint_section, "nrespa", i_val=pint_env%nrespa) 280 CALL section_vals_val_get(pint_section, "Temp", r_val=pint_env%kT) 281 CALL section_vals_val_get(pint_section, "T_TOL", & 282 r_val=pint_env%t_tol) 283 284 CALL section_vals_val_get(pint_section, "HARM_INT", i_val=pint_env%harm_integrator) 285 286 ALLOCATE (pint_env%propagator) 287 CALL section_vals_val_get(pint_section, "propagator", & 288 i_val=pint_env%propagator%prop_kind) 289 !Initialize simulation temperature depending on the propagator 290 !As well as the scaling factor for the physical potential 291 IF (pint_env%propagator%prop_kind == propagator_rpmd) THEN 292 pint_env%propagator%temp_phys2sim = REAL(pint_env%p, dp) 293 pint_env%propagator%physpotscale = 1.0_dp 294 ELSE 295 pint_env%propagator%temp_phys2sim = 1.0_dp 296 pint_env%propagator%physpotscale = 1.0_dp/REAL(pint_env%p, dp) 297 END IF 298 pint_env%propagator%temp_sim2phys = 1.0_dp/pint_env%propagator%temp_phys2sim 299 pint_env%kT = pint_env%kT*pint_env%propagator%temp_phys2sim 300 301 CALL section_vals_val_get(pint_section, "transformation", & 302 i_val=pint_env%transform) 303 304 NULLIFY (pint_env%tx, pint_env%tv, pint_env%tv_t, pint_env%tv_old, pint_env%tv_new, pint_env%tf) 305 306 pint_env%nnos = 0 307 pint_env%pimd_thermostat = thermostat_none 308 nose_section => section_vals_get_subs_vals(input, "MOTION%PINT%NOSE") 309 CALL section_vals_get(nose_section, explicit=explicit) 310 IF (explicit) THEN 311 IF (pint_env%propagator%prop_kind == propagator_rpmd) THEN 312 CPABORT("RPMD Propagator with Nose-thermostat not implemented!") 313 END IF 314 CALL section_vals_val_get(nose_section, "nnos", i_val=pint_env%nnos) 315 IF (pint_env%nnos > 0) THEN 316 pint_env%pimd_thermostat = thermostat_nose 317 ALLOCATE ( & 318 pint_env%tx(pint_env%nnos, pint_env%p, pint_env%ndim), & 319 pint_env%tv(pint_env%nnos, pint_env%p, pint_env%ndim), & 320 pint_env%tv_t(pint_env%nnos, pint_env%p, pint_env%ndim), & 321 pint_env%tv_old(pint_env%nnos, pint_env%p, pint_env%ndim), & 322 pint_env%tv_new(pint_env%nnos, pint_env%p, pint_env%ndim), & 323 pint_env%tf(pint_env%nnos, pint_env%p, pint_env%ndim)) 324 pint_env%tx = 0._dp 325 pint_env%tv = 0._dp 326 pint_env%tv_t = 0._dp 327 pint_env%tv_old = 0._dp 328 pint_env%tv_new = 0._dp 329 pint_env%tf = 0._dp 330 END IF 331 END IF 332 333 pint_env%beta = 1._dp/(pint_env%kT*pint_env%propagator%temp_sim2phys) 334!TODO 335! v_tol not in current input structure 336! should also probably be part of nose_section 337! CALL section_vals_val_get(transform_section,"v_tol_nose",r_val=pint_env%v_tol) 338!MK ... but we have to initialise v_tol 339 pint_env%v_tol = 0.0_dp ! to be fixed 340 341 NULLIFY (pint_env%randomG) 342 343 seed(:, :) = next_rng_seed() 344 CALL create_rng_stream(pint_env%randomG, & 345 name="pint_randomG", & 346 distribution_type=GAUSSIAN, & 347 extended_precision=.TRUE., & 348 seed=seed) 349 350 ALLOCATE (pint_env%e_pot_bead(pint_env%p)) 351 pint_env%e_pot_bead = 0._dp 352 pint_env%e_pot_h = 0._dp 353 pint_env%e_kin_beads = 0._dp 354 pint_env%e_pot_t = 0._dp 355 pint_env%e_gle = 0._dp 356 pint_env%e_pile = 0._dp 357 pint_env%e_piglet = 0._dp 358 pint_env%e_qtb = 0._dp 359 pint_env%e_kin_t = 0._dp 360 pint_env%energy(:) = 0.0_dp 361 362!TODO: rearrange to use standard nose hoover chain functions/data types 363 364 ALLOCATE ( & 365 pint_env%x(pint_env%p, pint_env%ndim), & 366 pint_env%v(pint_env%p, pint_env%ndim), & 367 pint_env%f(pint_env%p, pint_env%ndim), & 368 pint_env%external_f(pint_env%p, pint_env%ndim), & 369 pint_env%ux(pint_env%p, pint_env%ndim), & 370 pint_env%ux_t(pint_env%p, pint_env%ndim), & 371 pint_env%uv(pint_env%p, pint_env%ndim), & 372 pint_env%uv_t(pint_env%p, pint_env%ndim), & 373 pint_env%uv_new(pint_env%p, pint_env%ndim), & 374 pint_env%uf(pint_env%p, pint_env%ndim), & 375 pint_env%uf_h(pint_env%p, pint_env%ndim), & 376 pint_env%centroid(pint_env%ndim), & 377 pint_env%rtmp_ndim(pint_env%ndim), & 378 pint_env%rtmp_natom(pint_env%ndim/3), & 379 STAT=stat) 380 CPASSERT(stat == 0) 381 pint_env%x = 0._dp 382 pint_env%v = 0._dp 383 pint_env%f = 0._dp 384 pint_env%external_f = 0._dp 385 pint_env%ux = 0._dp 386 pint_env%ux_t = 0._dp 387 pint_env%uv = 0._dp 388 pint_env%uv_t = 0._dp 389 pint_env%uv_new = 0._dp 390 pint_env%uf = 0._dp 391 pint_env%uf_h = 0._dp 392 pint_env%centroid(:) = 0.0_dp 393 pint_env%rtmp_ndim = 0._dp 394 pint_env%rtmp_natom = 0._dp 395 pint_env%time_per_step = 0.0_dp 396 397 IF (pint_env%transform == transformation_stage) THEN 398 transform_section => section_vals_get_subs_vals(input, & 399 "MOTION%PINT%STAGING") 400 CALL staging_env_create(pint_env%staging_env, transform_section, & 401 p=pint_env%p, kT=pint_env%kT) 402 ELSE 403 transform_section => section_vals_get_subs_vals(input, & 404 "MOTION%PINT%NORMALMODE") 405 CALL normalmode_env_create(pint_env%normalmode_env, & 406 transform_section, p=pint_env%p, kT=pint_env%kT, propagator=pint_env%propagator%prop_kind) 407 IF (para_env%ionode) THEN 408 IF (pint_env%harm_integrator == integrate_numeric) THEN 409 IF (10.0_dp*pint_env%dt/REAL(pint_env%nrespa, dp) > & 410 twopi/(pint_env%p*SQRT(MAXVAL(pint_env%normalmode_env%lambda))* & 411 pint_env%normalmode_env%modefactor)) THEN 412 msg = "PINT WARNING| Number of RESPA steps to small "// & 413 "to integrate the harmonic springs." 414 CPWARN(msg) 415 END IF 416 END IF 417 END IF 418 END IF 419 ALLOCATE (pint_env%mass(pint_env%ndim)) 420 CALL f_env_add_defaults(f_env_id=pint_env%replicas%f_env_id, & 421 f_env=f_env) 422 CALL force_env_get(force_env=f_env%force_env, subsys=subsys) 423 CALL cp_subsys_get(subsys, particles=particles, gci=gci) 424 425!TODO length of pint_env%mass is redundant 426 idim = 0 427 DO iat = 1, pint_env%ndim/3 428 CALL get_atomic_kind(particles%els(iat)%atomic_kind, mass=mass) 429 DO idir = 1, 3 430 idim = idim + 1 431 pint_env%mass(idim) = mass 432 END DO 433 END DO 434 CALL f_env_rm_defaults(f_env, ierr) 435 CPASSERT(ierr == 0) 436 437 ALLOCATE (pint_env%Q(pint_env%p), & 438 pint_env%mass_beads(pint_env%p, pint_env%ndim), & 439 pint_env%mass_fict(pint_env%p, pint_env%ndim)) 440 IF (pint_env%transform == transformation_stage) THEN 441 CALL staging_init_masses(pint_env%staging_env, mass=pint_env%mass, & 442 mass_beads=pint_env%mass_beads, mass_fict=pint_env%mass_fict, & 443 Q=pint_env%Q) 444 ELSE 445 CALL normalmode_init_masses(pint_env%normalmode_env, & 446 mass=pint_env%mass, mass_beads=pint_env%mass_beads, & 447 mass_fict=pint_env%mass_fict, Q=pint_env%Q) 448 END IF 449 450 NULLIFY (pint_env%gle) 451 gle_section => section_vals_get_subs_vals(input, "MOTION%PINT%GLE") 452 CALL section_vals_get(gle_section, explicit=explicit) 453 IF (explicit) THEN 454 ALLOCATE (pint_env%gle) 455 CALL gle_init(pint_env%gle, dt=pint_env%dt/pint_env%nrespa, temp=pint_env%kT, & 456 section=gle_section) 457 IF (pint_env%pimd_thermostat == thermostat_none .AND. pint_env%gle%ndim .GT. 0) THEN 458 pint_env%pimd_thermostat = thermostat_gle 459 460 ! initialize a GLE with ALL degrees of freedom on node 0, 461 ! as it seems to me that here everything but force eval is replicated 462 pint_env%gle%loc_num_gle = pint_env%p*pint_env%ndim 463 pint_env%gle%glob_num_gle = pint_env%gle%loc_num_gle 464 ALLOCATE (pint_env%gle%map_info%index(pint_env%gle%loc_num_gle)) 465 CPASSERT(stat == 0) 466 DO itmp = 1, pint_env%gle%loc_num_gle 467 pint_env%gle%map_info%index(itmp) = itmp 468 ENDDO 469 CALL gle_thermo_create(pint_env%gle, pint_env%gle%loc_num_gle) 470 471 ! here we should have read a_mat and c_mat; 472 !we can therefore compute the matrices needed for the propagator 473 ! deterministic part of the propagator 474 CALL gle_matrix_exp((-pint_env%dt/pint_env%nrespa*0.5_dp)*pint_env%gle%a_mat, & 475 pint_env%gle%ndim, 15, 15, pint_env%gle%gle_t) 476 ! stochastic part 477 CALL gle_cholesky_stab(pint_env%gle%c_mat - MATMUL(pint_env%gle%gle_t, & 478 MATMUL(pint_env%gle%c_mat, TRANSPOSE(pint_env%gle%gle_t))), & 479 pint_env%gle%gle_s, pint_env%gle%ndim) 480 ! and initialize the additional momenta 481 CALL pint_gle_init(pint_env) 482 END IF 483 END IF 484 485 !Setup pile thermostat 486 NULLIFY (pint_env%pile_therm) 487 pile_section => section_vals_get_subs_vals(input, "MOTION%PINT%PILE") 488 CALL section_vals_get(pile_section, explicit=explicit) 489 IF (explicit) THEN 490 CALL cite_reference(Ceriotti2010) 491 CALL section_vals_val_get(pint_env%input, & 492 "MOTION%PINT%INIT%THERMOSTAT_SEED", & 493 i_val=pint_env%thermostat_rng_seed) 494 IF (pint_env%pimd_thermostat == thermostat_none) THEN 495 pint_env%pimd_thermostat = thermostat_pile 496 CALL pint_pile_init(pile_therm=pint_env%pile_therm, & 497 pint_env=pint_env, & 498 normalmode_env=pint_env%normalmode_env, & 499 section=pile_section) 500 ELSE 501 CPABORT("PILE thermostat can't be used with another thermostat.") 502 END IF 503 END IF 504 505 !Setup PIGLET thermostat 506 NULLIFY (pint_env%piglet_therm) 507 piglet_section => section_vals_get_subs_vals(input, "MOTION%PINT%PIGLET") 508 CALL section_vals_get(piglet_section, explicit=explicit) 509 IF (explicit) THEN 510 CALL cite_reference(Ceriotti2012) 511 CALL section_vals_val_get(pint_env%input, & 512 "MOTION%PINT%INIT%THERMOSTAT_SEED", & 513 i_val=pint_env%thermostat_rng_seed) 514 IF (pint_env%pimd_thermostat == thermostat_none) THEN 515 pint_env%pimd_thermostat = thermostat_piglet 516 CALL pint_piglet_create(pint_env%piglet_therm, & 517 pint_env, & 518 piglet_section) 519 CALL pint_piglet_init(pint_env%piglet_therm, & 520 pint_env, & 521 piglet_section, & 522 dt=pint_env%dt, para_env=para_env) 523 ELSE 524 CPABORT("PILE thermostat can't be used with another thermostat.") 525 END IF 526 END IF 527 528 !Setup qtb thermostat 529 NULLIFY (pint_env%qtb_therm) 530 qtb_section => section_vals_get_subs_vals(input, "MOTION%PINT%QTB") 531 CALL section_vals_get(qtb_section, explicit=explicit) 532 IF (explicit) THEN 533 CALL cite_reference(Brieuc2016) 534 CALL section_vals_val_get(pint_env%input, & 535 "MOTION%PINT%INIT%THERMOSTAT_SEED", & 536 i_val=pint_env%thermostat_rng_seed) 537 IF (pint_env%pimd_thermostat == thermostat_none) THEN 538 pint_env%pimd_thermostat = thermostat_qtb 539 CALL pint_qtb_init(qtb_therm=pint_env%qtb_therm, & 540 pint_env=pint_env, & 541 normalmode_env=pint_env%normalmode_env, & 542 section=qtb_section) 543 ELSE 544 CPABORT("QTB thermostat can't be used with another thermostat.") 545 END IF 546 END IF 547 548 !Initialize integrator scheme 549 CALL section_vals_val_get(pint_section, "HARM_INT", i_val=pint_env%harm_integrator) 550 IF (pint_env%harm_integrator == integrate_exact) THEN 551 IF (pint_env%pimd_thermostat == thermostat_nose) THEN 552 WRITE (UNIT=msg, FMT=*) "PINT WARNING| Nose Thermostat only avaliable in "// & 553 "the numeric harmonic integrator. Switching to numeric harmonic integrator." 554 CPWARN(msg) 555 pint_env%harm_integrator = integrate_numeric 556 END IF 557 IF (pint_env%pimd_thermostat == thermostat_gle) THEN 558 WRITE (UNIT=msg, FMT=*) "PINT WARNING| GLE Thermostat only avaliable in "// & 559 "the numeric harmonic integrator. Switching to numeric harmonic integrator." 560 CPWARN(msg) 561 pint_env%harm_integrator = integrate_numeric 562 END IF 563 ELSE IF (pint_env%harm_integrator == integrate_numeric) THEN 564 IF (pint_env%pimd_thermostat == thermostat_pile) THEN 565 WRITE (UNIT=msg, FMT=*) "PINT WARNING| PILE Thermostat only avaliable in "// & 566 "the exact harmonic integrator. Switching to exact harmonic integrator." 567 CPWARN(msg) 568 pint_env%harm_integrator = integrate_exact 569 END IF 570 IF (pint_env%pimd_thermostat == thermostat_piglet) THEN 571 WRITE (UNIT=msg, FMT=*) "PINT WARNING| PIGLET Thermostat only avaliable in "// & 572 "the exact harmonic integrator. Switching to exact harmonic integrator." 573 CPWARN(msg) 574 pint_env%harm_integrator = integrate_exact 575 END IF 576 IF (pint_env%pimd_thermostat == thermostat_qtb) THEN 577 WRITE (UNIT=msg, FMT=*) "PINT WARNING| QTB Thermostat only avaliable in "// & 578 "the exact harmonic integrator. Switching to exact harmonic integrator." 579 CPWARN(msg) 580 pint_env%harm_integrator = integrate_exact 581 END IF 582 END IF 583 584 IF (pint_env%harm_integrator == integrate_exact) THEN 585 IF (pint_env%nrespa /= 1) THEN 586 pint_env%nrespa = 1 587 WRITE (UNIT=msg, FMT=*) "PINT WARNING| Adjusting NRESPA to 1 for exact harmonic integration." 588 CPWARN(msg) 589 END IF 590 NULLIFY (pint_env%wsinex) 591 ALLOCATE (pint_env%wsinex(pint_env%p)) 592 NULLIFY (pint_env%iwsinex) 593 ALLOCATE (pint_env%iwsinex(pint_env%p)) 594 NULLIFY (pint_env%cosex) 595 ALLOCATE (pint_env%cosex(pint_env%p)) 596 dt = pint_env%dt/REAL(pint_env%nrespa, KIND=dp) 597 !Centroid mode shoud not be propagated 598 pint_env%wsinex(1) = 0.0_dp 599 pint_env%iwsinex(1) = dt 600 pint_env%cosex(1) = 1.0_dp 601 DO ibead = 2, pint_env%p 602 omega = SQRT(pint_env%normalmode_env%lambda(ibead)) 603 pint_env%wsinex(ibead) = SIN(omega*dt)*omega 604 pint_env%iwsinex(ibead) = SIN(omega*dt)/omega 605 pint_env%cosex(ibead) = COS(omega*dt) 606 END DO 607 END IF 608 609 CALL section_vals_val_get(pint_section, "FIX_CENTROID_POS", & 610 l_val=ltmp) 611 IF (ltmp .AND. (pint_env%transform .EQ. transformation_normal)) THEN 612 pint_env%first_propagated_mode = 2 613 ELSE 614 pint_env%first_propagated_mode = 1 615 END IF 616 617 ! Constraint information: 618 NULLIFY (pint_env%simpar) 619 constraint_section => section_vals_get_subs_vals(pint_env%input, & 620 "MOTION%CONSTRAINT") 621 CALL section_vals_get(constraint_section, explicit=explicit) 622 CALL create_simpar_type(pint_env%simpar) 623 pint_env%simpar%constraint = explicit 624 pint_env%kTcorr = 1.0_dp 625 IF (explicit) THEN 626 ! Staging not supported yet, since lowest mode is assumed 627 ! to be related to centroid 628 IF (pint_env%transform == transformation_stage) THEN 629 CPABORT("Centroid constraints are not supported for staging transformation") 630 END IF 631 632 ! Check thermostats that are not supported: 633 IF (pint_env%pimd_thermostat == thermostat_gle) THEN 634 WRITE (UNIT=msg, FMT=*) "GLE Thermostat not supported for "// & 635 "centroid constraints. Switch to NOSE for numeric integration." 636 CPABORT(msg) 637 END IF 638 ! Warn for NOSE 639 IF (pint_env%pimd_thermostat == thermostat_nose) THEN 640 WRITE (UNIT=msg, FMT=*) "PINT WARNING| Nose Thermostat set to "// & 641 "zero for constrained atoms. Careful interpretation of temperature." 642 CPWARN(msg) 643 WRITE (UNIT=msg, FMT=*) "PINT WARNING| Lagrange multipliers are "// & 644 "are printed every RESPA step and need to be treated carefully." 645 CPWARN(msg) 646 END IF 647 648 CALL section_vals_val_get(constraint_section, "SHAKE_TOLERANCE", & 649 r_val=pint_env%simpar%shake_tol) 650 pint_env%simpar%info_constraint = cp_print_key_unit_nr(pint_env%logger, & 651 constraint_section, & 652 "CONSTRAINT_INFO", & 653 extension=".shakeLog", & 654 log_filename=.FALSE.) 655 pint_env%simpar%lagrange_multipliers = cp_print_key_unit_nr(pint_env%logger, & 656 constraint_section, & 657 "LAGRANGE_MULTIPLIERS", & 658 extension=".LagrangeMultLog", & 659 log_filename=.FALSE.) 660 pint_env%simpar%dump_lm = & 661 BTEST(cp_print_key_should_output(pint_env%logger%iter_info, & 662 constraint_section, & 663 "LAGRANGE_MULTIPLIERS"), cp_p_file) 664 665 ! Determine constrained atoms: 666 pint_env%n_atoms_constraints = 0 667 DO ig = 1, gci%ncolv%ntot 668 ! Double counts, if the same atom is involved in different collective variables 669 pint_env%n_atoms_constraints = pint_env%n_atoms_constraints + SIZE(gci%colv_list(ig)%i_atoms) 670 END DO 671 672 ALLOCATE (pint_env%atoms_constraints(pint_env%n_atoms_constraints)) 673 icont = 0 674 DO ig = 1, gci%ncolv%ntot 675 DO iat = 1, SIZE(gci%colv_list(ig)%i_atoms) 676 icont = icont + 1 677 pint_env%atoms_constraints(icont) = gci%colv_list(ig)%i_atoms(iat) 678 END DO 679 END DO 680 681 ! Set the correction to the temperature due to the frozen degrees of freedom in NOSE: 682 CALL section_vals_val_get(pint_section, "kT_CORRECTION", & 683 l_val=ltmp) 684 IF (ltmp) THEN 685 pint_env%kTcorr = 1.0_dp + REAL(3*pint_env%n_atoms_constraints, dp)/(REAL(pint_env%ndim, dp)*REAL(pint_env%p, dp)) 686 END IF 687 END IF 688 689 CALL timestop(handle) 690 691 RETURN 692 END SUBROUTINE pint_create 693 694! *************************************************************************** 695!> \brief Retain a path integral environment 696!> \param pint_env the pint_env to retain 697!> \author Fawzi Mohamed 698! ************************************************************************************************** 699 SUBROUTINE pint_retain(pint_env) 700 TYPE(pint_env_type), POINTER :: pint_env 701 702 CHARACTER(len=*), PARAMETER :: routineN = 'pint_retain', routineP = moduleN//':'//routineN 703 704 CPASSERT(ASSOCIATED(pint_env)) 705 CPASSERT(pint_env%ref_count > 0) 706 pint_env%ref_count = pint_env%ref_count + 1 707 RETURN 708 END SUBROUTINE pint_retain 709 710! *************************************************************************** 711!> \brief Release a path integral environment 712!> \param pint_env the pint_env to release 713!> \par History 714!> Added normal mode transformation [hforbert] 715!> \author Fawzi Mohamed 716! ************************************************************************************************** 717 SUBROUTINE pint_release(pint_env) 718 TYPE(pint_env_type), POINTER :: pint_env 719 720 CHARACTER(len=*), PARAMETER :: routineN = 'pint_release', routineP = moduleN//':'//routineN 721 722 IF (ASSOCIATED(pint_env)) THEN 723 CPASSERT(pint_env%ref_count > 0) 724 pint_env%ref_count = pint_env%ref_count - 1 725 IF (pint_env%ref_count == 0) THEN 726 CALL rep_env_release(pint_env%replicas) 727 CALL section_vals_release(pint_env%input) 728 IF (ASSOCIATED(pint_env%staging_env)) THEN 729 CALL staging_release(pint_env%staging_env) 730 END IF 731 IF (ASSOCIATED(pint_env%normalmode_env)) THEN 732 CALL normalmode_release(pint_env%normalmode_env) 733 END IF 734 CALL delete_rng_stream(pint_env%randomG) 735 736 DEALLOCATE (pint_env%mass) 737 DEALLOCATE (pint_env%e_pot_bead) 738 739 DEALLOCATE (pint_env%x) 740 DEALLOCATE (pint_env%v) 741 DEALLOCATE (pint_env%f) 742 DEALLOCATE (pint_env%external_f) 743 DEALLOCATE (pint_env%mass_beads) 744 DEALLOCATE (pint_env%mass_fict) 745 DEALLOCATE (pint_env%ux) 746 DEALLOCATE (pint_env%ux_t) 747 DEALLOCATE (pint_env%uv) 748 DEALLOCATE (pint_env%uv_t) 749 DEALLOCATE (pint_env%uv_new) 750 DEALLOCATE (pint_env%uf) 751 DEALLOCATE (pint_env%uf_h) 752 DEALLOCATE (pint_env%centroid) 753 DEALLOCATE (pint_env%rtmp_ndim) 754 DEALLOCATE (pint_env%rtmp_natom) 755 DEALLOCATE (pint_env%propagator) 756 757 IF (pint_env%simpar%constraint) THEN 758 DEALLOCATE (pint_env%atoms_constraints) 759 END IF 760 CALL release_simpar_type(pint_env%simpar) 761 762 IF (pint_env%harm_integrator == integrate_exact) THEN 763 DEALLOCATE (pint_env%wsinex) 764 DEALLOCATE (pint_env%iwsinex) 765 DEALLOCATE (pint_env%cosex) 766 END IF 767 768 SELECT CASE (pint_env%pimd_thermostat) 769 CASE (thermostat_nose) 770 DEALLOCATE (pint_env%tx) 771 DEALLOCATE (pint_env%tv) 772 DEALLOCATE (pint_env%tv_t) 773 DEALLOCATE (pint_env%tv_old) 774 DEALLOCATE (pint_env%tv_new) 775 DEALLOCATE (pint_env%tf) 776 CASE (thermostat_gle) 777 CALL gle_dealloc(pint_env%gle) 778 CASE (thermostat_pile) 779 CALL pint_pile_release(pint_env%pile_therm) 780 CASE (thermostat_piglet) 781 CALL pint_piglet_release(pint_env%piglet_therm) 782 CASE (thermostat_qtb) 783 CALL pint_qtb_release(pint_env%qtb_therm) 784 END SELECT 785 786 DEALLOCATE (pint_env%Q) 787 788 DEALLOCATE (pint_env) 789 790 END IF 791 END IF 792 793 NULLIFY (pint_env) 794 795 RETURN 796 END SUBROUTINE pint_release 797 798! *************************************************************************** 799!> \brief Tests the path integral methods 800!> \param para_env parallel environment 801!> \param input the input to test 802!> \param input_declaration ... 803!> \author fawzi 804! ************************************************************************************************** 805 SUBROUTINE pint_test(para_env, input, input_declaration) 806 TYPE(cp_para_env_type), POINTER :: para_env 807 TYPE(section_vals_type), POINTER :: input 808 TYPE(section_type), POINTER :: input_declaration 809 810 CHARACTER(len=*), PARAMETER :: routineN = 'pint_test', routineP = moduleN//':'//routineN 811 812 INTEGER :: i, ib, idim, unit_nr 813 REAL(kind=dp) :: c, e_h, err 814 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: x1 815 TYPE(pint_env_type), POINTER :: pint_env 816 817 CPASSERT(ASSOCIATED(para_env)) 818 CPASSERT(ASSOCIATED(input)) 819 CPASSERT(para_env%ref_count > 0) 820 CPASSERT(input%ref_count > 0) 821 NULLIFY (pint_env) 822 unit_nr = cp_logger_get_default_io_unit() 823 CALL pint_create(pint_env, input, input_declaration, para_env) 824 IF (ASSOCIATED(pint_env)) THEN 825 ALLOCATE (x1(pint_env%ndim, pint_env%p)) 826 x1(:, :) = pint_env%x 827 CALL pint_x2u(pint_env) 828 pint_env%x = 0._dp 829 CALL pint_u2x(pint_env) 830 err = 0._dp 831 DO i = 1, pint_env%ndim 832 err = MAX(err, ABS(x1(1, i) - pint_env%x(1, i))) 833 END DO 834 IF (unit_nr > 0) WRITE (unit_nr, *) "diff_r1="//cp_to_string(err) 835 836 CALL pint_calc_uf_h(pint_env, e_h=e_h) 837 c = -pint_env%staging_env%w_p**2 838 pint_env%f = 0._dp 839 DO idim = 1, pint_env%ndim 840 DO ib = 1, pint_env%p 841 pint_env%f(ib, idim) = pint_env%f(ib, idim) + & 842 c*(2._dp*pint_env%x(ib, idim) & 843 - pint_env%x(MODULO(ib - 2, pint_env%p) + 1, idim) & 844 - pint_env%x(MODULO(ib, pint_env%p) + 1, idim)) 845 END DO 846 END DO 847 CALL pint_f2uf(pint_env) 848 err = 0._dp 849 DO idim = 1, pint_env%ndim 850 DO ib = 1, pint_env%p 851 err = MAX(err, ABS(pint_env%uf(ib, idim) - pint_env%uf_h(ib, idim))) 852 END DO 853 END DO 854 IF (unit_nr > 0) WRITE (unit_nr, *) "diff_f_h="//cp_to_string(err) 855 END IF 856 RETURN 857 END SUBROUTINE pint_test 858 859! *************************************************************************** 860!> \brief Perform a path integral simulation 861!> \param para_env parallel environment 862!> \param input the input to test 863!> \param input_declaration ... 864!> \param globenv ... 865!> \par History 866!> 2003-11 created [fawzi] 867!> 2009-12-14 globenv parameter added to handle soft exit 868!> requests [lwalewski] 869!> 2016-07-14 Modified to work with independent helium_env [cschran] 870!> \author Fawzi Mohamed 871! ************************************************************************************************** 872 SUBROUTINE do_pint_run(para_env, input, input_declaration, globenv) 873 TYPE(cp_para_env_type), POINTER :: para_env 874 TYPE(section_vals_type), POINTER :: input 875 TYPE(section_type), POINTER :: input_declaration 876 TYPE(global_environment_type), POINTER :: globenv 877 878 CHARACTER(len=*), PARAMETER :: routineN = 'do_pint_run', routineP = moduleN//':'//routineN 879 INTEGER, PARAMETER :: helium_only_mid = 1, & 880 int_pot_scan_mid = 4, & 881 solute_only_mid = 2, & 882 solute_with_helium_mid = 3 883 884 CHARACTER(len=default_string_length) :: stmp 885 INTEGER :: handle, mode 886 LOGICAL :: explicit, helium_only, int_pot_scan, & 887 solvent_present 888 TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env 889 TYPE(pint_env_type), POINTER :: pint_env 890 TYPE(section_vals_type), POINTER :: helium_section 891 892 CALL timeset(routineN, handle) 893 894 CPASSERT(ASSOCIATED(para_env)) 895 CPASSERT(ASSOCIATED(input)) 896 CPASSERT(para_env%ref_count > 0) 897 CPASSERT(input%ref_count > 0) 898 899 ! check if helium solvent is present 900 NULLIFY (helium_section) 901 helium_section => section_vals_get_subs_vals(input, & 902 "MOTION%PINT%HELIUM") 903 CALL section_vals_get(helium_section, explicit=explicit) 904 IF (explicit) THEN 905 CALL section_vals_val_get(helium_section, "_SECTION_PARAMETERS_", & 906 l_val=solvent_present) 907 ELSE 908 solvent_present = .FALSE. 909 END IF 910 911 ! check if there is anything but helium 912 IF (solvent_present) THEN 913 CALL section_vals_val_get(helium_section, "HELIUM_ONLY", & 914 l_val=helium_only) 915 ELSE 916 helium_only = .FALSE. 917 END IF 918 919 ! check wheather to perform solute-helium interaction pot scan 920 IF (solvent_present) THEN 921 CALL section_vals_val_get(helium_section, "INTERACTION_POT_SCAN", & 922 l_val=int_pot_scan) 923 ELSE 924 int_pot_scan = .FALSE. 925 END IF 926 927 ! input consistency check 928 IF (helium_only .AND. int_pot_scan) THEN 929 stmp = "Options HELIUM_ONLY and INTERACTION_POT_SCAN are exclusive" 930 CPABORT(stmp) 931 END IF 932 933 ! select mode of operation 934 mode = 0 935 IF (solvent_present) THEN 936 IF (helium_only) THEN 937 mode = helium_only_mid 938 ELSE 939 IF (int_pot_scan) THEN 940 mode = int_pot_scan_mid 941 ELSE 942 mode = solute_with_helium_mid 943 END IF 944 END IF 945 ELSE 946 mode = solute_only_mid 947 END IF 948 949 NULLIFY (pint_env) 950 951 ! perform the simulation according to the chosen mode 952 SELECT CASE (mode) 953 954 CASE (helium_only_mid) 955 CALL helium_create(helium_env, input) 956 CALL helium_init(helium_env, pint_env) 957 CALL helium_do_run(helium_env, globenv) 958 CALL helium_release(helium_env) 959 960 CASE (solute_only_mid) 961 CALL pint_create(pint_env, input, input_declaration, para_env) 962 CALL pint_init(pint_env) 963 CALL pint_do_run(pint_env, globenv) 964 CALL pint_release(pint_env) 965 966 CASE (int_pot_scan_mid) 967 CALL pint_create(pint_env, input, input_declaration, para_env) 968! TODO only initialization of positions is necessary, but rep_env_calc_e_f called 969! from within pint_init_f does something to the replica environments which can not be 970! avoided (has something to do with f_env_add_defaults) so leaving for now.. 971 CALL pint_init(pint_env) 972 CALL helium_create(helium_env, input, solute=pint_env) 973 CALL pint_run_scan(pint_env, helium_env) 974 CALL helium_release(helium_env) 975 CALL pint_release(pint_env) 976 977 CASE (solute_with_helium_mid) 978 CALL pint_create(pint_env, input, input_declaration, para_env) 979 ! init pint wihtout helium forces (they are not yet initialized) 980 CALL pint_init(pint_env) 981 ! init helium with solute's positions (they are already initialized) 982 CALL helium_create(helium_env, input, solute=pint_env) 983 CALL helium_init(helium_env, pint_env) 984 ! reinit pint forces with helium forces (they are now initialized) 985 CALL pint_init_f(pint_env, helium_env=helium_env) 986 987 CALL pint_do_run(pint_env, globenv, helium_env=helium_env) 988 CALL helium_release(helium_env) 989 CALL pint_release(pint_env) 990 991 CASE DEFAULT 992 CPABORT("Unknown mode ("//TRIM(ADJUSTL(cp_to_string(mode)))//")") 993 END SELECT 994 995 CALL timestop(handle) 996 997 RETURN 998 END SUBROUTINE do_pint_run 999 1000! *************************************************************************** 1001!> \brief Reads the restart, initializes the beads, etc. 1002!> \param pint_env ... 1003!> \par History 1004!> 11.2003 created [fawzi] 1005!> actually ASSIGN input pointer [hforbert] 1006!> 2010-12-16 turned into a wrapper routine [lwalewski] 1007!> \author Fawzi Mohamed 1008! ************************************************************************************************** 1009 SUBROUTINE pint_init(pint_env) 1010 1011 TYPE(pint_env_type), POINTER :: pint_env 1012 1013 CHARACTER(len=*), PARAMETER :: routineN = 'pint_init', routineP = moduleN//':'//routineN 1014 1015 CPASSERT(ASSOCIATED(pint_env)) 1016 CPASSERT(pint_env%ref_count > 0) 1017 1018 CALL pint_init_x(pint_env) 1019 CALL pint_init_v(pint_env) 1020 CALL pint_init_t(pint_env) 1021 CALL pint_init_f(pint_env) 1022 1023 RETURN 1024 END SUBROUTINE pint_init 1025 1026! *************************************************************************** 1027!> \brief Assign initial postions to the beads. 1028!> \param pint_env ... 1029!> \date 2010-12-15 1030!> \author Lukasz Walewski 1031!> \note Initialization is done in the following way: 1032!> 1. assign all beads with the same classical positions from 1033!> FORCE_EVAL (hot start) 1034!> 2. spread the beads around classical positions as if they were 1035!> free particles (if requested) 1036!> 3. replace positions generated in steps 1-2 with the explicit 1037!> ones if they are explicitly given in the input structure 1038!> 4. apply Gaussian noise to the positions generated so far (if 1039!> requested) 1040! ************************************************************************************************** 1041 SUBROUTINE pint_init_x(pint_env) 1042 1043 TYPE(pint_env_type), POINTER :: pint_env 1044 1045 CHARACTER(len=*), PARAMETER :: routineN = 'pint_init_x', routineP = moduleN//':'//routineN 1046 1047 CHARACTER(len=default_string_length) :: msg, tmp 1048 INTEGER :: ia, ib, ic, idim, input_seed, n_rep_val 1049 LOGICAL :: done_init, done_levy, done_rand, & 1050 explicit, levycorr, ltmp 1051 REAL(kind=dp) :: tcorr, var 1052 REAL(kind=dp), DIMENSION(3) :: x0 1053 REAL(kind=dp), DIMENSION(3, 2) :: seed 1054 REAL(kind=dp), DIMENSION(:), POINTER :: bx, r_vals 1055 TYPE(rng_stream_type), POINTER :: rng_gaussian 1056 TYPE(section_vals_type), POINTER :: input_section 1057 1058 CPASSERT(ASSOCIATED(pint_env)) 1059 CPASSERT(pint_env%ref_count > 0) 1060 1061 DO idim = 1, pint_env%ndim 1062 DO ib = 1, pint_env%p 1063 pint_env%x(ib, idim) = pint_env%replicas%r(idim, ib) 1064 END DO 1065 END DO 1066 1067 done_levy = .FALSE. 1068 CALL section_vals_val_get(pint_env%input, & 1069 "MOTION%PINT%INIT%LEVY_POS_SAMPLE", & 1070 l_val=ltmp) 1071 CALL section_vals_val_get(pint_env%input, & 1072 "MOTION%PINT%INIT%LEVY_TEMP_FACTOR", & 1073 r_val=tcorr) 1074 IF (ltmp) THEN 1075 1076 NULLIFY (bx) 1077 ALLOCATE (bx(3*pint_env%p)) 1078 NULLIFY (rng_gaussian) 1079 CALL section_vals_val_get(pint_env%input, & 1080 "MOTION%PINT%INIT%LEVY_SEED", i_val=input_seed) 1081 seed(:, :) = REAL(input_seed, KIND=dp) 1082! seed(:,:) = next_rng_seed() 1083 CALL create_rng_stream(rng_gaussian, & 1084 name="tmp_rng_gaussian", & 1085 distribution_type=GAUSSIAN, & 1086 extended_precision=.TRUE., & 1087 seed=seed) 1088 1089 CALL section_vals_val_get(pint_env%input, & 1090 "MOTION%PINT%INIT%LEVY_CORRELATED", & 1091 l_val=levycorr) 1092 1093 IF (levycorr) THEN 1094 1095 ! correlated Levy walk - the same path for all atoms 1096 x0 = (/0.0_dp, 0.0_dp, 0.0_dp/) 1097 CALL pint_levy_walk(x0, pint_env%p, 1.0_dp, bx, rng_gaussian) 1098 idim = 0 1099 DO ia = 1, pint_env%ndim/3 1100 var = SQRT(1.0_dp/(pint_env%kT*tcorr*pint_env%mass(3*ia))) 1101 DO ic = 1, 3 1102 idim = idim + 1 1103 DO ib = 1, pint_env%p 1104 pint_env%x(ib, idim) = pint_env%x(ib, idim) + bx(3*(ib - 1) + ic)*var 1105 END DO 1106 END DO 1107 END DO 1108 1109 ELSE 1110 1111 ! uncorrelated bead initialization - distinct Levy walk for each atom 1112 idim = 0 1113 DO ia = 1, pint_env%ndim/3 1114 x0(1) = pint_env%x(1, 3*(ia - 1) + 1) 1115 x0(2) = pint_env%x(1, 3*(ia - 1) + 2) 1116 x0(3) = pint_env%x(1, 3*(ia - 1) + 3) 1117 var = SQRT(1.0_dp/(pint_env%kT*tcorr*pint_env%mass(3*ia))) 1118 CALL pint_levy_walk(x0, pint_env%p, var, bx, rng_gaussian) 1119 DO ic = 1, 3 1120 idim = idim + 1 1121 DO ib = 1, pint_env%p 1122 pint_env%x(ib, idim) = pint_env%x(ib, idim) + bx(3*(ib - 1) + ic) 1123 END DO 1124 END DO 1125 END DO 1126 1127 END IF 1128 1129 CALL delete_rng_stream(rng_gaussian) 1130 DEALLOCATE (bx) 1131 done_levy = .TRUE. 1132 END IF 1133 1134 done_init = .FALSE. 1135 NULLIFY (input_section) 1136 input_section => section_vals_get_subs_vals(pint_env%input, & 1137 "MOTION%PINT%BEADS%COORD") 1138 CALL section_vals_get(input_section, explicit=explicit) 1139 IF (explicit) THEN 1140 CALL section_vals_val_get(input_section, "_DEFAULT_KEYWORD_", & 1141 n_rep_val=n_rep_val) 1142 IF (n_rep_val > 0) THEN 1143 CPASSERT(n_rep_val == 1) 1144 CALL section_vals_val_get(input_section, "_DEFAULT_KEYWORD_", & 1145 r_vals=r_vals) 1146 IF (SIZE(r_vals) /= pint_env%p*pint_env%ndim) & 1147 CPABORT("Invalid size of MOTION%PINT%BEADS%COORD") 1148 ic = 0 1149 DO idim = 1, pint_env%ndim 1150 DO ib = 1, pint_env%p 1151 ic = ic + 1 1152 pint_env%x(ib, idim) = r_vals(ic) 1153 END DO 1154 END DO 1155 done_init = .TRUE. 1156 END IF 1157 END IF 1158 1159 done_rand = .FALSE. 1160 CALL section_vals_val_get(pint_env%input, & 1161 "MOTION%PINT%INIT%RANDOMIZE_POS", & 1162 l_val=ltmp) 1163 IF (ltmp) THEN 1164 DO idim = 1, pint_env%ndim 1165 DO ib = 1, pint_env%p 1166 pint_env%x(ib, idim) = pint_env%x(ib, idim) + & 1167 next_random_number(rng_stream=pint_env%randomG, & 1168 variance=pint_env%beta/ & 1169 SQRT(12.0_dp*pint_env%mass(idim))) 1170 END DO 1171 END DO 1172 done_rand = .TRUE. 1173 END IF 1174 1175 WRITE (tmp, '(A)') "Bead positions initialization:" 1176 IF (done_init) THEN 1177 WRITE (msg, '(A,A)') TRIM(tmp), " input structure" 1178 ELSE IF (done_levy) THEN 1179 WRITE (msg, '(A,A)') TRIM(tmp), " Levy random walk" 1180 ELSE 1181 WRITE (msg, '(A,A)') TRIM(tmp), " hot start" 1182 END IF 1183 CALL pint_write_line(msg) 1184 1185 IF (done_levy) THEN 1186 WRITE (msg, '(A,F6.3)') "Levy walk at effective temperature: ", tcorr 1187 END IF 1188 1189 IF (done_rand) THEN 1190 WRITE (msg, '(A)') "Added gaussian noise to the positions of the beads." 1191 CALL pint_write_line(msg) 1192 END IF 1193 1194 RETURN 1195 END SUBROUTINE pint_init_x 1196 1197! *************************************************************************** 1198!> \brief Initialize velocities 1199!> \param pint_env the pint env in which you should initialize the 1200!> velocity 1201!> \par History 1202!> 2010-12-16 gathered all velocity-init code here [lwalewski] 1203!> 2011-04-05 added centroid velocity initialization [lwalewski] 1204!> 2011-12-19 removed optional parameter kT, target temperature is 1205!> now determined from the input directly [lwalewski] 1206!> \author fawzi 1207!> \note Initialization is done according to the following protocol: 1208!> 1. set all the velocities to FORCE_EVAL%SUBSYS%VELOCITY if present 1209!> 2. scale the velocities according to the actual temperature 1210!> (has no effect if vels not present in 1.) 1211!> 3. draw vels for the remaining dof from MB distribution 1212!> (all or non-centroid modes only depending on 1.) 1213!> 4. add random noise to the centroid vels if CENTROID_SPEED == T 1214!> 5. set the vels for all dof to 0.0 if VELOCITY_QUENCH == T 1215!> 6. set the vels according to the explicit values from the input 1216!> if present 1217! ************************************************************************************************** 1218 SUBROUTINE pint_init_v(pint_env) 1219 TYPE(pint_env_type), POINTER :: pint_env 1220 1221 CHARACTER(len=*), PARAMETER :: routineN = 'pint_init_v', routineP = moduleN//':'//routineN 1222 1223 CHARACTER(len=default_string_length) :: msg, stmp, stmp1, stmp2, unit_str 1224 INTEGER :: first_mode, i, ia, ib, ic, idim, ierr, & 1225 itmp, j, n_rep_val, nparticle, & 1226 nparticle_kind 1227 LOGICAL :: done_init, done_quench, done_scale, & 1228 done_sped, explicit, ltmp, vels_present 1229 REAL(kind=dp) :: actual_t, ek, factor, rtmp, target_t, & 1230 unit_conv 1231 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: vel 1232 REAL(kind=dp), DIMENSION(:), POINTER :: r_vals 1233 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds 1234 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 1235 TYPE(cell_type), POINTER :: cell 1236 TYPE(cp_logger_type), POINTER :: logger 1237 TYPE(cp_subsys_type), POINTER :: subsys 1238 TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles 1239 TYPE(f_env_type), POINTER :: f_env 1240 TYPE(global_constraint_type), POINTER :: gci 1241 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds 1242 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set 1243 TYPE(molecule_list_type), POINTER :: molecules 1244 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 1245 TYPE(particle_list_type), POINTER :: particles 1246 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 1247 TYPE(section_vals_type), POINTER :: input_section 1248 1249 CPASSERT(ASSOCIATED(pint_env)) 1250 CPASSERT(pint_env%ref_count > 0) 1251 1252 NULLIFY (logger) 1253 logger => cp_get_default_logger() 1254 1255 ! Get centroid constraint info, if needed 1256 ! Create a force environment which will be identical to 1257 ! the bead that is being processed by the processor. 1258 IF (pint_env%simpar%constraint) THEN 1259 NULLIFY (subsys, cell) 1260 NULLIFY (atomic_kinds, local_particles, particles) 1261 NULLIFY (local_molecules, molecules, molecule_kinds, gci) 1262 NULLIFY (atomic_kind_set, molecule_kind_set, particle_set, molecule_set) 1263 1264 CALL f_env_add_defaults(f_env_id=pint_env%replicas%f_env_id, f_env=f_env) 1265 CALL force_env_get(force_env=f_env%force_env, subsys=subsys) 1266 CALL f_env_rm_defaults(f_env, ierr) 1267 CPASSERT(ierr == 0) 1268 1269 ! Get gci and more from subsys 1270 CALL cp_subsys_get(subsys=subsys, & 1271 cell=cell, & 1272 atomic_kinds=atomic_kinds, & 1273 local_particles=local_particles, & 1274 particles=particles, & 1275 local_molecules=local_molecules, & 1276 molecules=molecules, & 1277 molecule_kinds=molecule_kinds, & 1278 gci=gci) 1279 1280 nparticle_kind = atomic_kinds%n_els 1281 atomic_kind_set => atomic_kinds%els 1282 molecule_kind_set => molecule_kinds%els 1283 nparticle = particles%n_els 1284 particle_set => particles%els 1285 molecule_set => molecules%els 1286 1287 ! Allocate work storage 1288 ALLOCATE (vel(3, nparticle)) 1289 vel(:, :) = 0.0_dp 1290 CALL getold(gci, local_molecules, molecule_set, & 1291 molecule_kind_set, particle_set, cell) 1292 END IF 1293 1294 ! read the velocities from the input file if they are given explicitly 1295 vels_present = .FALSE. 1296 NULLIFY (input_section) 1297 input_section => section_vals_get_subs_vals(pint_env%input, & 1298 "FORCE_EVAL%SUBSYS%VELOCITY") 1299 CALL section_vals_get(input_section, explicit=explicit) 1300 IF (explicit) THEN 1301 1302 CALL section_vals_val_get(input_section, "PINT_UNIT", & 1303 c_val=unit_str) 1304 unit_conv = cp_unit_to_cp2k(1.0_dp, TRIM(unit_str)) 1305 1306 ! assign all the beads with the same velocities from FORCE_EVAL%SUBSYS%VELOCITY 1307 NULLIFY (r_vals) 1308 CALL section_vals_val_get(input_section, "_DEFAULT_KEYWORD_", & 1309 n_rep_val=n_rep_val) 1310 stmp = "" 1311 WRITE (stmp, *) n_rep_val 1312 msg = "Invalid number of atoms in FORCE_EVAL%SUBSYS%VELOCITY ("// & 1313 TRIM(ADJUSTL(stmp))//")." 1314 IF (3*n_rep_val /= pint_env%ndim) & 1315 CPABORT(msg) 1316 DO ia = 1, pint_env%ndim/3 1317 CALL section_vals_val_get(input_section, "_DEFAULT_KEYWORD_", & 1318 i_rep_val=ia, r_vals=r_vals) 1319 itmp = SIZE(r_vals) 1320 stmp = "" 1321 WRITE (stmp, *) itmp 1322 msg = "Number of coordinates != 3 in FORCE_EVAL%SUBSYS%VELOCITY ("// & 1323 TRIM(ADJUSTL(stmp))//")." 1324 IF (itmp /= 3) & 1325 CPABORT(msg) 1326 DO ib = 1, pint_env%p 1327 DO ic = 1, 3 1328 idim = 3*(ia - 1) + ic 1329 pint_env%v(ib, idim) = r_vals(ic)*unit_conv 1330 END DO 1331 END DO 1332 END DO 1333 1334 vels_present = .TRUE. 1335 END IF 1336 1337 ! set the actual temperature... 1338 IF (vels_present) THEN 1339 ! ...from the initial velocities 1340 ek = 0.0_dp 1341 DO ia = 1, pint_env%ndim/3 1342 rtmp = 0.0_dp 1343 DO ic = 1, 3 1344 idim = 3*(ia - 1) + ic 1345 rtmp = rtmp + pint_env%v(1, idim)*pint_env%v(1, idim) 1346 END DO 1347 ek = ek + 0.5_dp*pint_env%mass(idim)*rtmp 1348 END DO 1349 actual_t = 2.0_dp*ek/pint_env%ndim 1350 ELSE 1351 ! ...using the temperature value from the input 1352 actual_t = pint_env%kT 1353 END IF 1354 1355 ! set the target temperature 1356 target_t = pint_env%kT 1357 CALL section_vals_val_get(pint_env%input, & 1358 "MOTION%PINT%INIT%VELOCITY_SCALE", & 1359 l_val=done_scale) 1360 IF (vels_present) THEN 1361 IF (done_scale) THEN 1362 ! rescale the velocities to match the target temperature 1363 rtmp = SQRT(target_t/actual_t) 1364 DO ia = 1, pint_env%ndim/3 1365 DO ib = 1, pint_env%p 1366 DO ic = 1, 3 1367 idim = 3*(ia - 1) + ic 1368 pint_env%v(ib, idim) = rtmp*pint_env%v(ib, idim) 1369 END DO 1370 END DO 1371 END DO 1372 ELSE 1373 target_t = actual_t 1374 END IF 1375 END IF 1376 1377 ! draw velocities from the M-B distribution... 1378 IF (vels_present) THEN 1379 ! ...for non-centroid modes only 1380 CALL pint_x2u(pint_env, ux=pint_env%uv, x=pint_env%v) 1381 first_mode = 2 1382 ELSE 1383 ! ...for all the modes 1384 first_mode = 1 1385 END IF 1386 DO idim = 1, SIZE(pint_env%uv, 2) 1387 DO ib = first_mode, SIZE(pint_env%uv, 1) 1388 pint_env%uv(ib, idim) = & 1389 next_random_number(rng_stream=pint_env%randomG, & 1390 variance=target_t/pint_env%mass_fict(ib, idim)) 1391 END DO 1392 END DO 1393 1394 ! add random component to the centroid velocity if requsted 1395 done_sped = .FALSE. 1396 CALL section_vals_val_get(pint_env%input, & 1397 "MOTION%PINT%INIT%CENTROID_SPEED", & 1398 l_val=ltmp) 1399 IF (ltmp) THEN 1400 CALL pint_u2x(pint_env, ux=pint_env%uv, x=pint_env%v) 1401 DO idim = 1, pint_env%ndim 1402 rtmp = next_random_number(rng_stream=pint_env%randomG, & 1403 variance=pint_env%mass(idim)*pint_env%kT) & 1404 /pint_env%mass(idim) 1405 DO ib = 1, pint_env%p 1406 pint_env%v(ib, idim) = pint_env%v(ib, idim) + rtmp 1407 END DO 1408 END DO 1409 CALL pint_x2u(pint_env, ux=pint_env%uv, x=pint_env%v) 1410 done_sped = .TRUE. 1411 END IF 1412 1413 ! quench (set to zero) velocities for all the modes if requested 1414 ! (disregard all the initialization done so far) 1415 done_quench = .FALSE. 1416 CALL section_vals_val_get(pint_env%input, & 1417 "MOTION%PINT%INIT%VELOCITY_QUENCH", & 1418 l_val=ltmp) 1419 IF (ltmp) THEN 1420 DO idim = 1, pint_env%ndim 1421 DO ib = 1, pint_env%p 1422 pint_env%v(ib, idim) = 0.0_dp 1423 END DO 1424 END DO 1425 CALL pint_x2u(pint_env, ux=pint_env%uv, x=pint_env%v) 1426 done_quench = .TRUE. 1427 END IF 1428 1429 ! set the velocities to the values from the input if they are explicit 1430 ! (disregard all the initialization done so far) 1431 done_init = .FALSE. 1432 NULLIFY (input_section) 1433 input_section => section_vals_get_subs_vals(pint_env%input, & 1434 "MOTION%PINT%BEADS%VELOCITY") 1435 CALL section_vals_get(input_section, explicit=explicit) 1436 IF (explicit) THEN 1437 CALL section_vals_val_get(input_section, "_DEFAULT_KEYWORD_", & 1438 n_rep_val=n_rep_val) 1439 IF (n_rep_val > 0) THEN 1440 CPASSERT(n_rep_val == 1) 1441 CALL section_vals_val_get(input_section, "_DEFAULT_KEYWORD_", & 1442 r_vals=r_vals) 1443 IF (SIZE(r_vals) /= pint_env%p*pint_env%ndim) & 1444 CPABORT("Invalid size of MOTION%PINT%BEAD%VELOCITY") 1445 itmp = 0 1446 DO idim = 1, pint_env%ndim 1447 DO ib = 1, pint_env%p 1448 itmp = itmp + 1 1449 pint_env%v(ib, idim) = r_vals(itmp) 1450 END DO 1451 END DO 1452 CALL pint_x2u(pint_env, ux=pint_env%uv, x=pint_env%v) 1453 done_init = .TRUE. 1454 END IF 1455 END IF 1456 1457 unit_conv = cp_unit_from_cp2k(1.0_dp, "K") 1458 WRITE (stmp1, '(F10.2)') target_t*pint_env%propagator%temp_sim2phys*unit_conv 1459 msg = "Bead velocities initialization:" 1460 IF (done_init) THEN 1461 msg = TRIM(msg)//" input structure" 1462 ELSE IF (done_quench) THEN 1463 msg = TRIM(msg)//" quenching (set to 0.0)" 1464 ELSE 1465 IF (vels_present) THEN 1466 msg = TRIM(ADJUSTL(msg))//" centroid +" 1467 END IF 1468 msg = TRIM(ADJUSTL(msg))//" Maxwell-Boltzmann at "//TRIM(ADJUSTL(stmp1))//" K." 1469 END IF 1470 CALL pint_write_line(msg) 1471 1472 IF (done_init .AND. done_quench) THEN 1473 msg = "WARNING: exclusive options requested (velocity restart and quenching)" 1474 CPWARN(msg) 1475 msg = "WARNING: velocity restart took precedence" 1476 CPWARN(msg) 1477 END IF 1478 1479 IF ((.NOT. done_init) .AND. (.NOT. done_quench)) THEN 1480 IF (vels_present .AND. done_scale) THEN 1481 WRITE (stmp1, '(F10.2)') actual_t*unit_conv 1482 WRITE (stmp2, '(F10.2)') target_t*unit_conv 1483 msg = "Scaled initial velocities from "//TRIM(ADJUSTL(stmp1))// & 1484 " to "//TRIM(ADJUSTL(stmp2))//" K as requested." 1485 CPWARN(msg) 1486 END IF 1487 IF (done_sped) THEN 1488 msg = "Added random component to the initial centroid velocities." 1489 CPWARN(msg) 1490 END IF 1491 END IF 1492 1493 ! Apply constraints to the initial velocities for centroid constraints 1494 IF (pint_env%simpar%constraint) THEN 1495 IF (pint_env%propagator%prop_kind == propagator_rpmd) THEN 1496 ! Multiply with 1/SQRT(n_beads) due to normal mode transformation in RPMD 1497 factor = SQRT(REAL(pint_env%p, dp)) 1498 ELSE 1499 ! lowest NM is centroid 1500 factor = 1.0_dp 1501 END IF 1502 1503 IF (pint_env%logger%para_env%ionode) THEN 1504 DO i = 1, nparticle 1505 DO j = 1, 3 1506 vel(j, i) = pint_env%uv(1, j + (i - 1)*3)/factor 1507 END DO 1508 END DO 1509 1510 ! Possibly update the target values 1511 CALL shake_update_targets(gci, local_molecules, molecule_set, & 1512 molecule_kind_set, pint_env%dt, & 1513 f_env%force_env%root_section) 1514 CALL rattle_control(gci, local_molecules, molecule_set, & 1515 molecule_kind_set, particle_set, & 1516 vel, pint_env%dt, pint_env%simpar%shake_tol, & 1517 pint_env%simpar%info_constraint, & 1518 pint_env%simpar%lagrange_multipliers, & 1519 .FALSE., & 1520 cell, mp_comm_self, & 1521 local_particles) 1522 END IF 1523 1524 ! Broadcast updated velocities to other nodes 1525 CALL mp_bcast(vel, & 1526 pint_env%logger%para_env%source, & 1527 pint_env%logger%para_env%group) 1528 1529 DO i = 1, nparticle 1530 DO j = 1, 3 1531 pint_env%uv(1, j + (i - 1)*3) = vel(j, i)*factor 1532 END DO 1533 END DO 1534 1535 END IF 1536 1537 RETURN 1538 END SUBROUTINE pint_init_v 1539 1540! *************************************************************************** 1541!> \brief Assign initial postions and velocities to the thermostats. 1542!> \param pint_env ... 1543!> \param kT ... 1544!> \date 2010-12-15 1545!> \author Lukasz Walewski 1546!> \note Extracted from pint_init 1547! ************************************************************************************************** 1548 SUBROUTINE pint_init_t(pint_env, kT) 1549 1550 TYPE(pint_env_type), POINTER :: pint_env 1551 REAL(kind=dp), INTENT(in), OPTIONAL :: kT 1552 1553 CHARACTER(len=*), PARAMETER :: routineN = 'pint_init_t', routineP = moduleN//':'//routineN 1554 1555 INTEGER :: ib, idim, ii, inos, n_rep_val 1556 LOGICAL :: explicit, gle_restart 1557 REAL(kind=dp) :: mykt 1558 REAL(kind=dp), DIMENSION(:), POINTER :: r_vals 1559 TYPE(section_vals_type), POINTER :: input_section 1560 1561 CPASSERT(ASSOCIATED(pint_env)) 1562 CPASSERT(pint_env%ref_count > 0) 1563 1564 IF (pint_env%pimd_thermostat == thermostat_nose) THEN 1565 1566 mykt = pint_env%kT 1567 IF (PRESENT(kT)) mykt = kT 1568 DO idim = 1, SIZE(pint_env%tv, 3) 1569 DO ib = 1, SIZE(pint_env%tv, 2) 1570 DO inos = 1, SIZE(pint_env%tv, 1) 1571 pint_env%tv(inos, ib, idim) = & 1572 next_random_number(rng_stream=pint_env%randomG, & 1573 variance=mykt/pint_env%Q(ib)) 1574 END DO 1575 END DO 1576 END DO 1577 1578 NULLIFY (input_section) 1579 input_section => section_vals_get_subs_vals(pint_env%input, & 1580 "MOTION%PINT%NOSE%COORD") 1581 CALL section_vals_get(input_section, explicit=explicit) 1582 IF (explicit) THEN 1583 CALL section_vals_val_get(input_section, "_DEFAULT_KEYWORD_", & 1584 n_rep_val=n_rep_val) 1585 IF (n_rep_val > 0) THEN 1586 CPASSERT(n_rep_val == 1) 1587 CALL section_vals_val_get(input_section, "_DEFAULT_KEYWORD_", & 1588 r_vals=r_vals) 1589 IF (SIZE(r_vals) /= pint_env%p*pint_env%ndim*pint_env%nnos) & 1590 CPABORT("Invalid size of MOTION%PINT%NOSE%COORD") 1591 ii = 0 1592 DO idim = 1, pint_env%ndim 1593 DO ib = 1, pint_env%p 1594 DO inos = 1, pint_env%nnos 1595 ii = ii + 1 1596 pint_env%tx(inos, ib, idim) = r_vals(ii) 1597 END DO 1598 END DO 1599 END DO 1600 END IF 1601 END IF 1602 1603 NULLIFY (input_section) 1604 input_section => section_vals_get_subs_vals(pint_env%input, & 1605 "MOTION%PINT%NOSE%VELOCITY") 1606 CALL section_vals_get(input_section, explicit=explicit) 1607 IF (explicit) THEN 1608 CALL section_vals_val_get(input_section, "_DEFAULT_KEYWORD_", & 1609 n_rep_val=n_rep_val) 1610 IF (n_rep_val > 0) THEN 1611 CPASSERT(n_rep_val == 1) 1612 CALL section_vals_val_get(input_section, "_DEFAULT_KEYWORD_", & 1613 r_vals=r_vals) 1614 IF (SIZE(r_vals) /= pint_env%p*pint_env%ndim*pint_env%nnos) & 1615 CPABORT("Invalid size of MOTION%PINT%NOSE%VELOCITY") 1616 ii = 0 1617 DO idim = 1, pint_env%ndim 1618 DO ib = 1, pint_env%p 1619 DO inos = 1, pint_env%nnos 1620 ii = ii + 1 1621 pint_env%tv(inos, ib, idim) = r_vals(ii) 1622 END DO 1623 END DO 1624 END DO 1625 END IF 1626 END IF 1627 1628 ELSEIF (pint_env%pimd_thermostat == thermostat_gle) THEN 1629 NULLIFY (input_section) 1630 input_section => section_vals_get_subs_vals(pint_env%input, & 1631 "MOTION%PINT%GLE") 1632 CALL section_vals_get(input_section, explicit=explicit) 1633 IF (explicit) THEN 1634 CALL restart_gle(pint_env%gle, input_section, save_mem=.FALSE., & 1635 restart=gle_restart) 1636 END IF 1637 END IF 1638 1639 RETURN 1640 END SUBROUTINE pint_init_t 1641 1642! *************************************************************************** 1643!> \brief Prepares the forces, etc. to perform an PIMD step 1644!> \param pint_env ... 1645!> \param helium_env ... 1646!> \par History 1647!> Added nh_energy calculation [hforbert] 1648!> Bug fixes for no thermostats [hforbert] 1649!> 2016-07-14 Modified to work with independent helium_env [cschran] 1650!> \author fawzi 1651! ************************************************************************************************** 1652 SUBROUTINE pint_init_f(pint_env, helium_env) 1653 TYPE(pint_env_type), POINTER :: pint_env 1654 TYPE(helium_solvent_p_type), DIMENSION(:), & 1655 OPTIONAL, POINTER :: helium_env 1656 1657 CHARACTER(len=*), PARAMETER :: routineN = 'pint_init_f', routineP = moduleN//':'//routineN 1658 1659 INTEGER :: ib, idim, inos 1660 REAL(kind=dp) :: e_h 1661 TYPE(cp_logger_type), POINTER :: logger 1662 1663 CPASSERT(ASSOCIATED(pint_env)) 1664 CPASSERT(pint_env%ref_count > 0) 1665 1666 NULLIFY (logger) 1667 logger => cp_get_default_logger() 1668 1669 ! initialize iteration info 1670 CALL cp_iterate(logger%iter_info, iter_nr=pint_env%first_step) 1671 CALL cp_iterate(pint_env%logger%iter_info, iter_nr=pint_env%first_step) 1672 1673 CALL pint_x2u(pint_env) 1674 CALL pint_calc_uf_h(pint_env=pint_env, e_h=e_h) 1675 CALL pint_calc_f(pint_env) 1676 1677 ! add helium forces to the solute's internal ones 1678 ! Assume that helium has been already initialized and helium_env(1) 1679 ! contains proper forces in force_avrg array at ionode 1680 IF (PRESENT(helium_env)) THEN 1681 IF (logger%para_env%ionode) THEN 1682 pint_env%f(:, :) = pint_env%f(:, :) + helium_env(1)%helium%force_avrg(:, :) 1683 END IF 1684 CALL mp_bcast(pint_env%f, & 1685 logger%para_env%source, & 1686 logger%para_env%group) 1687 END IF 1688 CALL pint_f2uf(pint_env) 1689 1690 ! set the centroid forces to 0 if FIX_CENTROID_POS 1691 IF (pint_env%first_propagated_mode .EQ. 2) THEN 1692 pint_env%uf(1, :) = 0.0_dp 1693 END IF 1694 1695 CALL pint_calc_e_kin_beads_u(pint_env) 1696 CALL pint_calc_e_vir(pint_env) 1697 DO idim = 1, SIZE(pint_env%uf_h, 2) 1698 DO ib = pint_env%first_propagated_mode, SIZE(pint_env%uf_h, 1) 1699 pint_env%uf(ib, idim) = REAL(pint_env%nrespa, dp)*pint_env%uf(ib, idim) 1700 END DO 1701 END DO 1702 1703 IF (pint_env%nnos > 0) THEN 1704 DO idim = 1, SIZE(pint_env%uf_h, 2) 1705 DO ib = 1, SIZE(pint_env%uf_h, 1) 1706 pint_env%tf(1, ib, idim) = (pint_env%mass_fict(ib, idim)* & 1707 pint_env%uv(ib, idim)**2 - pint_env%kT)/pint_env%Q(ib) 1708 END DO 1709 END DO 1710 1711 DO idim = 1, pint_env%ndim 1712 DO ib = 1, pint_env%p 1713 DO inos = 1, pint_env%nnos - 1 1714 pint_env%tf(inos + 1, ib, idim) = pint_env%tv(inos, ib, idim)**2 - & 1715 pint_env%kT/pint_env%Q(ib) 1716 END DO 1717 DO inos = 1, pint_env%nnos - 1 1718 pint_env%tf(inos, ib, idim) = pint_env%tf(inos, ib, idim) & 1719 - pint_env%tv(inos, ib, idim)*pint_env%tv(inos + 1, ib, idim) 1720 END DO 1721 END DO 1722 END DO 1723 CALL pint_calc_nh_energy(pint_env) 1724 END IF 1725 1726 RETURN 1727 END SUBROUTINE pint_init_f 1728 1729! *************************************************************************** 1730!> \brief Perform the PIMD simulation (main MD loop) 1731!> \param pint_env ... 1732!> \param globenv ... 1733!> \param helium_env ... 1734!> \par History 1735!> 2003-11 created [fawzi] 1736!> renamed from pint_run to pint_do_run because of conflicting name 1737!> of pint_run in input_constants [hforbert] 1738!> 2009-12-14 globenv parameter added to handle soft exit 1739!> requests [lwalewski] 1740!> 2016-07-14 Modified to work with independent helium_env [cschran] 1741!> \author Fawzi Mohamed 1742!> \note Everything should be read for an md step. 1743! ************************************************************************************************** 1744 SUBROUTINE pint_do_run(pint_env, globenv, helium_env) 1745 TYPE(pint_env_type), POINTER :: pint_env 1746 TYPE(global_environment_type), POINTER :: globenv 1747 TYPE(helium_solvent_p_type), DIMENSION(:), & 1748 OPTIONAL, POINTER :: helium_env 1749 1750 CHARACTER(len=*), PARAMETER :: routineN = 'pint_do_run', routineP = moduleN//':'//routineN 1751 1752 INTEGER :: k, step 1753 LOGICAL :: should_stop 1754 REAL(kind=dp) :: scal 1755 TYPE(cp_logger_type), POINTER :: logger 1756 TYPE(f_env_type), POINTER :: f_env 1757 1758 CPASSERT(ASSOCIATED(pint_env)) 1759 1760 ! initialize iteration info 1761 CALL cp_iterate(pint_env%logger%iter_info, iter_nr=pint_env%first_step) 1762 1763 ! iterate replica pint counter by accessing the globally saved 1764 ! force environment error/logger variables and setting them 1765 ! explicitly to the pimd "PINT" step value 1766 CALL f_env_add_defaults(f_env_id=pint_env%replicas%f_env_id, & 1767 f_env=f_env) 1768 NULLIFY (logger) 1769 logger => cp_get_default_logger() 1770 CALL cp_iterate(logger%iter_info, & 1771 iter_nr=pint_env%first_step) 1772 CALL f_env_rm_defaults(f_env) 1773 1774 pint_env%iter = pint_env%first_step 1775 1776 IF (PRESENT(helium_env)) THEN 1777 IF (ASSOCIATED(helium_env)) THEN 1778 ! set the properties accumulated over the whole MC process to 0 1779 DO k = 1, SIZE(helium_env) 1780 helium_env(k)%helium%proarea%accu(:) = 0.0_dp 1781 helium_env(k)%helium%prarea2%accu(:) = 0.0_dp 1782 helium_env(k)%helium%wnmber2%accu(:) = 0.0_dp 1783 helium_env(k)%helium%mominer%accu(:) = 0.0_dp 1784 IF (helium_env(k)%helium%rho_present) THEN 1785 helium_env(k)%helium%rho_accu(:, :, :, :) = 0.0_dp 1786 END IF 1787 IF (helium_env(k)%helium%rdf_present) THEN 1788 helium_env(k)%helium%rdf_accu(:, :) = 0.0_dp 1789 END IF 1790 END DO 1791 END IF 1792 END IF 1793 1794 ! write the properties at 0-th step 1795 CALL pint_calc_energy(pint_env) 1796 CALL pint_calc_total_action(pint_env) 1797 CALL pint_write_ener(pint_env) 1798 CALL pint_write_action(pint_env) 1799 CALL pint_write_centroids(pint_env) 1800 CALL pint_write_trajectory(pint_env) 1801 CALL pint_write_com(pint_env) 1802 CALL pint_write_rgyr(pint_env) 1803 1804 ! main PIMD loop 1805 DO step = 1, pint_env%num_steps 1806 1807 pint_env%iter = pint_env%iter + 1 1808 CALL cp_iterate(pint_env%logger%iter_info, & 1809 last=(step == pint_env%num_steps), & 1810 iter_nr=pint_env%iter) 1811 CALL cp_iterate(logger%iter_info, & 1812 last=(step == pint_env%num_steps), & 1813 iter_nr=pint_env%iter) 1814 pint_env%t = pint_env%t + pint_env%dt 1815 1816 IF (pint_env%t_tol > 0.0_dp) THEN 1817 IF (ABS(2._dp*pint_env%e_kin_beads/(pint_env%p*pint_env%ndim) & 1818 - pint_env%kT) > pint_env%t_tol) THEN 1819 scal = SQRT(pint_env%kT*(pint_env%p*pint_env%ndim)/(2.0_dp*pint_env%e_kin_beads)) 1820 pint_env%uv = scal*pint_env%uv 1821 CALL pint_init_f(pint_env, helium_env=helium_env) 1822 END IF 1823 END IF 1824 CALL pint_step(pint_env, helium_env=helium_env) 1825 1826 CALL pint_write_ener(pint_env) 1827 CALL pint_write_action(pint_env) 1828 CALL pint_write_centroids(pint_env) 1829 CALL pint_write_trajectory(pint_env) 1830 CALL pint_write_com(pint_env) 1831 CALL pint_write_rgyr(pint_env) 1832 1833 CALL write_restart(root_section=pint_env%input, & 1834 pint_env=pint_env, helium_env=helium_env) 1835 1836 ! exit from the main loop if soft exit has been requested 1837 CALL external_control(should_stop, "PINT", globenv=globenv) 1838 IF (should_stop) EXIT 1839 1840 END DO 1841 1842 ! remove iteration level 1843 CALL cp_rm_iter_level(pint_env%logger%iter_info, "PINT") 1844 1845 RETURN 1846 END SUBROUTINE pint_do_run 1847 1848! *************************************************************************** 1849!> \brief Performs a scan of the helium-solute interaction energy 1850!> \param pint_env ... 1851!> \param helium_env ... 1852!> \date 2013-11-26 1853!> \parm History 1854!> 2016-07-14 Modified to work with independent helium_env [cschran] 1855!> \author Lukasz Walewski 1856! ************************************************************************************************** 1857 SUBROUTINE pint_run_scan(pint_env, helium_env) 1858 TYPE(pint_env_type), POINTER :: pint_env 1859 TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env 1860 1861 CHARACTER(len=*), PARAMETER :: routineN = 'pint_run_scan', routineP = moduleN//':'//routineN 1862 1863 CHARACTER(len=default_string_length) :: comment 1864 INTEGER :: unit_nr 1865 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: DATA 1866 TYPE(section_vals_type), POINTER :: print_key 1867 1868 NULLIFY (pint_env%logger, print_key) 1869 pint_env%logger => cp_get_default_logger() 1870 1871 ! assume that ionode always has at least one helium_env 1872 IF (pint_env%logger%para_env%ionode) THEN 1873 print_key => section_vals_get_subs_vals(helium_env(1)%helium%input, & 1874 "MOTION%PINT%HELIUM%PRINT%RHO") 1875 END IF 1876 1877 ! perform the actual scan wrt the COM of the solute 1878 CALL helium_intpot_scan(pint_env, helium_env) 1879 1880 ! output the interaction potential into a cubefile 1881 ! assume that ionode always has at least one helium_env 1882 IF (pint_env%logger%para_env%ionode) THEN 1883 1884 unit_nr = cp_print_key_unit_nr( & 1885 pint_env%logger, & 1886 print_key, & 1887 middle_name="helium-pot", & 1888 extension=".cube", & 1889 file_position="REWIND", & 1890 do_backup=.FALSE.) 1891 1892 comment = "Solute - helium interaction potential" 1893 NULLIFY (DATA) 1894 DATA => helium_env(1)%helium%rho_inst(1, :, :, :) 1895 CALL helium_write_cubefile( & 1896 unit_nr, & 1897 comment, & 1898 helium_env(1)%helium%center - 0.5_dp* & 1899 (helium_env(1)%helium%rho_maxr - helium_env(1)%helium%rho_delr), & 1900 helium_env(1)%helium%rho_delr, & 1901 helium_env(1)%helium%rho_nbin, & 1902 DATA) 1903 1904 CALL m_flush(unit_nr) 1905 CALL cp_print_key_finished_output(unit_nr, pint_env%logger, print_key) 1906 1907 END IF 1908 1909 ! output solute positions 1910 CALL pint_write_centroids(pint_env) 1911 CALL pint_write_trajectory(pint_env) 1912 1913 RETURN 1914 END SUBROUTINE pint_run_scan 1915 1916! *************************************************************************** 1917!> \brief Does an PINT step (and nrespa harmonic evaluations) 1918!> \param pint_env ... 1919!> \param helium_env ... 1920!> \par History 1921!> various bug fixes [hforbert] 1922!> 10.2015 Added RPMD propagator and harmonic integrator [Felix Uhl] 1923!> 04.2016 Changed to work with helium_env [cschran] 1924!> 10.2018 Added centroid constraints [cschran+rperez] 1925!> \author fawzi 1926! ************************************************************************************************** 1927 SUBROUTINE pint_step(pint_env, helium_env) 1928 TYPE(pint_env_type), POINTER :: pint_env 1929 TYPE(helium_solvent_p_type), DIMENSION(:), & 1930 OPTIONAL, POINTER :: helium_env 1931 1932 CHARACTER(len=*), PARAMETER :: routineN = 'pint_step', routineP = moduleN//':'//routineN 1933 1934 INTEGER :: handle, i, ia, ib, idim, ierr, inos, & 1935 iresp, j, k, nparticle, nparticle_kind 1936 REAL(kind=dp) :: dt_temp, dti, dti2, dti22, e_h, factor, & 1937 rn, tdti, time_start, time_stop, tol 1938 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: pos, vel 1939 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: tmp 1940 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds 1941 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 1942 TYPE(cell_type), POINTER :: cell 1943 TYPE(cp_subsys_type), POINTER :: subsys 1944 TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles 1945 TYPE(f_env_type), POINTER :: f_env 1946 TYPE(global_constraint_type), POINTER :: gci 1947 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds 1948 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set 1949 TYPE(molecule_list_type), POINTER :: molecules 1950 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 1951 TYPE(particle_list_type), POINTER :: particles 1952 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 1953 1954 CALL timeset(routineN, handle) 1955 time_start = m_walltime() 1956 1957 CPASSERT(ASSOCIATED(pint_env)) 1958 1959 rn = REAL(pint_env%nrespa, dp) 1960 dti = pint_env%dt/rn 1961 dti2 = dti/2._dp 1962 tdti = 2.*dti 1963 dti22 = dti**2/2._dp 1964 1965 ! Get centroid constraint info, if needed 1966 ! Create a force environment which will be identical to 1967 ! the bead that is being processed by the processor. 1968 IF (pint_env%simpar%constraint) THEN 1969 NULLIFY (subsys, cell) 1970 NULLIFY (atomic_kinds, local_particles, particles) 1971 NULLIFY (local_molecules, molecules, molecule_kinds, gci) 1972 NULLIFY (atomic_kind_set, molecule_kind_set, particle_set, molecule_set) 1973 1974 CALL f_env_add_defaults(f_env_id=pint_env%replicas%f_env_id, f_env=f_env) 1975 CALL force_env_get(force_env=f_env%force_env, subsys=subsys) 1976 CALL f_env_rm_defaults(f_env, ierr) 1977 CPASSERT(ierr == 0) 1978 1979 ! Get gci and more from subsys 1980 CALL cp_subsys_get(subsys=subsys, & 1981 cell=cell, & 1982 atomic_kinds=atomic_kinds, & 1983 local_particles=local_particles, & 1984 particles=particles, & 1985 local_molecules=local_molecules, & 1986 molecules=molecules, & 1987 molecule_kinds=molecule_kinds, & 1988 gci=gci) 1989 1990 nparticle_kind = atomic_kinds%n_els 1991 atomic_kind_set => atomic_kinds%els 1992 molecule_kind_set => molecule_kinds%els 1993 nparticle = particles%n_els 1994 particle_set => particles%els 1995 molecule_set => molecules%els 1996 1997 ! Allocate work storage 1998 ALLOCATE (pos(3, nparticle)) 1999 ALLOCATE (vel(3, nparticle)) 2000 pos(:, :) = 0.0_dp 2001 vel(:, :) = 0.0_dp 2002 2003 IF (pint_env%propagator%prop_kind == propagator_rpmd) THEN 2004 ! Multiply with 1/SQRT(n_beads) due to normal mode transformation in RPMD 2005 factor = SQRT(REAL(pint_env%p, dp)) 2006 ELSE 2007 factor = 1.0_dp 2008 END IF 2009 2010 CALL getold(gci, local_molecules, molecule_set, & 2011 molecule_kind_set, particle_set, cell) 2012 END IF 2013 2014 SELECT CASE (pint_env%harm_integrator) 2015 CASE (integrate_numeric) 2016 2017 DO iresp = 1, pint_env%nrespa 2018 2019 ! integrate bead positions, first_propagated_mode = { 1, 2 } 2020 ! Nose needs an extra step 2021 IF (pint_env%pimd_thermostat == thermostat_nose) THEN 2022 2023 !Set thermostat action of constrained DoF to zero: 2024 IF (pint_env%simpar%constraint) THEN 2025 DO k = 1, pint_env%n_atoms_constraints 2026 ia = pint_env%atoms_constraints(k) 2027 DO j = 3*(ia - 1) + 1, 3*ia 2028 pint_env%tv(:, 1, j) = 0.0_dp 2029 END DO 2030 END DO 2031 END IF 2032 2033 DO i = pint_env%first_propagated_mode, pint_env%p 2034 pint_env%ux(i, :) = pint_env%ux(i, :) - & 2035 dti22*pint_env%uv(i, :)*pint_env%tv(1, i, :) 2036 END DO 2037 pint_env%tx = pint_env%tx + dti*pint_env%tv + dti22*pint_env%tf 2038 2039 END IF 2040 !Integrate position in harmonic springs (uf_h) and physical potential 2041 !(uf) 2042 DO i = pint_env%first_propagated_mode, pint_env%p 2043 pint_env%ux_t(i, :) = pint_env%ux(i, :) + & 2044 dti*pint_env%uv(i, :) + & 2045 dti22*(pint_env%uf_h(i, :) + & 2046 pint_env%uf(i, :)) 2047 END DO 2048 2049 ! apply thermostats to velocities 2050 SELECT CASE (pint_env%pimd_thermostat) 2051 CASE (thermostat_nose) 2052 pint_env%uv_t = pint_env%uv - dti2* & 2053 pint_env%uv*pint_env%tv(1, :, :) 2054 tmp => pint_env%tv_t 2055 pint_env%tv_t => pint_env%tv 2056 pint_env%tv => tmp 2057 pint_env%tv = pint_env%tv_old + tdti*pint_env%tf 2058 pint_env%tv_old = pint_env%tv_t 2059 pint_env%tv_t = pint_env%tv_t + dti2*pint_env%tf 2060 CASE DEFAULT 2061 pint_env%uv_t = pint_env%uv 2062 END SELECT 2063 2064 !Set thermostat action of constrained DoF to zero: 2065 IF (pint_env%simpar%constraint) THEN 2066 DO k = 1, pint_env%n_atoms_constraints 2067 ia = pint_env%atoms_constraints(k) 2068 DO j = 3*(ia - 1) + 1, 3*ia 2069 pint_env%tv(:, 1, j) = 0.0_dp 2070 pint_env%tv_t(:, 1, j) = 0.0_dp 2071 END DO 2072 END DO 2073 END IF 2074 2075 !Integrate harmonic velocities and physical velocities 2076 pint_env%uv_t = pint_env%uv_t + dti2*(pint_env%uf_h + pint_env%uf) 2077 2078 ! physical forces are only applied in first respa step. 2079 pint_env%uf = 0.0_dp 2080 ! calc harmonic forces at new pos 2081 pint_env%ux = pint_env%ux_t 2082 2083 ! Apply centroid constraints (SHAKE) 2084 IF (pint_env%simpar%constraint) THEN 2085 IF (pint_env%logger%para_env%ionode) THEN 2086 DO i = 1, nparticle 2087 DO j = 1, 3 2088 pos(j, i) = pint_env%ux(1, j + (i - 1)*3) 2089 vel(j, i) = pint_env%uv_t(1, j + (i - 1)*3) 2090 END DO 2091 END DO 2092 2093 ! Possibly update the target values 2094 CALL shake_update_targets(gci, local_molecules, molecule_set, & 2095 molecule_kind_set, dti, & 2096 f_env%force_env%root_section) 2097 CALL shake_control(gci, local_molecules, molecule_set, & 2098 molecule_kind_set, particle_set, & 2099 pos, vel, dti, pint_env%simpar%shake_tol, & 2100 pint_env%simpar%info_constraint, & 2101 pint_env%simpar%lagrange_multipliers, & 2102 pint_env%simpar%dump_lm, cell, & 2103 mp_comm_self, local_particles) 2104 END IF 2105 ! Positions and velocities of centroid were constrained by SHAKE 2106 CALL mp_bcast(pos, & 2107 pint_env%logger%para_env%source, & 2108 pint_env%logger%para_env%group) 2109 CALL mp_bcast(vel, & 2110 pint_env%logger%para_env%source, & 2111 pint_env%logger%para_env%group) 2112 ! Transform back to normal modes: 2113 DO i = 1, nparticle 2114 DO j = 1, 3 2115 pint_env%ux(1, j + (i - 1)*3) = pos(j, i) 2116 pint_env%uv_t(1, j + (i - 1)*3) = vel(j, i) 2117 END DO 2118 END DO 2119 2120 END IF 2121 2122 CALL pint_calc_uf_h(pint_env=pint_env, e_h=e_h) 2123 pint_env%uv_t = pint_env%uv_t + dti2*(pint_env%uf_h + pint_env%uf) 2124 2125 ! For last respa step include integration of physical and helium 2126 ! forces 2127 IF (iresp == pint_env%nrespa) THEN 2128 CALL pint_u2x(pint_env) 2129 CALL pint_calc_f(pint_env) 2130 ! perform helium step and add helium forces 2131 IF (PRESENT(helium_env)) THEN 2132 CALL helium_step(helium_env, pint_env) 2133 !Update force of solute in pint_env 2134 IF (pint_env%logger%para_env%ionode) THEN 2135 pint_env%f(:, :) = pint_env%f(:, :) + helium_env(1)%helium%force_avrg(:, :) 2136 END IF 2137 CALL mp_bcast(pint_env%f, & 2138 pint_env%logger%para_env%source, & 2139 pint_env%logger%para_env%group) 2140 END IF 2141 2142 CALL pint_f2uf(pint_env) 2143 ! set the centroid forces to 0 if FIX_CENTROID_POS 2144 IF (pint_env%first_propagated_mode .EQ. 2) THEN 2145 pint_env%uf(1, :) = 0.0_dp 2146 END IF 2147 !Scale physical forces and integrate velocities with physical 2148 !forces 2149 pint_env%uf = pint_env%uf*rn 2150 pint_env%uv_t = pint_env%uv_t + dti2*pint_env%uf 2151 2152 END IF 2153 2154 ! Apply second half of thermostats 2155 SELECT CASE (pint_env%pimd_thermostat) 2156 CASE (thermostat_nose) 2157 DO i = 1, 6 2158 tol = 0._dp 2159 pint_env%uv_new = pint_env%uv_t/(1.+dti2*pint_env%tv(1, :, :)) 2160 DO idim = 1, pint_env%ndim 2161 DO ib = 1, pint_env%p 2162 pint_env%tf(1, ib, idim) = (pint_env%mass_fict(ib, idim)* & 2163 pint_env%uv_new(ib, idim)**2 - pint_env%kT*pint_env%kTcorr)/ & 2164 pint_env%Q(ib) 2165 END DO 2166 END DO 2167 2168 !Set thermostat action of constrained DoF to zero: 2169 IF (pint_env%simpar%constraint) THEN 2170 DO k = 1, pint_env%n_atoms_constraints 2171 ia = pint_env%atoms_constraints(k) 2172 DO j = 3*(ia - 1) + 1, 3*ia 2173 pint_env%tf(:, 1, j) = 0.0_dp 2174 END DO 2175 END DO 2176 END IF 2177 2178 DO idim = 1, pint_env%ndim 2179 DO ib = 1, pint_env%p 2180 DO inos = 1, pint_env%nnos - 1 2181 pint_env%tv_new(inos, ib, idim) = & 2182 (pint_env%tv_t(inos, ib, idim) + dti2*pint_env%tf(inos, ib, idim))/ & 2183 (1._dp + dti2*pint_env%tv(inos + 1, ib, idim)) 2184 pint_env%tf(inos + 1, ib, idim) = & 2185 (pint_env%tv_new(inos, ib, idim)**2 - & 2186 pint_env%kT*pint_env%kTcorr/pint_env%Q(ib)) 2187 tol = MAX(tol, ABS(pint_env%tv(inos, ib, idim) & 2188 - pint_env%tv_new(inos, ib, idim))) 2189 END DO 2190 !Set thermostat action of constrained DoF to zero: 2191 IF (pint_env%simpar%constraint) THEN 2192 DO k = 1, pint_env%n_atoms_constraints 2193 ia = pint_env%atoms_constraints(k) 2194 DO j = 3*(ia - 1) + 1, 3*ia 2195 pint_env%tv_new(:, 1, j) = 0.0_dp 2196 pint_env%tf(:, 1, j) = 0.0_dp 2197 END DO 2198 END DO 2199 END IF 2200 2201 pint_env%tv_new(pint_env%nnos, ib, idim) = & 2202 pint_env%tv_t(pint_env%nnos, ib, idim) + & 2203 dti2*pint_env%tf(pint_env%nnos, ib, idim) 2204 tol = MAX(tol, ABS(pint_env%tv(pint_env%nnos, ib, idim) & 2205 - pint_env%tv_new(pint_env%nnos, ib, idim))) 2206 tol = MAX(tol, ABS(pint_env%uv(ib, idim) & 2207 - pint_env%uv_new(ib, idim))) 2208 !Set thermostat action of constrained DoF to zero: 2209 IF (pint_env%simpar%constraint) THEN 2210 DO k = 1, pint_env%n_atoms_constraints 2211 ia = pint_env%atoms_constraints(k) 2212 DO j = 3*(ia - 1) + 1, 3*ia 2213 pint_env%tv_new(:, 1, j) = 0.0_dp 2214 END DO 2215 END DO 2216 END IF 2217 2218 END DO 2219 END DO 2220 2221 pint_env%uv = pint_env%uv_new 2222 pint_env%tv = pint_env%tv_new 2223 IF (tol <= pint_env%v_tol) EXIT 2224 END DO 2225 2226 ! Apply centroid constraints (RATTLE) 2227 IF (pint_env%simpar%constraint) THEN 2228 IF (pint_env%logger%para_env%ionode) THEN 2229 ! Reset particle r, due to force calc: 2230 DO i = 1, nparticle 2231 DO j = 1, 3 2232 vel(j, i) = pint_env%uv(1, j + (i - 1)*3) 2233 particle_set(i)%r(j) = pint_env%ux(1, j + (i - 1)*3) 2234 END DO 2235 END DO 2236 2237 ! Small time step for all small integrations steps 2238 ! Big step for last RESPA 2239 IF (iresp == pint_env%nrespa) THEN 2240 dt_temp = dti 2241 ELSE 2242 dt_temp = dti*rn 2243 END IF 2244 CALL rattle_control(gci, local_molecules, molecule_set, & 2245 molecule_kind_set, particle_set, & 2246 vel, dt_temp, pint_env%simpar%shake_tol, & 2247 pint_env%simpar%info_constraint, & 2248 pint_env%simpar%lagrange_multipliers, & 2249 pint_env%simpar%dump_lm, cell, & 2250 mp_comm_self, local_particles) 2251 END IF 2252 ! Velocities of centroid were constrained by RATTLE 2253 ! Broadcast updated velocities to other nodes 2254 CALL mp_bcast(vel, & 2255 pint_env%logger%para_env%source, & 2256 pint_env%logger%para_env%group) 2257 2258 DO i = 1, nparticle 2259 DO j = 1, 3 2260 pint_env%uv(1, j + (i - 1)*3) = vel(j, i) 2261 END DO 2262 END DO 2263 END IF 2264 2265 DO inos = 1, pint_env%nnos - 1 2266 pint_env%tf(inos, :, :) = pint_env%tf(inos, :, :) & 2267 - pint_env%tv(inos, :, :)*pint_env%tv(inos + 1, :, :) 2268 END DO 2269 2270 CASE (thermostat_gle) 2271 CALL pint_gle_step(pint_env) 2272 pint_env%uv = pint_env%uv_t 2273 CASE DEFAULT 2274 pint_env%uv = pint_env%uv_t 2275 END SELECT 2276 END DO 2277 2278 CASE (integrate_exact) 2279 ! The Liouvillian splitting is as follows: 2280 ! 1. Thermostat 2281 ! 2. 0.5*physical integration 2282 ! 3. Exact harmonic integration ( + apply centroid constraints (SHAKE)) 2283 ! 4. 0.5*physical integration 2284 ! 5. Thermostat 2285 ! 6. Apply centroid constraints (RATTLE) 2286 2287 ! 1. Apply thermostats 2288 SELECT CASE (pint_env%pimd_thermostat) 2289 CASE (thermostat_pile) 2290 CALL pint_pile_step(vold=pint_env%uv, & 2291 vnew=pint_env%uv_t, & 2292 p=pint_env%p, & 2293 ndim=pint_env%ndim, & 2294 first_mode=pint_env%first_propagated_mode, & 2295 masses=pint_env%mass_fict, & 2296 pile_therm=pint_env%pile_therm) 2297 CASE (thermostat_piglet) 2298 CALL pint_piglet_step(vold=pint_env%uv, & 2299 vnew=pint_env%uv_t, & 2300 first_mode=pint_env%first_propagated_mode, & 2301 masses=pint_env%mass_fict, & 2302 piglet_therm=pint_env%piglet_therm) 2303 CASE (thermostat_qtb) 2304 CALL pint_qtb_step(vold=pint_env%uv, & 2305 vnew=pint_env%uv_t, & 2306 p=pint_env%p, & 2307 ndim=pint_env%ndim, & 2308 masses=pint_env%mass_fict, & 2309 qtb_therm=pint_env%qtb_therm) 2310 CASE DEFAULT 2311 pint_env%uv_t = pint_env%uv 2312 END SELECT 2313 2314 ! 2. 1/2*Physical integration 2315 pint_env%uv_t = pint_env%uv_t + dti2*pint_env%uf 2316 2317 ! 3. Exact harmonic integration 2318 IF (pint_env%first_propagated_mode == 1) THEN 2319 ! The centroid is integrated via standard velocity-verlet 2320 ! Commented out code is only there to show similarities to 2321 ! Numeric integrator 2322 pint_env%ux_t(1, :) = pint_env%ux(1, :) + & 2323 dti*pint_env%uv_t(1, :) !+ & 2324 ! dti22*pint_env%uf_h(1, :) 2325 !pint_env%uv_t(1, :) = pint_env%uv_t(1, :)+ & 2326 ! dti2*pint_env%uf_h(1, :) 2327 2328 ! Apply centroid constraints (SHAKE) 2329 IF (pint_env%simpar%constraint) THEN 2330 IF (pint_env%logger%para_env%ionode) THEN 2331 ! Transform positions and velocities to Cartesian coordinates: 2332 DO i = 1, nparticle 2333 DO j = 1, 3 2334 pos(j, i) = pint_env%ux_t(1, j + (i - 1)*3)/factor 2335 vel(j, i) = pint_env%uv_t(1, j + (i - 1)*3)/factor 2336 END DO 2337 END DO 2338 2339 ! Possibly update the target values 2340 CALL shake_update_targets(gci, local_molecules, molecule_set, & 2341 molecule_kind_set, dti, & 2342 f_env%force_env%root_section) 2343 CALL shake_control(gci, local_molecules, molecule_set, & 2344 molecule_kind_set, particle_set, & 2345 pos, vel, dti, pint_env%simpar%shake_tol, & 2346 pint_env%simpar%info_constraint, & 2347 pint_env%simpar%lagrange_multipliers, & 2348 pint_env%simpar%dump_lm, cell, & 2349 mp_comm_self, local_particles) 2350 END IF 2351 ! Positions and velocities of centroid were constrained by SHAKE 2352 CALL mp_bcast(pos, & 2353 pint_env%logger%para_env%source, & 2354 pint_env%logger%para_env%group) 2355 CALL mp_bcast(vel, & 2356 pint_env%logger%para_env%source, & 2357 pint_env%logger%para_env%group) 2358 ! Transform back to normal modes: 2359 DO i = 1, nparticle 2360 DO j = 1, 3 2361 pint_env%ux_t(1, j + (i - 1)*3) = pos(j, i)*factor 2362 pint_env%uv_t(1, j + (i - 1)*3) = vel(j, i)*factor 2363 END DO 2364 END DO 2365 2366 END IF 2367 2368 ELSE 2369 ! set velocities to zero for fixed centroids 2370 pint_env%ux_t(1, :) = pint_env%ux(1, :) 2371 pint_env%uv_t(1, :) = 0.0_dp 2372 END IF 2373 2374 DO i = 2, pint_env%p 2375 pint_env%ux_t(i, :) = pint_env%cosex(i)*pint_env%ux(i, :) & 2376 + pint_env%iwsinex(i)*pint_env%uv_t(i, :) 2377 pint_env%uv_t(i, :) = pint_env%cosex(i)*pint_env%uv_t(i, :) & 2378 - pint_env%wsinex(i)*pint_env%ux(i, :) 2379 END DO 2380 pint_env%ux = pint_env%ux_t 2381 2382 ! 4. 1/2*Physical integration 2383 pint_env%uf = 0.0_dp 2384 CALL pint_u2x(pint_env) 2385 !TODO apply bead-wise constraints 2386 CALL pint_calc_f(pint_env) 2387 ! perform helium step and add helium forces 2388 IF (PRESENT(helium_env)) THEN 2389 CALL helium_step(helium_env, pint_env) 2390 !Update force of solute in pint_env 2391 IF (pint_env%logger%para_env%ionode) THEN 2392 pint_env%f(:, :) = pint_env%f(:, :) + helium_env(1)%helium%force_avrg(:, :) 2393 END IF 2394 CALL mp_bcast(pint_env%f, & 2395 pint_env%logger%para_env%source, & 2396 pint_env%logger%para_env%group) 2397 END IF 2398 CALL pint_f2uf(pint_env) 2399 ! set the centroid forces to 0 if FIX_CENTROID_POS 2400 IF (pint_env%first_propagated_mode .EQ. 2) THEN 2401 pint_env%uf(1, :) = 0.0_dp 2402 END IF 2403 pint_env%uv_t = pint_env%uv_t + dti2*pint_env%uf 2404 2405 ! 5. Apply thermostats 2406 SELECT CASE (pint_env%pimd_thermostat) 2407 CASE (thermostat_pile) 2408 CALL pint_pile_step(vold=pint_env%uv_t, & 2409 vnew=pint_env%uv, & 2410 p=pint_env%p, & 2411 ndim=pint_env%ndim, & 2412 first_mode=pint_env%first_propagated_mode, & 2413 masses=pint_env%mass_fict, & 2414 pile_therm=pint_env%pile_therm) 2415 CASE (thermostat_piglet) 2416 CALL pint_piglet_step(vold=pint_env%uv_t, & 2417 vnew=pint_env%uv, & 2418 first_mode=pint_env%first_propagated_mode, & 2419 masses=pint_env%mass_fict, & 2420 piglet_therm=pint_env%piglet_therm) 2421 CASE (thermostat_qtb) 2422 CALL pint_qtb_step(vold=pint_env%uv_t, & 2423 vnew=pint_env%uv, & 2424 p=pint_env%p, & 2425 ndim=pint_env%ndim, & 2426 masses=pint_env%mass_fict, & 2427 qtb_therm=pint_env%qtb_therm) 2428 CASE DEFAULT 2429 pint_env%uv = pint_env%uv_t 2430 END SELECT 2431 2432 ! 6. Apply centroid constraints (RATTLE) 2433 IF (pint_env%simpar%constraint) THEN 2434 IF (pint_env%logger%para_env%ionode) THEN 2435 ! Transform positions and velocities to Cartesian coordinates: 2436 ! Reset particle r, due to force calc: 2437 DO i = 1, nparticle 2438 DO j = 1, 3 2439 vel(j, i) = pint_env%uv(1, j + (i - 1)*3)/factor 2440 particle_set(i)%r(j) = pint_env%ux(1, j + (i - 1)*3)/factor 2441 END DO 2442 END DO 2443 CALL rattle_control(gci, local_molecules, & 2444 molecule_set, molecule_kind_set, & 2445 particle_set, vel, dti, & 2446 pint_env%simpar%shake_tol, & 2447 pint_env%simpar%info_constraint, & 2448 pint_env%simpar%lagrange_multipliers, & 2449 pint_env%simpar%dump_lm, cell, & 2450 mp_comm_self, local_particles) 2451 END IF 2452 ! Velocities of centroid were constrained by RATTLE 2453 ! Broadcast updated velocities to other nodes 2454 CALL mp_bcast(vel, & 2455 pint_env%logger%para_env%source, & 2456 pint_env%logger%para_env%group) 2457 2458 ! Transform back to normal modes: 2459 ! Multiply with SQRT(n_beads) due to normal mode transformation 2460 DO i = 1, nparticle 2461 DO j = 1, 3 2462 pint_env%uv(1, j + (i - 1)*3) = vel(j, i)*factor 2463 END DO 2464 END DO 2465 2466 !TODO apply bead-wise constraints 2467 END IF 2468 2469 END SELECT 2470 2471 IF (pint_env%simpar%constraint) DEALLOCATE (pos, vel) 2472 2473 ! calculate the energy components 2474 CALL pint_calc_energy(pint_env) 2475 CALL pint_calc_total_action(pint_env) 2476 2477 ! check that the number of PINT steps matches 2478 ! the number of force evaluations done so far 2479!TODO make this check valid if we start from ITERATION != 0 2480! CALL f_env_add_defaults(f_env_id=pint_env%replicas%f_env_id,& 2481! f_env=f_env,new_error=new_error) 2482! NULLIFY(logger) 2483! logger => cp_error_get_logger(new_error) 2484! IF(logger%iter_info%iteration(2)/=pint_env%iter+1)& 2485! CPABORT("md & force_eval lost sychro") 2486! CALL f_env_rm_defaults(f_env,new_error,ierr) 2487 2488 time_stop = m_walltime() 2489 pint_env%time_per_step = time_stop - time_start 2490 CALL pint_write_step_info(pint_env) 2491 CALL timestop(handle) 2492 2493 RETURN 2494 END SUBROUTINE pint_step 2495 2496! *************************************************************************** 2497!> \brief Calculate the energy components (private wrapper function) 2498!> \param pint_env ... 2499!> \date 2011-01-07 2500!> \author Lukasz Walewski 2501! ************************************************************************************************** 2502 SUBROUTINE pint_calc_energy(pint_env) 2503 2504 TYPE(pint_env_type), POINTER :: pint_env 2505 2506 CHARACTER(len=*), PARAMETER :: routineN = 'pint_calc_energy', & 2507 routineP = moduleN//':'//routineN 2508 2509 REAL(KIND=dp) :: e_h 2510 2511 CALL pint_calc_e_kin_beads_u(pint_env) 2512 CALL pint_calc_e_vir(pint_env) 2513 2514 CALL pint_calc_uf_h(pint_env, e_h) 2515 pint_env%e_pot_h = e_h 2516 2517 SELECT CASE (pint_env%pimd_thermostat) 2518 CASE (thermostat_nose) 2519 CALL pint_calc_nh_energy(pint_env) 2520 CASE (thermostat_gle) 2521 CALL pint_calc_gle_energy(pint_env) 2522 CASE (thermostat_pile) 2523 CALL pint_calc_pile_energy(pint_env) 2524 CASE (thermostat_qtb) 2525 CALL pint_calc_qtb_energy(pint_env) 2526 CASE (thermostat_piglet) 2527 CALL pint_calc_piglet_energy(pint_env) 2528 END SELECT 2529 2530 pint_env%energy(e_kin_thermo_id) = & 2531 (0.5_dp*REAL(pint_env%p, dp)*REAL(pint_env%ndim, dp)*pint_env%kT - & 2532 pint_env%e_pot_h)*pint_env%propagator%temp_sim2phys 2533 2534 pint_env%energy(e_potential_id) = SUM(pint_env%e_pot_bead) 2535 2536 pint_env%energy(e_conserved_id) = & 2537 pint_env%energy(e_potential_id)*pint_env%propagator%physpotscale + & 2538 pint_env%e_pot_h + & 2539 pint_env%e_kin_beads + & 2540 pint_env%e_pot_t + & 2541 pint_env%e_kin_t + & 2542 pint_env%e_gle + pint_env%e_pile + pint_env%e_piglet + pint_env%e_qtb 2543 2544 pint_env%energy(e_potential_id) = & 2545 pint_env%energy(e_potential_id)/REAL(pint_env%p, dp) 2546 2547 RETURN 2548 END SUBROUTINE pint_calc_energy 2549 2550! *************************************************************************** 2551!> \brief Calculate the harmonic force in the u basis 2552!> \param pint_env the path integral environment in which the harmonic 2553!> forces should be calculated 2554!> \param e_h ... 2555!> \par History 2556!> Added normal mode transformation [hforbert] 2557!> \author fawzi 2558! ************************************************************************************************** 2559 SUBROUTINE pint_calc_uf_h(pint_env, e_h) 2560 TYPE(pint_env_type), POINTER :: pint_env 2561 REAL(KIND=dp), INTENT(OUT) :: e_h 2562 2563 CHARACTER(len=*), PARAMETER :: routineN = 'pint_calc_uf_h', routineP = moduleN//':'//routineN 2564 2565 IF (pint_env%transform == transformation_stage) THEN 2566 CALL staging_calc_uf_h(pint_env%staging_env, & 2567 pint_env%mass_beads, & 2568 pint_env%ux, & 2569 pint_env%uf_h, & 2570 pint_env%e_pot_h) 2571 ELSE 2572 CALL normalmode_calc_uf_h(pint_env%normalmode_env, & 2573 pint_env%mass_beads, & 2574 pint_env%ux, & 2575 pint_env%uf_h, & 2576 pint_env%e_pot_h) 2577 END IF 2578 e_h = pint_env%e_pot_h 2579 pint_env%uf_h = pint_env%uf_h/pint_env%mass_fict 2580 RETURN 2581 END SUBROUTINE pint_calc_uf_h 2582 2583! *************************************************************************** 2584!> \brief calculates the force (and energy) in each bead, returns the sum 2585!> of the potential energy 2586!> \param pint_env path integral environment on which you want to calculate 2587!> the forces 2588!> \param x positions at which you want to evaluate the forces 2589!> \param f the forces 2590!> \param e potential energy on each bead 2591!> \par History 2592!> 2009-06-15 moved helium calls out from here [lwalewski] 2593!> \author fawzi 2594! ************************************************************************************************** 2595 SUBROUTINE pint_calc_f(pint_env, x, f, e) 2596 TYPE(pint_env_type), POINTER :: pint_env 2597 REAL(kind=dp), DIMENSION(:, :), INTENT(in), & 2598 OPTIONAL, TARGET :: x 2599 REAL(kind=dp), DIMENSION(:, :), INTENT(out), & 2600 OPTIONAL, TARGET :: f 2601 REAL(kind=dp), DIMENSION(:), INTENT(out), & 2602 OPTIONAL, TARGET :: e 2603 2604 CHARACTER(len=*), PARAMETER :: routineN = 'pint_calc_f', routineP = moduleN//':'//routineN 2605 2606 INTEGER :: ib, idim 2607 REAL(kind=dp), DIMENSION(:), POINTER :: my_e 2608 REAL(kind=dp), DIMENSION(:, :), POINTER :: my_f, my_x 2609 2610 CPASSERT(ASSOCIATED(pint_env)) 2611 CPASSERT(pint_env%ref_count > 0) 2612 my_x => pint_env%x 2613 IF (PRESENT(x)) my_x => x 2614 my_f => pint_env%f 2615 IF (PRESENT(f)) my_f => f 2616 my_e => pint_env%e_pot_bead 2617 IF (PRESENT(e)) my_e => e 2618 DO idim = 1, pint_env%ndim 2619 DO ib = 1, pint_env%p 2620 pint_env%replicas%r(idim, ib) = my_x(ib, idim) 2621 END DO 2622 END DO 2623 CALL rep_env_calc_e_f(pint_env%replicas, calc_f=.TRUE.) 2624 DO idim = 1, pint_env%ndim 2625 DO ib = 1, pint_env%p 2626 !ljw: is that fine ? - idim <-> ib 2627 my_f(ib, idim) = pint_env%replicas%f(idim, ib) 2628 END DO 2629 END DO 2630 my_e = pint_env%replicas%f(SIZE(pint_env%replicas%f, 1), :) 2631 2632 RETURN 2633 END SUBROUTINE pint_calc_f 2634 2635! *************************************************************************** 2636!> \brief Calculate the kinetic energy of the beads (in the u variables) 2637!> \param pint_env ... 2638!> \param uv ... 2639!> \param e_k ... 2640!> \par History 2641!> Bug fix to give my_uv a default location if not given in call [hforbert] 2642!> \author fawzi 2643! ************************************************************************************************** 2644 SUBROUTINE pint_calc_e_kin_beads_u(pint_env, uv, e_k) 2645 TYPE(pint_env_type), POINTER :: pint_env 2646 REAL(kind=dp), DIMENSION(:, :), INTENT(in), & 2647 OPTIONAL, TARGET :: uv 2648 REAL(kind=dp), INTENT(out), OPTIONAL :: e_k 2649 2650 CHARACTER(len=*), PARAMETER :: routineN = 'pint_calc_e_kin_beads_u', & 2651 routineP = moduleN//':'//routineN 2652 2653 INTEGER :: ib, idim 2654 REAL(kind=dp) :: res 2655 REAL(kind=dp), DIMENSION(:, :), POINTER :: my_uv 2656 2657 CPASSERT(ASSOCIATED(pint_env)) 2658 CPASSERT(pint_env%ref_count > 0) 2659 res = -1.0_dp 2660 my_uv => pint_env%uv 2661 IF (PRESENT(uv)) my_uv => uv 2662 res = 0._dp 2663 DO idim = 1, pint_env%ndim 2664 DO ib = 1, pint_env%p 2665 res = res + pint_env%mass_fict(ib, idim)*my_uv(ib, idim)**2 2666 END DO 2667 END DO 2668 res = res*0.5 2669 IF (.NOT. PRESENT(uv)) pint_env%e_kin_beads = res 2670 IF (PRESENT(e_k)) e_k = res 2671 RETURN 2672 END SUBROUTINE pint_calc_e_kin_beads_u 2673 2674! *************************************************************************** 2675!> \brief Calculate the virial estimator of the real (quantum) kinetic energy 2676!> \param pint_env ... 2677!> \param e_vir ... 2678!> \author hforbert 2679!> \note This subroutine modifies pint_env%energy(e_kin_virial_id) global 2680!> variable [lwalewski] 2681! ************************************************************************************************** 2682 SUBROUTINE pint_calc_e_vir(pint_env, e_vir) 2683 TYPE(pint_env_type), POINTER :: pint_env 2684 REAL(kind=dp), INTENT(out), OPTIONAL :: e_vir 2685 2686 CHARACTER(len=*), PARAMETER :: routineN = 'pint_calc_e_vir', & 2687 routineP = moduleN//':'//routineN 2688 2689 INTEGER :: ib, idim 2690 REAL(kind=dp) :: res, xcentroid 2691 2692 CPASSERT(ASSOCIATED(pint_env)) 2693 CPASSERT(pint_env%ref_count > 0) 2694 res = -1.0_dp 2695 res = 0._dp 2696 DO idim = 1, pint_env%ndim 2697 ! calculate the centroid 2698 xcentroid = 0._dp 2699 DO ib = 1, pint_env%p 2700 xcentroid = xcentroid + pint_env%x(ib, idim) 2701 END DO 2702 xcentroid = xcentroid/REAL(pint_env%p, dp) 2703 DO ib = 1, pint_env%p 2704 res = res + (pint_env%x(ib, idim) - xcentroid)*pint_env%f(ib, idim) 2705 END DO 2706 END DO 2707 res = 0.5_dp*(REAL(pint_env%ndim, dp)* & 2708 (pint_env%kT*pint_env%propagator%temp_sim2phys) - res/REAL(pint_env%p, dp)) 2709 pint_env%energy(e_kin_virial_id) = res 2710 IF (PRESENT(e_vir)) e_vir = res 2711 RETURN 2712 END SUBROUTINE pint_calc_e_vir 2713 2714! *************************************************************************** 2715!> \brief calculates the energy (potential and kinetic) of the Nose-Hoover 2716!> chain thermostats 2717!> \param pint_env the path integral environment 2718!> \author fawzi 2719! ************************************************************************************************** 2720 SUBROUTINE pint_calc_nh_energy(pint_env) 2721 TYPE(pint_env_type), POINTER :: pint_env 2722 2723 CHARACTER(len=*), PARAMETER :: routineN = 'pint_calc_nh_energy', & 2724 routineP = moduleN//':'//routineN 2725 2726 INTEGER :: ib, idim, inos 2727 REAL(kind=dp) :: ekin, epot 2728 2729 CPASSERT(ASSOCIATED(pint_env)) 2730 CPASSERT(pint_env%ref_count > 0) 2731 ekin = 0._dp 2732 DO idim = 1, pint_env%ndim 2733 DO ib = 1, pint_env%p 2734 DO inos = 1, pint_env%nnos 2735 ekin = ekin + pint_env%Q(ib)*pint_env%tv(inos, ib, idim)**2 2736 END DO 2737 END DO 2738 END DO 2739 pint_env%e_kin_t = 0.5_dp*ekin 2740 epot = 0._dp 2741 DO idim = 1, pint_env%ndim 2742 DO ib = 1, pint_env%p 2743 DO inos = 1, pint_env%nnos 2744 epot = epot + pint_env%tx(inos, ib, idim) 2745 END DO 2746 END DO 2747 END DO 2748 pint_env%e_pot_t = pint_env%kT*epot 2749 RETURN 2750 END SUBROUTINE pint_calc_nh_energy 2751 2752! *************************************************************************** 2753!> \brief calculates the total link action of the PI system (excluding helium) 2754!> \param pint_env the path integral environment 2755!> \return ... 2756!> \author Felix Uhl 2757! ************************************************************************************************** 2758 FUNCTION pint_calc_total_link_action(pint_env) RESULT(link_action) 2759 TYPE(pint_env_type), POINTER :: pint_env 2760 REAL(KIND=dp) :: link_action 2761 2762 CHARACTER(len=*), PARAMETER :: routineN = 'pint_calc_total_link_action', & 2763 routineP = moduleN//':'//routineN 2764 2765 INTEGER :: iatom, ibead, idim, indx 2766 REAL(KIND=dp) :: hb2m, tau, tmp_link_action 2767 REAL(KIND=dp), DIMENSION(3) :: r 2768 2769 !tau = 1/(k_B T p) 2770 tau = pint_env%beta/REAL(pint_env%p, dp) 2771 2772 CPASSERT(ASSOCIATED(pint_env)) 2773 CPASSERT(pint_env%ref_count > 0) 2774 link_action = 0.0_dp 2775 DO iatom = 1, pint_env%ndim/3 2776 ! hbar / (2.0*m) 2777 hb2m = 1.0_dp/pint_env%mass((iatom - 1)*3 + 1) 2778 tmp_link_action = 0.0_dp 2779 DO ibead = 1, pint_env%p - 1 2780 DO idim = 1, 3 2781 indx = (iatom - 1)*3 + idim 2782 r(idim) = pint_env%x(ibead, indx) - pint_env%x(ibead + 1, indx) 2783 END DO 2784 tmp_link_action = tmp_link_action + (r(1)*r(1) + r(2)*r(2) + r(3)*r(3)) 2785 END DO 2786 DO idim = 1, 3 2787 indx = (iatom - 1)*3 + idim 2788 r(idim) = pint_env%x(pint_env%p, indx) - pint_env%x(1, indx) 2789 END DO 2790 tmp_link_action = tmp_link_action + (r(1)*r(1) + r(2)*r(2) + r(3)*r(3)) 2791 link_action = link_action + tmp_link_action/hb2m 2792 END DO 2793 2794 link_action = link_action/(2.0_dp*tau) 2795 2796 RETURN 2797 2798 END FUNCTION pint_calc_total_link_action 2799 2800! *************************************************************************** 2801!> \brief calculates the potential action of the PI system (excluding helium) 2802!> \param pint_env the path integral environment 2803!> \return ... 2804!> \author Felix Uhl 2805! ************************************************************************************************** 2806 FUNCTION pint_calc_total_pot_action(pint_env) RESULT(pot_action) 2807 TYPE(pint_env_type), POINTER :: pint_env 2808 REAL(KIND=dp) :: pot_action 2809 2810 CHARACTER(len=*), PARAMETER :: routineN = 'pint_calc_total_pot_action', & 2811 routineP = moduleN//':'//routineN 2812 2813 REAL(KIND=dp) :: tau 2814 2815 CPASSERT(ASSOCIATED(pint_env)) 2816 CPASSERT(pint_env%ref_count > 0) 2817 2818 tau = pint_env%beta/REAL(pint_env%p, dp) 2819 pot_action = tau*SUM(pint_env%e_pot_bead) 2820 2821 RETURN 2822 2823 END FUNCTION pint_calc_total_pot_action 2824 2825! *************************************************************************** 2826!> \brief calculates the total action of the PI system (excluding helium) 2827!> \param pint_env the path integral environment 2828!> \author Felix Uhl 2829! ************************************************************************************************** 2830 SUBROUTINE pint_calc_total_action(pint_env) 2831 TYPE(pint_env_type), POINTER :: pint_env 2832 2833 CHARACTER(len=*), PARAMETER :: routineN = 'pint_calc_total_action', & 2834 routineP = moduleN//':'//routineN 2835 2836 CPASSERT(ASSOCIATED(pint_env)) 2837 CPASSERT(pint_env%ref_count > 0) 2838 2839 pint_env%pot_action = pint_calc_total_pot_action(pint_env) 2840 pint_env%link_action = pint_calc_total_link_action(pint_env) 2841 2842 RETURN 2843 END SUBROUTINE pint_calc_total_action 2844 2845END MODULE pint_methods 2846