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