1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Provides integrator routines (velocity verlet) for all the 8!> ensemble types 9!> \par History 10!> JGH (15-Mar-2001) : Pass logical for box change to force routine 11!> Harald Forbert (Apr-2001): added path integral routine nvt_pimd 12!> CJM (15-Apr-2001) : added coef integrators and energy routines 13!> Joost VandeVondele (Juli-2003): simple version of isokinetic ensemble 14!> Teodoro Laino [tlaino] 10.2007 - University of Zurich: Generalization to 15!> different kind of thermostats 16!> Teodoro Laino [tlaino] 11.2007 - Metadynamics: now part of the MD modules 17!> Marcella Iannuzzi 02.2008 - Collecting common code (VV and creation of 18!> a temporary type) 19!> Teodoro Laino [tlaino] 02.2008 - Splitting integrator module and keeping in 20!> integrator only the INTEGRATORS 21!> Lianheng Tong [LT] 12.2013 - Added regions to Langevin MD 22!> \author CJM 23! ************************************************************************************************** 24MODULE integrator 25 USE atomic_kind_list_types, ONLY: atomic_kind_list_type 26 USE atomic_kind_types, ONLY: atomic_kind_type,& 27 get_atomic_kind,& 28 get_atomic_kind_set 29 USE barostat_types, ONLY: barostat_type 30 USE cell_types, ONLY: cell_type,& 31 init_cell,& 32 parse_cell_line 33 USE constraint, ONLY: rattle_control,& 34 shake_control,& 35 shake_roll_control,& 36 shake_update_targets 37 USE constraint_fxd, ONLY: fix_atom_control 38 USE constraint_util, ONLY: getold,& 39 pv_constraint 40 USE cp_control_types, ONLY: dft_control_type 41 USE cp_log_handling, ONLY: cp_to_string 42 USE cp_para_types, ONLY: cp_para_env_type 43 USE cp_parser_methods, ONLY: parser_get_next_line,& 44 parser_read_line 45 USE cp_subsys_types, ONLY: cp_subsys_get,& 46 cp_subsys_type 47 USE cp_units, ONLY: cp_unit_to_cp2k 48 USE distribution_1d_types, ONLY: distribution_1d_type 49 USE eigenvalueproblems, ONLY: diagonalise 50 USE extended_system_dynamics, ONLY: shell_scale_comv 51 USE extended_system_types, ONLY: npt_info_type 52 USE force_env_methods, ONLY: force_env_calc_energy_force 53 USE force_env_types, ONLY: force_env_get,& 54 force_env_type 55 USE global_types, ONLY: global_environment_type 56 USE input_constants, ONLY: ehrenfest,& 57 npe_f_ensemble,& 58 npe_i_ensemble 59 USE integrator_utils, ONLY: & 60 allocate_old, allocate_tmp, damp_v, damp_veps, deallocate_old, get_s_ds, & 61 old_variables_type, rattle_roll_setup, set, tmp_variables_type, update_dealloc_tmp, & 62 update_pv, update_veps, variable_timestep, vv_first, vv_second 63 USE kinds, ONLY: dp 64 USE mathlib, ONLY: matmul_3x3,& 65 transpose_3d 66 USE md_environment_types, ONLY: get_md_env,& 67 md_environment_type,& 68 set_md_env 69 USE metadynamics, ONLY: metadyn_integrator,& 70 metadyn_velocities_colvar 71 USE molecule_kind_list_types, ONLY: molecule_kind_list_type 72 USE molecule_kind_types, ONLY: molecule_kind_type 73 USE molecule_list_types, ONLY: molecule_list_type 74 USE molecule_types, ONLY: global_constraint_type,& 75 molecule_type 76 USE particle_list_types, ONLY: particle_list_type 77 USE particle_types, ONLY: particle_type,& 78 update_particle_set 79 USE physcon, ONLY: bohr,& 80 femtoseconds 81 USE qmmm_util, ONLY: apply_qmmm_walls_reflective 82 USE qmmmx_update, ONLY: qmmmx_update_force_env 83 USE qs_environment_types, ONLY: get_qs_env 84 USE reftraj_types, ONLY: reftraj_type 85 USE reftraj_util, ONLY: compute_msd_reftraj 86 USE rt_propagation_methods, ONLY: propagation_step 87 USE rt_propagation_output, ONLY: rt_prop_output 88 USE rt_propagation_types, ONLY: rt_prop_type 89 USE shell_opt, ONLY: optimize_shell_core 90 USE simpar_types, ONLY: simpar_type 91 USE thermal_region_types, ONLY: thermal_region_type,& 92 thermal_regions_type 93 USE thermostat_methods, ONLY: apply_thermostat_baro,& 94 apply_thermostat_particles,& 95 apply_thermostat_shells 96 USE thermostat_types, ONLY: thermostat_type 97 USE virial_methods, ONLY: virial_evaluate 98 USE virial_types, ONLY: virial_type 99#include "../base/base_uses.f90" 100 101 IMPLICIT NONE 102 103 PRIVATE 104 105 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'integrator' 106 107 PUBLIC :: isokin, langevin, nve, nvt, npt_i, npt_f, nve_respa 108 PUBLIC :: nph_uniaxial_damped, nph_uniaxial, nvt_adiabatic, reftraj 109 110CONTAINS 111 112! ************************************************************************************************** 113!> \brief Langevin integrator for particle positions & momenta (Brownian dynamics) 114!> \param md_env ... 115!> \par Literature 116!> - A. Ricci and G. Ciccotti, Mol. Phys. 101, 1927-1931 (2003) 117!> - For langevin regions: 118!> - L. Kantorovich, Phys. Rev. B 78, 094304 (2008) 119!> - L. Kantorovich and N. Rompotis, Phys. Rev. B 78, 094305 (2008) 120!> \par History 121!> - Created (01.07.2005,MK) 122!> - Added support for only performing Langevin MD on a region of atoms 123!> (01.12.2013, LT) 124!> \author Matthias Krack 125! ************************************************************************************************** 126 SUBROUTINE langevin(md_env) 127 128 TYPE(md_environment_type), POINTER :: md_env 129 130 CHARACTER(len=*), PARAMETER :: routineN = 'langevin', routineP = moduleN//':'//routineN 131 132 INTEGER :: iparticle, iparticle_kind, iparticle_local, iparticle_reg, ireg, nparticle, & 133 nparticle_kind, nparticle_local, nshell 134 INTEGER, POINTER :: itimes 135 LOGICAL, ALLOCATABLE, DIMENSION(:) :: do_langevin 136 REAL(KIND=dp) :: c, c1, c2, c3, c4, dm, dt, gam, mass, & 137 noisy_gamma_region, reg_temp, sigma 138 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: var_w 139 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: pos, vel, w 140 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds 141 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 142 TYPE(atomic_kind_type), POINTER :: atomic_kind 143 TYPE(cell_type), POINTER :: cell 144 TYPE(cp_para_env_type), POINTER :: para_env 145 TYPE(cp_subsys_type), POINTER :: subsys 146 TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles 147 TYPE(force_env_type), POINTER :: force_env 148 TYPE(global_constraint_type), POINTER :: gci 149 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds 150 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set 151 TYPE(molecule_list_type), POINTER :: molecules 152 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 153 TYPE(particle_list_type), POINTER :: particles 154 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 155 TYPE(simpar_type), POINTER :: simpar 156 TYPE(thermal_region_type), POINTER :: thermal_region 157 TYPE(thermal_regions_type), POINTER :: thermal_regions 158 TYPE(virial_type), POINTER :: virial 159 160 NULLIFY (cell, para_env, gci, force_env) 161 NULLIFY (atomic_kinds, local_particles, subsys, local_molecules, molecule_kinds, molecules) 162 NULLIFY (molecule_kind_set, molecule_set, particles, particle_set, simpar, virial) 163 NULLIFY (thermal_region, thermal_regions, itimes) 164 165 CALL get_md_env(md_env=md_env, simpar=simpar, force_env=force_env, & 166 para_env=para_env, thermal_regions=thermal_regions, & 167 itimes=itimes) 168 169 dt = simpar%dt 170 gam = simpar%gamma + simpar%shadow_gamma 171 nshell = 0 172 173 CALL force_env_get(force_env=force_env, subsys=subsys, cell=cell) 174 175 ! Do some checks on coordinates and box 176 CALL apply_qmmm_walls_reflective(force_env) 177 178 CALL cp_subsys_get(subsys=subsys, & 179 atomic_kinds=atomic_kinds, & 180 gci=gci, & 181 local_particles=local_particles, & 182 local_molecules=local_molecules, & 183 molecules=molecules, & 184 molecule_kinds=molecule_kinds, & 185 nshell=nshell, & 186 particles=particles, & 187 virial=virial) 188 IF (nshell /= 0) & 189 CPABORT("Langevin dynamics is not yet implemented for core-shell models") 190 191 nparticle_kind = atomic_kinds%n_els 192 atomic_kind_set => atomic_kinds%els 193 molecule_kind_set => molecule_kinds%els 194 195 nparticle = particles%n_els 196 particle_set => particles%els 197 molecule_set => molecules%els 198 199 ! Setup the langevin regions information 200 ALLOCATE (do_langevin(nparticle)) 201 IF (simpar%do_thermal_region) THEN 202 DO iparticle = 1, nparticle 203 do_langevin(iparticle) = thermal_regions%do_langevin(iparticle) 204 END DO 205 ELSE 206 do_langevin(1:nparticle) = .TRUE. 207 END IF 208 209 ! Allocate the temperature dependent variance (var_w) of the 210 ! random variable for each atom. It may be different for different 211 ! atoms because of the possibility of Langevin regions, and var_w 212 ! for each region should depend on the temperature defined in the 213 ! region 214 ! RZK explains: sigma is the variance of the Wiener process associated 215 ! with the stochastic term, sigma = m*var_w = m*(2*k_B*T*gamma*dt), 216 ! noisy_gamma adds excessive noise that is not balanced by the damping term 217 ALLOCATE (var_w(nparticle)) 218 var_w(1:nparticle) = simpar%var_w 219 IF (simpar%do_thermal_region) THEN 220 DO ireg = 1, thermal_regions%nregions 221 thermal_region => thermal_regions%thermal_region(ireg) 222 noisy_gamma_region = thermal_region%noisy_gamma_region 223 DO iparticle_reg = 1, thermal_region%npart 224 iparticle = thermal_region%part_index(iparticle_reg) 225 reg_temp = thermal_region%temp_expected 226 var_w(iparticle) = 2.0_dp*reg_temp*simpar%dt*(simpar%gamma + noisy_gamma_region) 227 END DO 228 END DO 229 END IF 230 231 ! Allocate work storage 232 ALLOCATE (pos(3, nparticle)) 233 pos(:, :) = 0.0_dp 234 235 ALLOCATE (vel(3, nparticle)) 236 vel(:, :) = 0.0_dp 237 238 ALLOCATE (w(3, nparticle)) 239 w(:, :) = 0.0_dp 240 241 IF (simpar%constraint) CALL getold(gci, local_molecules, molecule_set, & 242 molecule_kind_set, particle_set, cell) 243 244 ! Generate random variables 245 DO iparticle_kind = 1, nparticle_kind 246 atomic_kind => atomic_kind_set(iparticle_kind) 247 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass) 248 nparticle_local = local_particles%n_el(iparticle_kind) 249 DO iparticle_local = 1, nparticle_local 250 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) 251 IF (do_langevin(iparticle)) THEN 252 sigma = var_w(iparticle)*mass 253 ASSOCIATE (rng_stream=>local_particles%local_particle_set(iparticle_kind)% & 254 rng(iparticle_local)) 255 w(1, iparticle) = rng_stream%stream%next(variance=sigma) 256 w(2, iparticle) = rng_stream%stream%next(variance=sigma) 257 w(3, iparticle) = rng_stream%stream%next(variance=sigma) 258 END ASSOCIATE 259 END IF 260 END DO 261 END DO 262 263 DEALLOCATE (var_w) 264 265 ! Apply fix atom constraint 266 CALL fix_atom_control(force_env, w) 267 268 ! Velocity Verlet (first part) 269 c = EXP(-0.25_dp*dt*gam) 270 c2 = c*c 271 c4 = c2*c2 272 c1 = dt*c2 273 274 DO iparticle_kind = 1, nparticle_kind 275 atomic_kind => atomic_kind_set(iparticle_kind) 276 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass) 277 nparticle_local = local_particles%n_el(iparticle_kind) 278 dm = 0.5_dp*dt/mass 279 c3 = dm/c2 280 DO iparticle_local = 1, nparticle_local 281 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) 282 IF (do_langevin(iparticle)) THEN 283 vel(:, iparticle) = particle_set(iparticle)%v(:) + & 284 c3*particle_set(iparticle)%f(:) 285 pos(:, iparticle) = particle_set(iparticle)%r(:) + & 286 c1*particle_set(iparticle)%v(:) + & 287 c*dm*(dt*particle_set(iparticle)%f(:) + & 288 w(:, iparticle)) 289 ELSE 290 vel(:, iparticle) = particle_set(iparticle)%v(:) + & 291 dm*particle_set(iparticle)%f(:) 292 pos(:, iparticle) = particle_set(iparticle)%r(:) + & 293 dt*particle_set(iparticle)%v(:) + & 294 dm*dt*particle_set(iparticle)%f(:) 295 END IF 296 END DO 297 END DO 298 299 IF (simpar%constraint) THEN 300 ! Possibly update the target values 301 CALL shake_update_targets(gci, local_molecules, molecule_set, & 302 molecule_kind_set, dt, force_env%root_section) 303 304 CALL shake_control(gci, local_molecules, molecule_set, molecule_kind_set, & 305 particle_set, pos, vel, dt, simpar%shake_tol, & 306 simpar%info_constraint, simpar%lagrange_multipliers, & 307 simpar%dump_lm, cell, para_env%group, local_particles) 308 END IF 309 310 ! Broadcast the new particle positions 311 CALL update_particle_set(particle_set, para_env%group, pos=pos) 312 313 DEALLOCATE (pos) 314 315 ! Update forces 316 CALL force_env_calc_energy_force(force_env) 317 318 ! Metadynamics 319 CALL metadyn_integrator(force_env, itimes, vel) 320 321 ! Update Verlet (second part) 322 DO iparticle_kind = 1, nparticle_kind 323 atomic_kind => atomic_kind_set(iparticle_kind) 324 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass) 325 dm = 0.5_dp*dt/mass 326 c3 = dm/c2 327 nparticle_local = local_particles%n_el(iparticle_kind) 328 DO iparticle_local = 1, nparticle_local 329 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) 330 IF (do_langevin(iparticle)) THEN 331 vel(1, iparticle) = vel(1, iparticle) + c3*particle_set(iparticle)%f(1) 332 vel(2, iparticle) = vel(2, iparticle) + c3*particle_set(iparticle)%f(2) 333 vel(3, iparticle) = vel(3, iparticle) + c3*particle_set(iparticle)%f(3) 334 vel(1, iparticle) = c4*vel(1, iparticle) + c2*w(1, iparticle)/mass 335 vel(2, iparticle) = c4*vel(2, iparticle) + c2*w(2, iparticle)/mass 336 vel(3, iparticle) = c4*vel(3, iparticle) + c2*w(3, iparticle)/mass 337 ELSE 338 vel(1, iparticle) = vel(1, iparticle) + dm*particle_set(iparticle)%f(1) 339 vel(2, iparticle) = vel(2, iparticle) + dm*particle_set(iparticle)%f(2) 340 vel(3, iparticle) = vel(3, iparticle) + dm*particle_set(iparticle)%f(3) 341 END IF 342 END DO 343 END DO 344 345 IF (simpar%temperature_annealing) THEN 346 simpar%temp_ext = simpar%temp_ext*simpar%f_temperature_annealing 347 simpar%var_w = simpar%var_w*simpar%f_temperature_annealing 348 END IF 349 350 IF (simpar%constraint) THEN 351 CALL rattle_control(gci, local_molecules, molecule_set, molecule_kind_set, & 352 particle_set, vel, dt, simpar%shake_tol, & 353 simpar%info_constraint, simpar%lagrange_multipliers, & 354 simpar%dump_lm, cell, para_env%group, local_particles) 355 END IF 356 357 ! Broadcast the new particle velocities 358 CALL update_particle_set(particle_set, para_env%group, vel=vel) 359 360 DEALLOCATE (vel) 361 362 DEALLOCATE (w) 363 364 DEALLOCATE (do_langevin) 365 366 ! Update virial 367 IF (simpar%constraint) CALL pv_constraint(gci, local_molecules, molecule_set, & 368 molecule_kind_set, particle_set, virial, para_env%group) 369 370 CALL virial_evaluate(atomic_kind_set, particle_set, local_particles, & 371 virial, para_env%group) 372 373 END SUBROUTINE langevin 374 375! ************************************************************************************************** 376!> \brief nve integrator for particle positions & momenta 377!> \param md_env ... 378!> \param globenv ... 379!> \par History 380!> - the local particle lists are used instead of pnode (Sep. 2003,MK) 381!> - usage of fragments retrieved from the force environment (Oct. 2003,MK) 382!> \author CJM 383! ************************************************************************************************** 384 SUBROUTINE nve(md_env, globenv) 385 386 TYPE(md_environment_type), POINTER :: md_env 387 TYPE(global_environment_type), POINTER :: globenv 388 389 CHARACTER(len=*), PARAMETER :: routineN = 'nve', routineP = moduleN//':'//routineN 390 391 INTEGER :: i_iter, n_iter, nparticle, & 392 nparticle_kind, nshell 393 INTEGER, POINTER :: itimes 394 LOGICAL :: deallocate_vel, ehrenfest_md, & 395 shell_adiabatic, shell_check_distance, & 396 shell_present 397 REAL(KIND=dp) :: dt 398 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: v_old 399 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds 400 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 401 TYPE(cell_type), POINTER :: cell 402 TYPE(cp_para_env_type), POINTER :: para_env 403 TYPE(cp_subsys_type), POINTER :: subsys 404 TYPE(dft_control_type), POINTER :: dft_control 405 TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles 406 TYPE(force_env_type), POINTER :: force_env 407 TYPE(global_constraint_type), POINTER :: gci 408 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds 409 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set 410 TYPE(molecule_list_type), POINTER :: molecules 411 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 412 TYPE(particle_list_type), POINTER :: core_particles, particles, & 413 shell_particles 414 TYPE(particle_type), DIMENSION(:), POINTER :: core_particle_set, particle_set, & 415 shell_particle_set 416 TYPE(rt_prop_type), POINTER :: rtp 417 TYPE(simpar_type), POINTER :: simpar 418 TYPE(thermostat_type), POINTER :: thermostat_coeff, thermostat_shell 419 TYPE(tmp_variables_type), POINTER :: tmp 420 TYPE(virial_type), POINTER :: virial 421 422 NULLIFY (thermostat_coeff, tmp) 423 NULLIFY (subsys, simpar, para_env, cell, gci, force_env, virial) 424 NULLIFY (atomic_kinds, local_particles, molecules, molecule_kind_set, molecule_set, particle_set) 425 NULLIFY (shell_particles, shell_particle_set, core_particles, & 426 core_particle_set, thermostat_shell, dft_control, itimes) 427 CALL get_md_env(md_env=md_env, simpar=simpar, force_env=force_env, & 428 thermostat_coeff=thermostat_coeff, thermostat_shell=thermostat_shell, & 429 para_env=para_env, ehrenfest_md=ehrenfest_md, itimes=itimes) 430 dt = simpar%dt 431 CALL force_env_get(force_env=force_env, subsys=subsys, cell=cell) 432 433 ! Do some checks on coordinates and box 434 CALL apply_qmmm_walls_reflective(force_env) 435 436 CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, & 437 particles=particles, local_molecules=local_molecules, molecules=molecules, & 438 molecule_kinds=molecule_kinds, gci=gci, virial=virial) 439 440 nparticle_kind = atomic_kinds%n_els 441 atomic_kind_set => atomic_kinds%els 442 molecule_kind_set => molecule_kinds%els 443 444 nparticle = particles%n_els 445 particle_set => particles%els 446 molecule_set => molecules%els 447 448 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, & 449 shell_present=shell_present, shell_adiabatic=shell_adiabatic, & 450 shell_check_distance=shell_check_distance) 451 452 IF (shell_present) THEN 453 CALL cp_subsys_get(subsys=subsys, shell_particles=shell_particles, & 454 core_particles=core_particles) 455 shell_particle_set => shell_particles%els 456 nshell = SIZE(shell_particles%els) 457 458 IF (shell_adiabatic) THEN 459 core_particle_set => core_particles%els 460 END IF 461 END IF 462 463 CALL allocate_tmp(md_env, tmp, nparticle, nshell, shell_adiabatic) 464 465 ! Apply thermostat over the full set of shells if required 466 IF (shell_adiabatic) THEN 467 CALL apply_thermostat_shells(thermostat_shell, atomic_kind_set, particle_set, & 468 local_particles, para_env%group, shell_particle_set=shell_particle_set, & 469 core_particle_set=core_particle_set) 470 END IF 471 472 IF (simpar%constraint) CALL getold(gci, local_molecules, molecule_set, & 473 molecule_kind_set, particle_set, cell) 474 475 ! Velocity Verlet (first part) 476 CALL vv_first(tmp, atomic_kind_set, local_particles, particle_set, & 477 core_particle_set, shell_particle_set, nparticle_kind, shell_adiabatic, dt) 478 479 IF (simpar%variable_dt) CALL variable_timestep(md_env, tmp, dt, simpar, para_env, atomic_kind_set, & 480 local_particles, particle_set, core_particle_set, shell_particle_set, & 481 nparticle_kind, shell_adiabatic) 482 483 IF (simpar%constraint) THEN 484 ! Possibly update the target values 485 CALL shake_update_targets(gci, local_molecules, molecule_set, & 486 molecule_kind_set, dt, force_env%root_section) 487 488 CALL shake_control(gci, local_molecules, molecule_set, & 489 molecule_kind_set, particle_set, tmp%pos, tmp%vel, dt, simpar%shake_tol, & 490 simpar%info_constraint, simpar%lagrange_multipliers, simpar%dump_lm, & 491 cell, para_env%group, local_particles) 492 END IF 493 494 ! Broadcast the new particle positions and deallocate pos part of temporary 495 CALL update_dealloc_tmp(tmp, particle_set, shell_particle_set, & 496 core_particle_set, para_env, shell_adiabatic, pos=.TRUE.) 497 498 IF (shell_adiabatic .AND. shell_check_distance) THEN 499 CALL optimize_shell_core(force_env, particle_set, & 500 shell_particle_set, core_particle_set, globenv, tmp=tmp, check=.TRUE.) 501 END IF 502 503 ! Update forces 504 ! In case of ehrenfest dynamics, velocities need to be iterated 505 IF (ehrenfest_md) THEN 506 ALLOCATE (v_old(3, SIZE(tmp%vel, 2))) 507 v_old(:, :) = tmp%vel 508 CALL vv_second(tmp, atomic_kind_set, local_particles, particle_set, & 509 core_particle_set, shell_particle_set, nparticle_kind, shell_adiabatic, dt) 510 CALL update_dealloc_tmp(tmp, particle_set, shell_particle_set, & 511 core_particle_set, para_env, shell_adiabatic, vel=.TRUE., & 512 should_deall_vel=.FALSE.) 513 tmp%vel = v_old 514 CALL get_qs_env(force_env%qs_env, dft_control=dft_control) 515 n_iter = dft_control%rtp_control%max_iter 516 ELSE 517 n_iter = 1 518 END IF 519 520 DO i_iter = 1, n_iter 521 522 IF (ehrenfest_md) THEN 523 CALL get_qs_env(qs_env=force_env%qs_env, rtp=rtp) 524 rtp%iter = i_iter 525 tmp%vel = v_old 526 CALL propagation_step(force_env%qs_env, rtp, dft_control%rtp_control) 527 END IF 528 529 ![NB] let nve work with force mixing which does not have consistent energies and forces 530 CALL force_env_calc_energy_force(force_env, require_consistent_energy_force=.FALSE.) 531 532 IF (ehrenfest_md) THEN 533 CALL rt_prop_output(force_env%qs_env, ehrenfest, delta_iter=force_env%qs_env%rtp%delta_iter) 534 ENDIF 535 536 ! Metadynamics 537 CALL metadyn_integrator(force_env, itimes, tmp%vel) 538 539 ! Velocity Verlet (second part) 540 CALL vv_second(tmp, atomic_kind_set, local_particles, particle_set, & 541 core_particle_set, shell_particle_set, nparticle_kind, shell_adiabatic, dt) 542 543 IF (simpar%constraint) CALL rattle_control(gci, local_molecules, molecule_set, & 544 molecule_kind_set, particle_set, tmp%vel, dt, simpar%shake_tol, & 545 simpar%info_constraint, simpar%lagrange_multipliers, simpar%dump_lm, & 546 cell, para_env%group, local_particles) 547 548 ! Apply thermostat over the full set of shell if required 549 IF (shell_adiabatic) THEN 550 CALL apply_thermostat_shells(thermostat_shell, atomic_kind_set, particle_set, & 551 local_particles, para_env%group, vel=tmp%vel, & 552 shell_vel=tmp%shell_vel, core_vel=tmp%core_vel) 553 END IF 554 555 IF (simpar%annealing) THEN 556 tmp%vel(:, :) = tmp%vel(:, :)*simpar%f_annealing 557 IF (shell_adiabatic) THEN 558 CALL shell_scale_comv(atomic_kind_set, local_particles, particle_set, & 559 tmp%vel, tmp%shell_vel, tmp%core_vel) 560 END IF 561 END IF 562 563 IF (ehrenfest_md) deallocate_vel = force_env%qs_env%rtp%converged 564 IF (i_iter .EQ. n_iter) deallocate_vel = .TRUE. 565 ! Broadcast the new particle velocities and deallocate the full temporary 566 CALL update_dealloc_tmp(tmp, particle_set, shell_particle_set, & 567 core_particle_set, para_env, shell_adiabatic, vel=.TRUE., & 568 should_deall_vel=deallocate_vel) 569 IF (ehrenfest_md) THEN 570 IF (force_env%qs_env%rtp%converged) EXIT 571 END IF 572 573 END DO 574 575 ! Update virial 576 IF (simpar%constraint) CALL pv_constraint(gci, local_molecules, & 577 molecule_set, molecule_kind_set, particle_set, virial, para_env%group) 578 579 CALL virial_evaluate(atomic_kind_set, particle_set, & 580 local_particles, virial, para_env%group) 581 582 END SUBROUTINE nve 583 584! ************************************************************************************************** 585!> \brief simplest version of the isokinetic gaussian thermostat 586!> \param md_env ... 587!> \par History 588!> - Created [2004-07] 589!> \author Joost VandeVondele 590!> \note 591!> - time reversible, and conserves the kinetic energy to machine precision 592!> - is not yet supposed to work for e.g. constraints, our the extended version 593!> of this thermostat 594!> see: 595!> - Zhang F. , JCP 106, 6102 (1997) 596!> - Minary P. et al, JCP 118, 2510 (2003) 597! ************************************************************************************************** 598 SUBROUTINE isokin(md_env) 599 600 TYPE(md_environment_type), POINTER :: md_env 601 602 CHARACTER(len=*), PARAMETER :: routineN = 'isokin', routineP = moduleN//':'//routineN 603 604 INTEGER :: nparticle, nparticle_kind, nshell 605 INTEGER, POINTER :: itimes 606 LOGICAL :: shell_adiabatic, shell_present 607 REAL(KIND=dp) :: dt 608 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds 609 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 610 TYPE(cp_para_env_type), POINTER :: para_env 611 TYPE(cp_subsys_type), POINTER :: subsys 612 TYPE(distribution_1d_type), POINTER :: local_particles 613 TYPE(force_env_type), POINTER :: force_env 614 TYPE(particle_list_type), POINTER :: core_particles, particles, & 615 shell_particles 616 TYPE(particle_type), DIMENSION(:), POINTER :: core_particle_set, particle_set, & 617 shell_particle_set 618 TYPE(simpar_type), POINTER :: simpar 619 TYPE(tmp_variables_type), POINTER :: tmp 620 621 NULLIFY (force_env, tmp, simpar, itimes) 622 NULLIFY (atomic_kinds, para_env, subsys, local_particles) 623 NULLIFY (core_particles, particles, shell_particles) 624 NULLIFY (core_particle_set, particle_set, shell_particle_set) 625 626 CALL get_md_env(md_env=md_env, simpar=simpar, force_env=force_env, & 627 para_env=para_env, itimes=itimes) 628 629 dt = simpar%dt 630 631 CALL force_env_get(force_env=force_env, subsys=subsys) 632 633 ! Do some checks on coordinates and box 634 CALL apply_qmmm_walls_reflective(force_env) 635 636 IF (simpar%constraint) THEN 637 CPABORT("Constraints not yet implemented") 638 END IF 639 640 CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds, & 641 local_particles=local_particles, & 642 particles=particles) 643 644 nparticle_kind = atomic_kinds%n_els 645 atomic_kind_set => atomic_kinds%els 646 nparticle = particles%n_els 647 particle_set => particles%els 648 649 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, & 650 shell_present=shell_present, shell_adiabatic=shell_adiabatic) 651 652 IF (shell_present) THEN 653 CALL cp_subsys_get(subsys=subsys, shell_particles=shell_particles, & 654 core_particles=core_particles) 655 shell_particle_set => shell_particles%els 656 nshell = SIZE(shell_particles%els) 657 658 IF (shell_adiabatic) THEN 659 core_particle_set => core_particles%els 660 END IF 661 END IF 662 663 CALL allocate_tmp(md_env, tmp, nparticle, nshell, shell_adiabatic) 664 665 ! compute s,ds 666 CALL get_s_ds(tmp, nparticle_kind, atomic_kind_set, local_particles, particle_set, & 667 dt, para_env) 668 669 ! Velocity Verlet (first part) 670 tmp%scale_v(1:3) = SQRT(1.0_dp/tmp%ds) 671 tmp%poly_v(1:3) = 2.0_dp*tmp%s/SQRT(tmp%ds)/dt 672 CALL vv_first(tmp, atomic_kind_set, local_particles, particle_set, & 673 core_particle_set, shell_particle_set, nparticle_kind, & 674 shell_adiabatic, dt) 675 676 IF (simpar%variable_dt) CALL variable_timestep(md_env, tmp, dt, simpar, para_env, atomic_kind_set, & 677 local_particles, particle_set, core_particle_set, shell_particle_set, & 678 nparticle_kind, shell_adiabatic) 679 680 ! Broadcast the new particle positions and deallocate the pos components of temporary 681 CALL update_dealloc_tmp(tmp, particle_set, shell_particle_set, & 682 core_particle_set, para_env, shell_adiabatic, pos=.TRUE.) 683 684 CALL force_env_calc_energy_force(force_env) 685 686 ! Metadynamics 687 CALL metadyn_integrator(force_env, itimes, tmp%vel) 688 689 ! compute s,ds 690 CALL get_s_ds(tmp, nparticle_kind, atomic_kind_set, local_particles, particle_set, & 691 dt, para_env, tmpv=.TRUE.) 692 693 ! Velocity Verlet (second part) 694 tmp%scale_v(1:3) = SQRT(1.0_dp/tmp%ds) 695 tmp%poly_v(1:3) = 2.0_dp*tmp%s/SQRT(tmp%ds)/dt 696 CALL vv_second(tmp, atomic_kind_set, local_particles, particle_set, & 697 core_particle_set, shell_particle_set, nparticle_kind, & 698 shell_adiabatic, dt) 699 700 IF (simpar%annealing) tmp%vel(:, :) = tmp%vel(:, :)*simpar%f_annealing 701 702 ! Broadcast the new particle velocities and deallocate the temporary 703 CALL update_dealloc_tmp(tmp, particle_set, shell_particle_set, & 704 core_particle_set, para_env, shell_adiabatic, vel=.TRUE.) 705 706 END SUBROUTINE isokin 707! ************************************************************************************************** 708!> \brief nvt adiabatic integrator for particle positions & momenta 709!> \param md_env ... 710!> \param globenv ... 711!> \par History 712!> - the local particle lists are used instead of pnode (Sep. 2003,MK) 713!> - usage of fragments retrieved from the force environment (Oct. 2003,MK) 714!> \author CJM 715! ************************************************************************************************** 716 SUBROUTINE nvt_adiabatic(md_env, globenv) 717 718 TYPE(md_environment_type), POINTER :: md_env 719 TYPE(global_environment_type), POINTER :: globenv 720 721 CHARACTER(len=*), PARAMETER :: routineN = 'nvt_adiabatic', routineP = moduleN//':'//routineN 722 723 INTEGER :: ivar, nparticle, nparticle_kind, nshell 724 INTEGER, POINTER :: itimes 725 LOGICAL :: shell_adiabatic, shell_check_distance, & 726 shell_present 727 REAL(KIND=dp) :: dt 728 REAL(KIND=dp), DIMENSION(:), POINTER :: rand 729 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds 730 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 731 TYPE(cell_type), POINTER :: cell 732 TYPE(cp_para_env_type), POINTER :: para_env 733 TYPE(cp_subsys_type), POINTER :: subsys 734 TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles 735 TYPE(force_env_type), POINTER :: force_env 736 TYPE(global_constraint_type), POINTER :: gci 737 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds 738 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set 739 TYPE(molecule_list_type), POINTER :: molecules 740 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 741 TYPE(particle_list_type), POINTER :: core_particles, particles, & 742 shell_particles 743 TYPE(particle_type), DIMENSION(:), POINTER :: core_particle_set, particle_set, & 744 shell_particle_set 745 TYPE(simpar_type), POINTER :: simpar 746 TYPE(thermostat_type), POINTER :: thermostat_coeff, thermostat_fast, & 747 thermostat_shell, thermostat_slow 748 TYPE(tmp_variables_type), POINTER :: tmp 749 TYPE(virial_type), POINTER :: virial 750 751 NULLIFY (gci, force_env, thermostat_coeff, tmp, & 752 thermostat_fast, thermostat_slow, thermostat_shell, cell, shell_particles, & 753 shell_particle_set, core_particles, core_particle_set, rand) 754 NULLIFY (para_env, subsys, local_molecules, local_particles, molecule_kinds, & 755 molecules, molecule_kind_set, molecule_set, atomic_kinds, particles) 756 NULLIFY (simpar, itimes) 757 758 CALL get_md_env(md_env=md_env, simpar=simpar, force_env=force_env, & 759 thermostat_fast=thermostat_fast, thermostat_slow=thermostat_slow, & 760 thermostat_coeff=thermostat_coeff, thermostat_shell=thermostat_shell, & 761 para_env=para_env, itimes=itimes) 762 dt = simpar%dt 763 764 CALL force_env_get(force_env=force_env, subsys=subsys, cell=cell) 765 766 ! Do some checks on coordinates and box 767 CALL apply_qmmm_walls_reflective(force_env) 768 769 CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, & 770 particles=particles, local_molecules=local_molecules, molecules=molecules, & 771 molecule_kinds=molecule_kinds, gci=gci, virial=virial) 772 773 nparticle_kind = atomic_kinds%n_els 774 atomic_kind_set => atomic_kinds%els 775 molecule_kind_set => molecule_kinds%els 776 777 nparticle = particles%n_els 778 particle_set => particles%els 779 molecule_set => molecules%els 780 781 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, & 782 shell_present=shell_present, shell_adiabatic=shell_adiabatic, & 783 shell_check_distance=shell_check_distance) 784 785 IF (ASSOCIATED(force_env%meta_env)) THEN 786 ! Allocate random number for Langevin Thermostat acting on COLVARS 787 IF (force_env%meta_env%langevin) THEN 788 ALLOCATE (rand(force_env%meta_env%n_colvar)) 789 rand(:) = 0.0_dp 790 ENDIF 791 ENDIF 792 793 ! Allocate work storage for positions and velocities 794 IF (shell_present) THEN 795 CALL cp_subsys_get(subsys=subsys, shell_particles=shell_particles, & 796 core_particles=core_particles) 797 shell_particle_set => shell_particles%els 798 nshell = SIZE(shell_particles%els) 799 800 IF (shell_adiabatic) THEN 801 core_particle_set => core_particles%els 802 END IF 803 END IF 804 805 CALL allocate_tmp(md_env, tmp, nparticle, nshell, shell_adiabatic) 806 807 ! Apply Thermostat over the full set of particles 808 IF (shell_adiabatic) THEN 809! CALL apply_thermostat_particles(thermostat_part, molecule_kind_set, molecule_set,& 810! particle_set, local_molecules, para_env%group, shell_adiabatic=shell_adiabatic,& 811! shell_particle_set=shell_particle_set, core_particle_set=core_particle_set) 812 813 CALL apply_thermostat_shells(thermostat_shell, atomic_kind_set, particle_set, & 814 local_particles, para_env%group, shell_particle_set=shell_particle_set, & 815 core_particle_set=core_particle_set) 816 ELSE 817 CALL apply_thermostat_particles(thermostat_fast, force_env, molecule_kind_set, molecule_set, & 818 particle_set, local_molecules, local_particles, para_env%group) 819 820 CALL apply_thermostat_particles(thermostat_slow, force_env, molecule_kind_set, molecule_set, & 821 particle_set, local_molecules, local_particles, para_env%group) 822 END IF 823 824 IF (simpar%constraint) CALL getold(gci, local_molecules, molecule_set, & 825 molecule_kind_set, particle_set, cell) 826 827 ! *** Velocity Verlet for Langeving *** v(t)--> v(t+1/2) 828 IF (ASSOCIATED(force_env%meta_env)) THEN 829 IF (force_env%meta_env%langevin) THEN 830 DO ivar = 1, force_env%meta_env%n_colvar 831 rand(ivar) = force_env%meta_env%rng(ivar)%next() 832 ENDDO 833 CALL metadyn_velocities_colvar(force_env, rand) 834 ENDIF 835 ENDIF 836 837 ! Velocity Verlet (first part) 838 CALL vv_first(tmp, atomic_kind_set, local_particles, particle_set, & 839 core_particle_set, shell_particle_set, nparticle_kind, shell_adiabatic, dt) 840 841 IF (simpar%variable_dt) CALL variable_timestep(md_env, tmp, dt, simpar, para_env, atomic_kind_set, & 842 local_particles, particle_set, core_particle_set, shell_particle_set, & 843 nparticle_kind, shell_adiabatic) 844 845 IF (simpar%constraint) THEN 846 ! Possibly update the target values 847 CALL shake_update_targets(gci, local_molecules, molecule_set, & 848 molecule_kind_set, dt, force_env%root_section) 849 850 CALL shake_control(gci, local_molecules, molecule_set, & 851 molecule_kind_set, particle_set, tmp%pos, tmp%vel, dt, simpar%shake_tol, & 852 simpar%info_constraint, simpar%lagrange_multipliers, simpar%dump_lm, & 853 cell, para_env%group, local_particles) 854 END IF 855 856 ! Broadcast the new particle positions and deallocate pos components of temporary 857 CALL update_dealloc_tmp(tmp, particle_set, shell_particle_set, & 858 core_particle_set, para_env, shell_adiabatic, pos=.TRUE.) 859 860 IF (shell_adiabatic .AND. shell_check_distance) THEN 861 CALL optimize_shell_core(force_env, particle_set, & 862 shell_particle_set, core_particle_set, globenv, tmp=tmp, check=.TRUE.) 863 END IF 864 865 ! Update forces 866 CALL force_env_calc_energy_force(force_env) 867 868 ! Metadynamics 869 CALL metadyn_integrator(force_env, itimes, tmp%vel, rand=rand) 870 871 ! Velocity Verlet (second part) 872 CALL vv_second(tmp, atomic_kind_set, local_particles, particle_set, & 873 core_particle_set, shell_particle_set, nparticle_kind, shell_adiabatic, dt) 874 875 IF (simpar%constraint) CALL rattle_control(gci, local_molecules, molecule_set, & 876 molecule_kind_set, particle_set, tmp%vel, dt, simpar%shake_tol, & 877 simpar%info_constraint, simpar%lagrange_multipliers, simpar%dump_lm, & 878 cell, para_env%group, local_particles) 879 880 ! Apply Thermostat over the full set of particles 881 IF (shell_adiabatic) THEN 882 ! CALL apply_thermostat_particles(thermostat_part,molecule_kind_set, molecule_set, & 883 ! particle_set, local_molecules, para_env%group, shell_adiabatic=shell_adiabatic,& 884 ! vel= tmp%vel, shell_vel= tmp%shell_vel, core_vel= tmp%core_vel) 885 886 CALL apply_thermostat_shells(thermostat_shell, atomic_kind_set, particle_set, & 887 local_particles, para_env%group, vel=tmp%vel, shell_vel=tmp%shell_vel, & 888 core_vel=tmp%core_vel) 889 ELSE 890 CALL apply_thermostat_particles(thermostat_slow, force_env, molecule_kind_set, molecule_set, & 891 particle_set, local_molecules, local_particles, para_env%group, vel=tmp%vel) 892 893 CALL apply_thermostat_particles(thermostat_fast, force_env, molecule_kind_set, molecule_set, & 894 particle_set, local_molecules, local_particles, para_env%group, vel=tmp%vel) 895 END IF 896 897 ! Broadcast the new particle velocities and deallocate temporary 898 CALL update_dealloc_tmp(tmp, particle_set, shell_particle_set, & 899 core_particle_set, para_env, shell_adiabatic, vel=.TRUE.) 900 901 IF (ASSOCIATED(force_env%meta_env)) THEN 902 IF (force_env%meta_env%langevin) THEN 903 DEALLOCATE (rand) 904 ENDIF 905 ENDIF 906 907 ! Update constraint virial 908 IF (simpar%constraint) CALL pv_constraint(gci, local_molecules, & 909 molecule_set, molecule_kind_set, particle_set, virial, para_env%group) 910 911 ! ** Evaluate Virial 912 CALL virial_evaluate(atomic_kind_set, particle_set, & 913 local_particles, virial, para_env%group) 914 915 END SUBROUTINE nvt_adiabatic 916 917! ************************************************************************************************** 918!> \brief nvt integrator for particle positions & momenta 919!> \param md_env ... 920!> \param globenv ... 921!> \par History 922!> - the local particle lists are used instead of pnode (Sep. 2003,MK) 923!> - usage of fragments retrieved from the force environment (Oct. 2003,MK) 924!> \author CJM 925! ************************************************************************************************** 926 SUBROUTINE nvt(md_env, globenv) 927 928 TYPE(md_environment_type), POINTER :: md_env 929 TYPE(global_environment_type), POINTER :: globenv 930 931 CHARACTER(len=*), PARAMETER :: routineN = 'nvt', routineP = moduleN//':'//routineN 932 933 INTEGER :: ivar, nparticle, nparticle_kind, nshell 934 INTEGER, POINTER :: itimes 935 LOGICAL :: shell_adiabatic, shell_check_distance, & 936 shell_present 937 REAL(KIND=dp) :: dt 938 REAL(KIND=dp), DIMENSION(:), POINTER :: rand 939 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds 940 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 941 TYPE(cell_type), POINTER :: cell 942 TYPE(cp_para_env_type), POINTER :: para_env 943 TYPE(cp_subsys_type), POINTER :: subsys 944 TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles 945 TYPE(force_env_type), POINTER :: force_env 946 TYPE(global_constraint_type), POINTER :: gci 947 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds 948 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set 949 TYPE(molecule_list_type), POINTER :: molecules 950 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 951 TYPE(particle_list_type), POINTER :: core_particles, particles, & 952 shell_particles 953 TYPE(particle_type), DIMENSION(:), POINTER :: core_particle_set, particle_set, & 954 shell_particle_set 955 TYPE(simpar_type), POINTER :: simpar 956 TYPE(thermostat_type), POINTER :: thermostat_coeff, thermostat_part, & 957 thermostat_shell 958 TYPE(tmp_variables_type), POINTER :: tmp 959 TYPE(virial_type), POINTER :: virial 960 961 NULLIFY (gci, force_env, thermostat_coeff, tmp, & 962 thermostat_part, thermostat_shell, cell, shell_particles, & 963 shell_particle_set, core_particles, core_particle_set, rand) 964 NULLIFY (para_env, subsys, local_molecules, local_particles, molecule_kinds, & 965 molecules, molecule_kind_set, molecule_set, atomic_kinds, particles) 966 NULLIFY (simpar, thermostat_coeff, thermostat_part, thermostat_shell, itimes) 967 968 CALL get_md_env(md_env=md_env, simpar=simpar, force_env=force_env, & 969 thermostat_part=thermostat_part, thermostat_coeff=thermostat_coeff, & 970 thermostat_shell=thermostat_shell, para_env=para_env, & 971 itimes=itimes) 972 dt = simpar%dt 973 974 CALL force_env_get(force_env=force_env, subsys=subsys, cell=cell) 975 976 ! Do some checks on coordinates and box 977 CALL apply_qmmm_walls_reflective(force_env) 978 979 CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, & 980 particles=particles, local_molecules=local_molecules, molecules=molecules, & 981 molecule_kinds=molecule_kinds, gci=gci, virial=virial) 982 983 nparticle_kind = atomic_kinds%n_els 984 atomic_kind_set => atomic_kinds%els 985 molecule_kind_set => molecule_kinds%els 986 987 nparticle = particles%n_els 988 particle_set => particles%els 989 molecule_set => molecules%els 990 991 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, & 992 shell_present=shell_present, shell_adiabatic=shell_adiabatic, & 993 shell_check_distance=shell_check_distance) 994 995 IF (ASSOCIATED(force_env%meta_env)) THEN 996 ! Allocate random number for Langevin Thermostat acting on COLVARS 997 IF (force_env%meta_env%langevin) THEN 998 ALLOCATE (rand(force_env%meta_env%n_colvar)) 999 rand(:) = 0.0_dp 1000 ENDIF 1001 ENDIF 1002 1003 ! Allocate work storage for positions and velocities 1004 IF (shell_present) THEN 1005 CALL cp_subsys_get(subsys=subsys, shell_particles=shell_particles, & 1006 core_particles=core_particles) 1007 shell_particle_set => shell_particles%els 1008 nshell = SIZE(shell_particles%els) 1009 1010 IF (shell_adiabatic) THEN 1011 core_particle_set => core_particles%els 1012 END IF 1013 END IF 1014 1015 CALL allocate_tmp(md_env, tmp, nparticle, nshell, shell_adiabatic) 1016 1017 ! Apply Thermostat over the full set of particles 1018 IF (shell_adiabatic) THEN 1019 CALL apply_thermostat_particles(thermostat_part, force_env, molecule_kind_set, molecule_set, & 1020 particle_set, local_molecules, local_particles, para_env%group, shell_adiabatic=shell_adiabatic, & 1021 shell_particle_set=shell_particle_set, core_particle_set=core_particle_set) 1022 1023 CALL apply_thermostat_shells(thermostat_shell, atomic_kind_set, particle_set, & 1024 local_particles, para_env%group, shell_particle_set=shell_particle_set, & 1025 core_particle_set=core_particle_set) 1026 ELSE 1027 CALL apply_thermostat_particles(thermostat_part, force_env, molecule_kind_set, molecule_set, & 1028 particle_set, local_molecules, local_particles, para_env%group) 1029 END IF 1030 1031 IF (simpar%constraint) CALL getold(gci, local_molecules, molecule_set, & 1032 molecule_kind_set, particle_set, cell) 1033 1034 ! *** Velocity Verlet for Langeving *** v(t)--> v(t+1/2) 1035 IF (ASSOCIATED(force_env%meta_env)) THEN 1036 IF (force_env%meta_env%langevin) THEN 1037 DO ivar = 1, force_env%meta_env%n_colvar 1038 rand(ivar) = force_env%meta_env%rng(ivar)%next() 1039 ENDDO 1040 CALL metadyn_velocities_colvar(force_env, rand) 1041 ENDIF 1042 ENDIF 1043 1044 ! Velocity Verlet (first part) 1045 CALL vv_first(tmp, atomic_kind_set, local_particles, particle_set, & 1046 core_particle_set, shell_particle_set, nparticle_kind, shell_adiabatic, dt) 1047 1048 IF (simpar%variable_dt) CALL variable_timestep(md_env, tmp, dt, simpar, para_env, atomic_kind_set, & 1049 local_particles, particle_set, core_particle_set, shell_particle_set, & 1050 nparticle_kind, shell_adiabatic) 1051 1052 IF (simpar%constraint) THEN 1053 ! Possibly update the target values 1054 CALL shake_update_targets(gci, local_molecules, molecule_set, & 1055 molecule_kind_set, dt, force_env%root_section) 1056 1057 CALL shake_control(gci, local_molecules, molecule_set, & 1058 molecule_kind_set, particle_set, tmp%pos, tmp%vel, dt, simpar%shake_tol, & 1059 simpar%info_constraint, simpar%lagrange_multipliers, simpar%dump_lm, & 1060 cell, para_env%group, local_particles) 1061 END IF 1062 1063 ! Broadcast the new particle positions and deallocate pos components of temporary 1064 CALL update_dealloc_tmp(tmp, particle_set, shell_particle_set, & 1065 core_particle_set, para_env, shell_adiabatic, pos=.TRUE.) 1066 1067 IF (shell_adiabatic .AND. shell_check_distance) THEN 1068 CALL optimize_shell_core(force_env, particle_set, & 1069 shell_particle_set, core_particle_set, globenv, tmp=tmp, check=.TRUE.) 1070 END IF 1071 1072 ![ADAPT] update input structure with new coordinates, make new labels 1073 CALL qmmmx_update_force_env(force_env, force_env%root_section) 1074 1075 ![NB] recreate pointers changed by creation of new subsys in qmmm_update_force_mixing_env 1076 ![NB] ugly hack, which is why adaptivity isn't implemented in most other ensembles 1077 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1078 CALL force_env_get(force_env=force_env, subsys=subsys, cell=cell) 1079 1080 CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, & 1081 particles=particles, local_molecules=local_molecules, molecules=molecules, & 1082 molecule_kinds=molecule_kinds, gci=gci, virial=virial) 1083 1084 nparticle_kind = atomic_kinds%n_els 1085 atomic_kind_set => atomic_kinds%els 1086 molecule_kind_set => molecule_kinds%els 1087 1088 nparticle = particles%n_els 1089 particle_set => particles%els 1090 molecule_set => molecules%els 1091 1092 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, & 1093 shell_present=shell_present, shell_adiabatic=shell_adiabatic, & 1094 shell_check_distance=shell_check_distance) 1095 1096 ! Allocate work storage for positions and velocities 1097 IF (shell_present) THEN 1098 CALL cp_subsys_get(subsys=subsys, shell_particles=shell_particles, & 1099 core_particles=core_particles) 1100 shell_particle_set => shell_particles%els 1101 nshell = SIZE(shell_particles%els) 1102 1103 IF (shell_adiabatic) THEN 1104 core_particle_set => core_particles%els 1105 END IF 1106 END IF 1107 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1108 1109 ! Update forces 1110 ![NB] let nvt work with force mixing which does not have consistent energies and forces 1111 CALL force_env_calc_energy_force(force_env, require_consistent_energy_force=.FALSE.) 1112 1113 ! Metadynamics 1114 CALL metadyn_integrator(force_env, itimes, tmp%vel, rand=rand) 1115 1116 ! Velocity Verlet (second part) 1117 CALL vv_second(tmp, atomic_kind_set, local_particles, particle_set, & 1118 core_particle_set, shell_particle_set, nparticle_kind, shell_adiabatic, dt) 1119 1120 IF (simpar%constraint) CALL rattle_control(gci, local_molecules, molecule_set, & 1121 molecule_kind_set, particle_set, tmp%vel, dt, simpar%shake_tol, & 1122 simpar%info_constraint, simpar%lagrange_multipliers, simpar%dump_lm, & 1123 cell, para_env%group, local_particles) 1124 1125 ! Apply Thermostat over the full set of particles 1126 IF (shell_adiabatic) THEN 1127 CALL apply_thermostat_particles(thermostat_part, force_env, molecule_kind_set, molecule_set, & 1128 particle_set, local_molecules, local_particles, para_env%group, shell_adiabatic=shell_adiabatic, & 1129 vel=tmp%vel, shell_vel=tmp%shell_vel, core_vel=tmp%core_vel) 1130 1131 CALL apply_thermostat_shells(thermostat_shell, atomic_kind_set, particle_set, & 1132 local_particles, para_env%group, vel=tmp%vel, shell_vel=tmp%shell_vel, & 1133 core_vel=tmp%core_vel) 1134 ELSE 1135 CALL apply_thermostat_particles(thermostat_part, force_env, molecule_kind_set, molecule_set, & 1136 particle_set, local_molecules, local_particles, para_env%group, vel=tmp%vel) 1137 END IF 1138 1139 ! Broadcast the new particle velocities and deallocate temporary 1140 CALL update_dealloc_tmp(tmp, particle_set, shell_particle_set, & 1141 core_particle_set, para_env, shell_adiabatic, vel=.TRUE.) 1142 1143 IF (ASSOCIATED(force_env%meta_env)) THEN 1144 IF (force_env%meta_env%langevin) THEN 1145 DEALLOCATE (rand) 1146 ENDIF 1147 ENDIF 1148 1149 ! Update constraint virial 1150 IF (simpar%constraint) CALL pv_constraint(gci, local_molecules, & 1151 molecule_set, molecule_kind_set, particle_set, virial, para_env%group) 1152 1153 ! ** Evaluate Virial 1154 CALL virial_evaluate(atomic_kind_set, particle_set, & 1155 local_particles, virial, para_env%group) 1156 1157 END SUBROUTINE nvt 1158 1159! ************************************************************************************************** 1160!> \brief npt_i integrator for particle positions & momenta 1161!> isotropic box changes 1162!> \param md_env ... 1163!> \param globenv ... 1164!> \par History 1165!> none 1166!> \author CJM 1167! ************************************************************************************************** 1168 SUBROUTINE npt_i(md_env, globenv) 1169 1170 TYPE(md_environment_type), POINTER :: md_env 1171 TYPE(global_environment_type), POINTER :: globenv 1172 1173 CHARACTER(len=*), PARAMETER :: routineN = 'npt_i', routineP = moduleN//':'//routineN 1174 REAL(KIND=dp), PARAMETER :: e2 = 1.0_dp/6.0_dp, e4 = e2/20.0_dp, & 1175 e6 = e4/42.0_dp, e8 = e6/72.0_dp 1176 1177 INTEGER :: iroll, ivar, nparticle, nparticle_kind, & 1178 nshell 1179 INTEGER, POINTER :: itimes 1180 LOGICAL :: first, first_time, shell_adiabatic, & 1181 shell_check_distance, shell_present 1182 REAL(KIND=dp) :: dt, infree, kin, roll_tol, roll_tol_thrs 1183 REAL(KIND=dp), DIMENSION(3) :: vector_r, vector_v 1184 REAL(KIND=dp), DIMENSION(3, 3) :: pv_kin 1185 REAL(KIND=dp), DIMENSION(:), POINTER :: rand 1186 REAL(KIND=dp), SAVE :: eps_0 1187 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds 1188 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 1189 TYPE(cell_type), POINTER :: cell 1190 TYPE(cp_para_env_type), POINTER :: para_env 1191 TYPE(cp_subsys_type), POINTER :: subsys 1192 TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles 1193 TYPE(force_env_type), POINTER :: force_env 1194 TYPE(global_constraint_type), POINTER :: gci 1195 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds 1196 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set 1197 TYPE(molecule_list_type), POINTER :: molecules 1198 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 1199 TYPE(npt_info_type), POINTER :: npt(:, :) 1200 TYPE(old_variables_type), POINTER :: old 1201 TYPE(particle_list_type), POINTER :: core_particles, particles, & 1202 shell_particles 1203 TYPE(particle_type), DIMENSION(:), POINTER :: core_particle_set, particle_set, & 1204 shell_particle_set 1205 TYPE(simpar_type), POINTER :: simpar 1206 TYPE(thermostat_type), POINTER :: thermostat_baro, thermostat_part, & 1207 thermostat_shell 1208 TYPE(tmp_variables_type), POINTER :: tmp 1209 TYPE(virial_type), POINTER :: virial 1210 1211 NULLIFY (gci, thermostat_baro, thermostat_part, thermostat_shell, force_env) 1212 NULLIFY (atomic_kinds, cell, para_env, subsys, local_molecules, local_particles) 1213 NULLIFY (molecule_kinds, molecules, molecule_kind_set, npt) 1214 NULLIFY (core_particles, particles, shell_particles, tmp, old) 1215 NULLIFY (core_particle_set, particle_set, shell_particle_set) 1216 NULLIFY (simpar, virial, rand, itimes) 1217 1218 CALL get_md_env(md_env=md_env, simpar=simpar, force_env=force_env, & 1219 thermostat_part=thermostat_part, thermostat_baro=thermostat_baro, & 1220 thermostat_shell=thermostat_shell, npt=npt, first_time=first_time, & 1221 para_env=para_env, itimes=itimes) 1222 dt = simpar%dt 1223 infree = 1.0_dp/REAL(simpar%nfree, KIND=dp) 1224 1225 CALL force_env_get(force_env=force_env, subsys=subsys, cell=cell) 1226 1227 ! Do some checks on coordinates and box 1228 CALL apply_qmmm_walls_reflective(force_env) 1229 1230 CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, & 1231 particles=particles, local_molecules=local_molecules, molecules=molecules, & 1232 gci=gci, molecule_kinds=molecule_kinds, virial=virial) 1233 1234 nparticle_kind = atomic_kinds%n_els 1235 atomic_kind_set => atomic_kinds%els 1236 molecule_kind_set => molecule_kinds%els 1237 1238 nparticle = particles%n_els 1239 particle_set => particles%els 1240 molecule_set => molecules%els 1241 1242 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, & 1243 shell_present=shell_present, shell_adiabatic=shell_adiabatic, & 1244 shell_check_distance=shell_check_distance) 1245 1246 IF (first_time) THEN 1247 CALL virial_evaluate(atomic_kind_set, particle_set, & 1248 local_particles, virial, para_env%group) 1249 END IF 1250 1251 ! Allocate work storage for positions and velocities 1252 CALL allocate_old(old, particle_set, npt) 1253 1254 IF (ASSOCIATED(force_env%meta_env)) THEN 1255 ! Allocate random number for Langevin Thermostat acting on COLVARS 1256 IF (force_env%meta_env%langevin) THEN 1257 ALLOCATE (rand(force_env%meta_env%n_colvar)) 1258 rand(:) = 0.0_dp 1259 ENDIF 1260 ENDIF 1261 1262 IF (shell_present) THEN 1263 CALL cp_subsys_get(subsys=subsys, & 1264 shell_particles=shell_particles, core_particles=core_particles) 1265 shell_particle_set => shell_particles%els 1266 nshell = SIZE(shell_particles%els) 1267 IF (shell_adiabatic) THEN 1268 core_particle_set => core_particles%els 1269 END IF 1270 END IF 1271 1272 CALL allocate_tmp(md_env, tmp, nparticle, nshell, shell_adiabatic) 1273 1274 ! Initialize eps_0 the first time through 1275 IF (first_time) eps_0 = npt(1, 1)%eps 1276 1277 ! Apply thermostat to barostat 1278 CALL apply_thermostat_baro(thermostat_baro, npt, para_env%group) 1279 1280 ! Apply Thermostat over the full set of particles 1281 IF (simpar%ensemble /= npe_i_ensemble) THEN 1282 IF (shell_adiabatic) THEN 1283 CALL apply_thermostat_particles(thermostat_part, force_env, molecule_kind_set, molecule_set, & 1284 particle_set, local_molecules, local_particles, para_env%group, shell_adiabatic=shell_adiabatic, & 1285 shell_particle_set=shell_particle_set, core_particle_set=core_particle_set) 1286 1287 ELSE 1288 CALL apply_thermostat_particles(thermostat_part, force_env, molecule_kind_set, molecule_set, & 1289 particle_set, local_molecules, local_particles, para_env%group) 1290 END IF 1291 END IF 1292 1293 ! Apply Thermostat over the core-shell motion 1294 CALL apply_thermostat_shells(thermostat_shell, atomic_kind_set, particle_set, & 1295 local_particles, para_env%group, shell_particle_set=shell_particle_set, & 1296 core_particle_set=core_particle_set) 1297 1298 IF (simpar%constraint) THEN 1299 ! Possibly update the target values 1300 CALL shake_update_targets(gci, local_molecules, molecule_set, & 1301 molecule_kind_set, dt, force_env%root_section) 1302 END IF 1303 1304 ! setting up for ROLL: saving old variables 1305 IF (simpar%constraint) THEN 1306 roll_tol_thrs = simpar%roll_tol 1307 iroll = 1 1308 CALL set(old, atomic_kind_set, particle_set, local_particles, cell, npt, 'F') 1309 CALL getold(gci, local_molecules, molecule_set, & 1310 molecule_kind_set, particle_set, cell) 1311 ELSE 1312 roll_tol_thrs = EPSILON(0.0_dp) 1313 ENDIF 1314 roll_tol = -roll_tol_thrs 1315 1316 ! *** Velocity Verlet for Langeving *** v(t)--> v(t+1/2) 1317 IF (ASSOCIATED(force_env%meta_env)) THEN 1318 IF (force_env%meta_env%langevin) THEN 1319 DO ivar = 1, force_env%meta_env%n_colvar 1320 rand(ivar) = force_env%meta_env%rng(ivar)%next() 1321 ENDDO 1322 CALL metadyn_velocities_colvar(force_env, rand) 1323 ENDIF 1324 ENDIF 1325 1326 SR: DO WHILE (ABS(roll_tol) >= roll_tol_thrs) ! SHAKE-ROLL LOOP 1327 1328 IF (simpar%constraint) THEN 1329 CALL set(old, atomic_kind_set, particle_set, local_particles, cell, npt, 'B') 1330 END IF 1331 1332 CALL update_pv(gci, simpar, atomic_kind_set, particle_set, & 1333 local_molecules, molecule_set, molecule_kind_set, & 1334 local_particles, kin, pv_kin, virial, para_env%group) 1335 CALL update_veps(cell, npt, simpar, pv_kin, kin, virial, infree) 1336 1337 tmp%arg_r(1) = (0.5_dp*npt(1, 1)%v*dt)* & 1338 (0.5_dp*npt(1, 1)%v*dt) 1339 tmp%poly_r(1:3) = 1.0_dp + e2*tmp%arg_r(1) + e4*tmp%arg_r(1)*tmp%arg_r(1) + & 1340 e6*tmp%arg_r(1)**3 + e8*tmp%arg_r(1)**4 1341 1342 tmp%arg_v(1) = (0.25_dp*npt(1, 1)%v*dt* & 1343 (1.0_dp + 3.0_dp*infree))*(0.25_dp*npt(1, 1)%v* & 1344 dt*(1.0_dp + 3.0_dp*infree)) 1345 tmp%poly_v(1:3) = 1.0_dp + e2*tmp%arg_v(1) + e4*tmp%arg_v(1)*tmp%arg_v(1) + & 1346 e6*tmp%arg_v(1)**3 + e8*tmp%arg_v(1)**4 1347 1348 tmp%scale_r(1:3) = EXP(0.5_dp*dt*npt(1, 1)%v) 1349 tmp%scale_v(1:3) = EXP(-0.25_dp*dt*npt(1, 1)%v* & 1350 (1.0_dp + 3.0_dp*infree)) 1351 1352 ! first half of velocity verlet 1353 CALL vv_first(tmp, atomic_kind_set, local_particles, particle_set, & 1354 core_particle_set, shell_particle_set, nparticle_kind, & 1355 shell_adiabatic, dt) 1356 1357 IF (simpar%variable_dt) CALL variable_timestep(md_env, tmp, dt, simpar, para_env, & 1358 atomic_kind_set, local_particles, particle_set, core_particle_set, & 1359 shell_particle_set, nparticle_kind, shell_adiabatic, npt=npt) 1360 1361 roll_tol = 0.0_dp 1362 vector_r(:) = tmp%scale_r(:)*tmp%poly_r(:) 1363 vector_v(:) = tmp%scale_v(:)*tmp%poly_v(:) 1364 1365 IF (simpar%constraint) CALL shake_roll_control(gci, local_molecules, & 1366 molecule_set, molecule_kind_set, particle_set, tmp%pos, tmp%vel, dt, simpar, & 1367 roll_tol, iroll, vector_r, vector_v, para_env%group, cell=cell, & 1368 local_particles=local_particles) 1369 END DO SR 1370 1371 ! Update eps: 1372 npt(:, :)%eps = npt(:, :)%eps + dt*npt(:, :)%v 1373 1374 ! Update h_mat 1375 cell%hmat(:, :) = cell%hmat(:, :)*EXP(npt(1, 1)%eps - eps_0) 1376 1377 eps_0 = npt(1, 1)%eps 1378 1379 ! Update the inverse 1380 CALL init_cell(cell) 1381 1382 ! Broadcast the new particle positions and deallocate the pos components of temporary 1383 CALL update_dealloc_tmp(tmp, particle_set, shell_particle_set, & 1384 core_particle_set, para_env, shell_adiabatic, pos=.TRUE.) 1385 1386 IF (shell_adiabatic .AND. shell_check_distance) THEN 1387 CALL optimize_shell_core(force_env, particle_set, & 1388 shell_particle_set, core_particle_set, globenv, tmp=tmp, check=.TRUE.) 1389 END IF 1390 1391 ! Update forces 1392 CALL force_env_calc_energy_force(force_env) 1393 1394 ! Metadynamics 1395 CALL metadyn_integrator(force_env, itimes, tmp%vel, rand=rand) 1396 1397 ! Velocity Verlet (second part) 1398 CALL vv_second(tmp, atomic_kind_set, local_particles, particle_set, & 1399 core_particle_set, shell_particle_set, nparticle_kind, & 1400 shell_adiabatic, dt) 1401 1402 IF (simpar%constraint) THEN 1403 roll_tol_thrs = simpar%roll_tol 1404 first = .TRUE. 1405 iroll = 1 1406 CALL set(old, atomic_kind_set, particle_set, tmp%vel, local_particles, cell, npt, 'F') 1407 ELSE 1408 roll_tol_thrs = EPSILON(0.0_dp) 1409 ENDIF 1410 roll_tol = -roll_tol_thrs 1411 1412 RR: DO WHILE (ABS(roll_tol) >= roll_tol_thrs) ! RATTLE-ROLL LOOP 1413 roll_tol = 0.0_dp 1414 IF (simpar%constraint) CALL rattle_roll_setup(old, gci, atomic_kind_set, & 1415 particle_set, local_particles, molecule_kind_set, molecule_set, & 1416 local_molecules, tmp%vel, dt, cell, npt, simpar, virial, vector_v, & 1417 roll_tol, iroll, infree, first, para_env) 1418 1419 CALL update_pv(gci, simpar, atomic_kind_set, tmp%vel, particle_set, & 1420 local_molecules, molecule_set, molecule_kind_set, & 1421 local_particles, kin, pv_kin, virial, para_env%group) 1422 CALL update_veps(cell, npt, simpar, pv_kin, kin, virial, infree) 1423 END DO RR 1424 1425 ! Apply Thermostat over the full set of particles 1426 IF (simpar%ensemble /= npe_i_ensemble) THEN 1427 IF (shell_adiabatic) THEN 1428 CALL apply_thermostat_particles(thermostat_part, force_env, molecule_kind_set, molecule_set, & 1429 particle_set, local_molecules, local_particles, para_env%group, shell_adiabatic=shell_adiabatic, & 1430 vel=tmp%vel, shell_vel=tmp%shell_vel, core_vel=tmp%core_vel) 1431 ELSE 1432 CALL apply_thermostat_particles(thermostat_part, force_env, molecule_kind_set, molecule_set, & 1433 particle_set, local_molecules, local_particles, para_env%group, vel=tmp%vel) 1434 END IF 1435 END IF 1436 1437 ! Apply Thermostat over the core-shell motion 1438 IF (ASSOCIATED(thermostat_shell)) THEN 1439 CALL apply_thermostat_shells(thermostat_shell, atomic_kind_set, particle_set, & 1440 local_particles, para_env%group, vel=tmp%vel, shell_vel=tmp%shell_vel, & 1441 core_vel=tmp%core_vel) 1442 END IF 1443 1444 ! Apply Thermostat to Barostat 1445 CALL apply_thermostat_baro(thermostat_baro, npt, para_env%group) 1446 1447 ! Annealing of particle velocities is only possible when no thermostat is active 1448 IF (simpar%ensemble == npe_i_ensemble .AND. simpar%annealing) THEN 1449 tmp%vel(:, :) = tmp%vel(:, :)*simpar%f_annealing 1450 IF (shell_adiabatic) THEN 1451 CALL shell_scale_comv(atomic_kind_set, local_particles, particle_set, & 1452 tmp%vel, tmp%shell_vel, tmp%core_vel) 1453 END IF 1454 END IF 1455 ! Annealing of CELL velocities is only possible when no thermostat is active 1456 IF (simpar%ensemble == npe_i_ensemble .AND. simpar%annealing_cell) THEN 1457 npt(1, 1)%v = npt(1, 1)%v*simpar%f_annealing_cell 1458 END IF 1459 1460 ! Broadcast the new particle velocities and deallocate temporary 1461 CALL update_dealloc_tmp(tmp, particle_set, shell_particle_set, & 1462 core_particle_set, para_env, shell_adiabatic, vel=.TRUE.) 1463 1464 ! Update constraint virial 1465 IF (simpar%constraint) CALL pv_constraint(gci, local_molecules, & 1466 molecule_set, molecule_kind_set, particle_set, virial, para_env%group) 1467 1468 CALL virial_evaluate(atomic_kind_set, particle_set, & 1469 local_particles, virial, para_env%group) 1470 1471 ! Deallocate old variables 1472 CALL deallocate_old(old) 1473 1474 IF (ASSOCIATED(force_env%meta_env)) THEN 1475 IF (force_env%meta_env%langevin) THEN 1476 DEALLOCATE (rand) 1477 ENDIF 1478 ENDIF 1479 1480 IF (first_time) THEN 1481 first_time = .FALSE. 1482 CALL set_md_env(md_env, first_time=first_time) 1483 END IF 1484 1485 END SUBROUTINE npt_i 1486 1487! ************************************************************************************************** 1488!> \brief uses coordinates in a file and generates frame after frame of these 1489!> \param md_env ... 1490!> \par History 1491!> - 04.2005 created [Joost VandeVondele] 1492!> - modified to make it more general [MI] 1493!> \note 1494!> it can be used to compute some properties on already available trajectories 1495! ************************************************************************************************** 1496 SUBROUTINE reftraj(md_env) 1497 TYPE(md_environment_type), POINTER :: md_env 1498 1499 CHARACTER(len=*), PARAMETER :: routineN = 'reftraj', routineP = moduleN//':'//routineN 1500 1501 CHARACTER(LEN=10) :: AA 1502 INTEGER :: cell_itimes, I, nparticle, Nread, & 1503 trj_itimes 1504 INTEGER, POINTER :: itimes 1505 LOGICAL :: init, my_end, test_ok 1506 REAL(KIND=dp) :: cell_time, h(3, 3), trj_epot, trj_time, & 1507 vol 1508 REAL(KIND=dp), POINTER :: time 1509 TYPE(cell_type), POINTER :: cell 1510 TYPE(cp_para_env_type), POINTER :: para_env 1511 TYPE(cp_subsys_type), POINTER :: subsys 1512 TYPE(force_env_type), POINTER :: force_env 1513 TYPE(particle_list_type), POINTER :: particles 1514 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 1515 TYPE(reftraj_type), POINTER :: reftraj_env 1516 TYPE(simpar_type), POINTER :: simpar 1517 1518 NULLIFY (reftraj_env, particle_set, particles, force_env, subsys, simpar, para_env, cell, itimes, time) 1519 CALL get_md_env(md_env=md_env, init=init, reftraj=reftraj_env, force_env=force_env, & 1520 para_env=para_env, simpar=simpar) 1521 1522 CALL force_env_get(force_env=force_env, cell=cell, subsys=subsys) 1523 reftraj_env%isnap = reftraj_env%isnap + reftraj_env%info%stride 1524 1525 ! Do some checks on coordinates and box 1526 CALL apply_qmmm_walls_reflective(force_env) 1527 CALL cp_subsys_get(subsys=subsys, particles=particles) 1528 nparticle = particles%n_els 1529 particle_set => particles%els 1530 1531 ! SnapShots read from an external file (parsers calls are buffered! please 1532 ! don't put any additional MPI call!) [tlaino] 1533 CALL parser_read_line(reftraj_env%info%traj_parser, 1) 1534 READ (reftraj_env%info%traj_parser%input_line, FMT="(I8)") nread 1535 CPASSERT(nread == nparticle) 1536 CALL parser_read_line(reftraj_env%info%traj_parser, 1) 1537 test_ok = .FALSE. 1538 READ (reftraj_env%info%traj_parser%input_line, FMT="(T6,I8,T23,F12.3,T41,F20.10)", ERR=999) & 1539 trj_itimes, trj_time, trj_epot 1540 test_ok = .TRUE. 1541999 IF (.NOT. test_ok) THEN 1542 ! Handling properly the error when reading the title of an XYZ 1543 CALL get_md_env(md_env, itimes=itimes) 1544 trj_itimes = itimes 1545 trj_time = 0.0_dp 1546 trj_epot = 0.0_dp 1547 END IF 1548 DO i = 1, nread - 1 1549 CALL parser_read_line(reftraj_env%info%traj_parser, 1) 1550 READ (reftraj_env%info%traj_parser%input_line(1:LEN_TRIM(reftraj_env%info%traj_parser%input_line)), *) AA, particle_set(i)%r 1551 particle_set(i)%r = particle_set(i)%r*bohr 1552 END DO 1553 ! End of file is properly addressed in the previous call.. 1554 ! Let's check directly (providing some info) also for the last 1555 ! line of this frame.. 1556 CALL parser_read_line(reftraj_env%info%traj_parser, 1, at_end=my_end) 1557 READ (unit=reftraj_env%info%traj_parser%input_line, fmt=*) AA, particle_set(i)%r 1558 particle_set(i)%r(1) = cp_unit_to_cp2k(particle_set(i)%r(1), "angstrom") 1559 particle_set(i)%r(2) = cp_unit_to_cp2k(particle_set(i)%r(2), "angstrom") 1560 particle_set(i)%r(3) = cp_unit_to_cp2k(particle_set(i)%r(3), "angstrom") 1561 1562 ! Check if we reached the end of the file and provide some info.. 1563 IF (my_end) THEN 1564 IF (reftraj_env%isnap /= (simpar%nsteps - 1)) & 1565 CALL cp_abort(__LOCATION__, & 1566 "Reached the end of the Trajectory frames in the TRAJECTORY file. Number of "// & 1567 "missing frames ("//cp_to_string((simpar%nsteps - 1) - reftraj_env%isnap)//").") 1568 END IF 1569 1570 IF (reftraj_env%info%variable_volume) THEN 1571 CALL parser_get_next_line(reftraj_env%info%cell_parser, 1, at_end=my_end) 1572 CALL parse_cell_line(reftraj_env%info%cell_parser%input_line, cell_itimes, cell_time, h, vol) 1573 CPASSERT(trj_itimes == cell_itimes) 1574 ! Check if we reached the end of the file and provide some info.. 1575 IF (my_end) THEN 1576 IF (reftraj_env%isnap /= (simpar%nsteps - 1)) & 1577 CALL cp_abort(__LOCATION__, & 1578 "Reached the end of the cell info frames in the CELL file. Number of "// & 1579 "missing frames ("//cp_to_string((simpar%nsteps - 1) - reftraj_env%isnap)//").") 1580 END IF 1581 END IF 1582 1583 IF (init) THEN 1584 reftraj_env%time0 = trj_time 1585 reftraj_env%epot0 = trj_epot 1586 reftraj_env%itimes0 = trj_itimes 1587 END IF 1588 1589 IF (trj_itimes /= 0.0_dp .AND. trj_time /= 0.0_dp) simpar%dt = (trj_time/femtoseconds)/REAL(trj_itimes, KIND=dp) 1590 1591 reftraj_env%epot = trj_epot 1592 reftraj_env%itimes = trj_itimes 1593 reftraj_env%time = trj_time/femtoseconds 1594 CALL get_md_env(md_env, t=time) 1595 time = reftraj_env%time 1596 1597 IF (reftraj_env%info%variable_volume) THEN 1598 cell%hmat = h 1599 CALL init_cell(cell) 1600 END IF 1601 1602 ![ADAPT] update input structure with new coordinates, make new labels 1603 CALL qmmmx_update_force_env(force_env, force_env%root_section) 1604 ! no pointers into force_env%subsys to update 1605 1606 ! Task to perform on the reference trajectory 1607 ! Compute energy and forces 1608 ![NB] let reftraj work with force mixing which does not have consistent energies and forces 1609 CALL force_env_calc_energy_force(force_env, calc_force=reftraj_env%info%eval_forces, eval_energy_forces=reftraj_env%info%eval_EF, & 1610 require_consistent_energy_force=.FALSE.) 1611 1612 ! Metadynamics 1613 CALL metadyn_integrator(force_env, trj_itimes) 1614 1615 ! Compute MSD with respect to a reference configuration 1616 IF (reftraj_env%info%msd) THEN 1617 CALL compute_msd_reftraj(reftraj_env, md_env, particle_set) 1618 END IF 1619 1620 ! Skip according the stride both Trajectory and Cell (if possible) 1621 CALL parser_get_next_line(reftraj_env%info%traj_parser, (reftraj_env%info%stride - 1)*(nparticle + 2)) 1622 IF (reftraj_env%info%variable_volume) THEN 1623 CALL parser_get_next_line(reftraj_env%info%cell_parser, (reftraj_env%info%stride - 1)) 1624 END IF 1625 END SUBROUTINE reftraj 1626 1627! ************************************************************************************************** 1628!> \brief nph_uniaxial integrator (non-Hamiltonian version) 1629!> for particle positions & momenta undergoing 1630!> uniaxial stress ( in x-direction of orthorhombic cell) 1631!> due to a shock compression: 1632!> Reed et. al. Physical Review Letters 90, 235503 (2003). 1633!> \param md_env ... 1634!> \par History 1635!> none 1636!> \author CJM 1637! ************************************************************************************************** 1638 SUBROUTINE nph_uniaxial(md_env) 1639 1640 TYPE(md_environment_type), POINTER :: md_env 1641 1642 CHARACTER(len=*), PARAMETER :: routineN = 'nph_uniaxial', routineP = moduleN//':'//routineN 1643 REAL(dp), PARAMETER :: e2 = 1._dp/6._dp, e4 = e2/20._dp, & 1644 e6 = e4/42._dp, e8 = e6/72._dp 1645 1646 INTEGER :: iroll, nparticle, nparticle_kind, nshell 1647 INTEGER, POINTER :: itimes 1648 LOGICAL :: first, first_time, shell_adiabatic, & 1649 shell_present 1650 REAL(KIND=dp) :: dt, infree, kin, roll_tol, roll_tol_thrs 1651 REAL(KIND=dp), DIMENSION(3) :: vector_r, vector_v 1652 REAL(KIND=dp), DIMENSION(3, 3) :: pv_kin 1653 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds 1654 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 1655 TYPE(cell_type), POINTER :: cell 1656 TYPE(cp_para_env_type), POINTER :: para_env 1657 TYPE(cp_subsys_type), POINTER :: subsys 1658 TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles 1659 TYPE(force_env_type), POINTER :: force_env 1660 TYPE(global_constraint_type), POINTER :: gci 1661 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds 1662 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set 1663 TYPE(molecule_list_type), POINTER :: molecules 1664 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 1665 TYPE(npt_info_type), POINTER :: npt(:, :) 1666 TYPE(old_variables_type), POINTER :: old 1667 TYPE(particle_list_type), POINTER :: core_particles, particles, & 1668 shell_particles 1669 TYPE(particle_type), DIMENSION(:), POINTER :: core_particle_set, particle_set, & 1670 shell_particle_set 1671 TYPE(simpar_type), POINTER :: simpar 1672 TYPE(tmp_variables_type), POINTER :: tmp 1673 TYPE(virial_type), POINTER :: virial 1674 1675 NULLIFY (gci, force_env) 1676 NULLIFY (atomic_kinds, cell, para_env, subsys, local_molecules, local_particles) 1677 NULLIFY (molecule_kinds, molecules, molecule_kind_set, npt) 1678 NULLIFY (core_particles, particles, shell_particles, tmp, old) 1679 NULLIFY (core_particle_set, particle_set, shell_particle_set) 1680 NULLIFY (simpar, virial, itimes) 1681 1682 CALL get_md_env(md_env=md_env, simpar=simpar, force_env=force_env, npt=npt, & 1683 first_time=first_time, para_env=para_env, itimes=itimes) 1684 dt = simpar%dt 1685 infree = 1.0_dp/REAL(simpar%nfree, dp) 1686 1687 CALL force_env_get(force_env, subsys=subsys, cell=cell) 1688 1689 ! Do some checks on coordinates and box 1690 CALL apply_qmmm_walls_reflective(force_env) 1691 1692 CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, & 1693 particles=particles, local_molecules=local_molecules, molecules=molecules, gci=gci, & 1694 molecule_kinds=molecule_kinds, virial=virial) 1695 1696 nparticle_kind = atomic_kinds%n_els 1697 atomic_kind_set => atomic_kinds%els 1698 molecule_kind_set => molecule_kinds%els 1699 1700 nparticle = particles%n_els 1701 particle_set => particles%els 1702 molecule_set => molecules%els 1703 1704 IF (first_time) THEN 1705 CALL virial_evaluate(atomic_kind_set, particle_set, & 1706 local_particles, virial, para_env%group) 1707 END IF 1708 1709 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, & 1710 shell_present=shell_present, shell_adiabatic=shell_adiabatic) 1711 1712 ! Allocate work storage for positions and velocities 1713 CALL allocate_old(old, particle_set, npt) 1714 1715 IF (shell_present) THEN 1716 CALL cp_subsys_get(subsys=subsys, & 1717 shell_particles=shell_particles, core_particles=core_particles) 1718 shell_particle_set => shell_particles%els 1719 nshell = SIZE(shell_particles%els) 1720 IF (shell_adiabatic) THEN 1721 core_particle_set => core_particles%els 1722 END IF 1723 END IF 1724 1725 CALL allocate_tmp(md_env, tmp, nparticle, nshell, shell_adiabatic) 1726 1727 IF (simpar%constraint) THEN 1728 ! Possibly update the target values 1729 CALL shake_update_targets(gci, local_molecules, molecule_set, & 1730 molecule_kind_set, dt, force_env%root_section) 1731 END IF 1732 1733 ! setting up for ROLL: saving old variables 1734 IF (simpar%constraint) THEN 1735 roll_tol_thrs = simpar%roll_tol 1736 iroll = 1 1737 CALL set(old, atomic_kind_set, particle_set, local_particles, cell, npt, 'F') 1738 CALL getold(gci, local_molecules, molecule_set, & 1739 molecule_kind_set, particle_set, cell) 1740 ELSE 1741 roll_tol_thrs = EPSILON(0.0_dp) 1742 ENDIF 1743 roll_tol = -roll_tol_thrs 1744 1745 SR: DO WHILE (ABS(roll_tol) >= roll_tol_thrs) ! SHAKE-ROLL LOOP 1746 1747 IF (simpar%constraint) THEN 1748 CALL set(old, atomic_kind_set, particle_set, local_particles, cell, npt, 'B') 1749 END IF 1750 CALL update_pv(gci, simpar, atomic_kind_set, particle_set, & 1751 local_molecules, molecule_set, molecule_kind_set, & 1752 local_particles, kin, pv_kin, virial, para_env%group) 1753 CALL update_veps(cell, npt, simpar, pv_kin, kin, virial, infree) 1754 1755 tmp%arg_r(1) = (0.5_dp*npt(1, 1)%v*dt)* & 1756 (0.5_dp*npt(1, 1)%v*dt) 1757 tmp%poly_r(1) = 1._dp + e2*tmp%arg_r(1) + e4*tmp%arg_r(1)*tmp%arg_r(1) + & 1758 e6*tmp%arg_r(1)**3 + e8*tmp%arg_r(1)**4 1759 tmp%poly_r(2) = 1.0_dp 1760 tmp%poly_r(3) = 1.0_dp 1761 1762 tmp%arg_v(1) = (0.25_dp*npt(1, 1)%v*dt* & 1763 (1._dp + infree))*(0.25_dp*npt(1, 1)%v* & 1764 dt*(1._dp + infree)) 1765 tmp%arg_v(2) = (0.25_dp*npt(1, 1)%v*dt*infree)* & 1766 (0.25_dp*npt(1, 1)%v*dt*infree) 1767 tmp%poly_v(1) = 1._dp + e2*tmp%arg_v(1) + e4*tmp%arg_v(1)*tmp%arg_v(1) + & 1768 e6*tmp%arg_v(1)**3 + e8*tmp%arg_v(1)**4 1769 tmp%poly_v(2) = 1._dp + e2*tmp%arg_v(2) + e4*tmp%arg_v(2)*tmp%arg_v(2) + & 1770 e6*tmp%arg_v(2)**3 + e8*tmp%arg_v(2)**4 1771 tmp%poly_v(3) = 1._dp + e2*tmp%arg_v(2) + e4*tmp%arg_v(2)*tmp%arg_v(2) + & 1772 e6*tmp%arg_v(2)**3 + e8*tmp%arg_v(2)**4 1773 1774 tmp%scale_r(1) = EXP(0.5_dp*dt*npt(1, 1)%v) 1775 tmp%scale_r(2) = 1.0_dp 1776 tmp%scale_r(3) = 1.0_dp 1777 1778 tmp%scale_v(1) = EXP(-0.25_dp*dt*npt(1, 1)%v* & 1779 (1._dp + infree)) 1780 tmp%scale_v(2) = EXP(-0.25_dp*dt*npt(1, 1)%v*infree) 1781 tmp%scale_v(3) = EXP(-0.25_dp*dt*npt(1, 1)%v*infree) 1782 1783 ! first half of velocity verlet 1784 CALL vv_first(tmp, atomic_kind_set, local_particles, particle_set, & 1785 core_particle_set, shell_particle_set, nparticle_kind, & 1786 shell_adiabatic, dt) 1787 1788 IF (simpar%variable_dt) CALL variable_timestep(md_env, tmp, dt, simpar, para_env, & 1789 atomic_kind_set, local_particles, particle_set, core_particle_set, & 1790 shell_particle_set, nparticle_kind, shell_adiabatic, npt=npt) 1791 1792 roll_tol = 0._dp 1793 vector_r(:) = 0._dp 1794 vector_v(:) = tmp%scale_v(:)*tmp%poly_v(:) 1795 vector_r(1) = tmp%scale_r(1)*tmp%poly_r(1) 1796 1797 IF (simpar%constraint) CALL shake_roll_control(gci, local_molecules, & 1798 molecule_set, molecule_kind_set, particle_set, tmp%pos, tmp%vel, dt, simpar, & 1799 roll_tol, iroll, vector_r, vector_v, para_env%group, cell=cell, & 1800 local_particles=local_particles) 1801 END DO SR 1802 1803 ! Update h_mat 1804 cell%hmat(1, 1) = cell%hmat(1, 1)*tmp%scale_r(1)*tmp%scale_r(1) 1805 1806 ! Update the cell 1807 CALL init_cell(cell) 1808 1809 ! Broadcast the new particle positions and deallocate the pos component of temporary 1810 CALL update_dealloc_tmp(tmp, particle_set, shell_particle_set, & 1811 core_particle_set, para_env, shell_adiabatic, pos=.TRUE.) 1812 1813 ! Update forces (and stress) 1814 CALL force_env_calc_energy_force(force_env) 1815 1816 ! Metadynamics 1817 CALL metadyn_integrator(force_env, itimes, tmp%vel) 1818 1819 ! Velocity Verlet (second part) 1820 CALL vv_second(tmp, atomic_kind_set, local_particles, particle_set, & 1821 core_particle_set, shell_particle_set, nparticle_kind, & 1822 shell_adiabatic, dt) 1823 1824 IF (simpar%constraint) THEN 1825 roll_tol_thrs = simpar%roll_tol 1826 first = .TRUE. 1827 iroll = 1 1828 CALL set(old, atomic_kind_set, particle_set, tmp%vel, local_particles, cell, npt, 'F') 1829 ELSE 1830 roll_tol_thrs = EPSILON(0.0_dp) 1831 ENDIF 1832 roll_tol = -roll_tol_thrs 1833 1834 RR: DO WHILE (ABS(roll_tol) >= roll_tol_thrs) ! RATTLE-ROLL LOOP 1835 roll_tol = 0._dp 1836 IF (simpar%constraint) CALL rattle_roll_setup(old, gci, atomic_kind_set, & 1837 particle_set, local_particles, molecule_kind_set, molecule_set, & 1838 local_molecules, tmp%vel, dt, cell, npt, simpar, virial, vector_v, & 1839 roll_tol, iroll, infree, first, para_env) 1840 1841 CALL update_pv(gci, simpar, atomic_kind_set, tmp%vel, particle_set, & 1842 local_molecules, molecule_set, molecule_kind_set, & 1843 local_particles, kin, pv_kin, virial, para_env%group) 1844 CALL update_veps(cell, npt, simpar, pv_kin, kin, virial, infree) 1845 END DO RR 1846 1847 IF (simpar%annealing) tmp%vel(:, :) = tmp%vel(:, :)*simpar%f_annealing 1848 1849 ! Broadcast the new particle velocities and deallocate the temporary 1850 CALL update_dealloc_tmp(tmp, particle_set, shell_particle_set, & 1851 core_particle_set, para_env, shell_adiabatic, vel=.TRUE.) 1852 1853 ! Update constraint virial 1854 IF (simpar%constraint) CALL pv_constraint(gci, local_molecules, & 1855 molecule_set, molecule_kind_set, particle_set, virial, para_env%group) 1856 1857 CALL virial_evaluate(atomic_kind_set, particle_set, & 1858 local_particles, virial, para_env%group) 1859 1860 ! Deallocate old variables 1861 CALL deallocate_old(old) 1862 1863 IF (first_time) THEN 1864 first_time = .FALSE. 1865 CALL set_md_env(md_env, first_time=first_time) 1866 END IF 1867 1868 END SUBROUTINE nph_uniaxial 1869 1870! ************************************************************************************************** 1871!> \brief nph_uniaxial integrator (non-Hamiltonian version) 1872!> for particle positions & momenta undergoing 1873!> uniaxial stress ( in x-direction of orthorhombic cell) 1874!> due to a shock compression: 1875!> Reed et. al. Physical Review Letters 90, 235503 (2003). 1876!> Added damping (e.g. thermostat to barostat) 1877!> \param md_env ... 1878!> \par History 1879!> none 1880!> \author CJM 1881! ************************************************************************************************** 1882 SUBROUTINE nph_uniaxial_damped(md_env) 1883 1884 TYPE(md_environment_type), POINTER :: md_env 1885 1886 CHARACTER(len=*), PARAMETER :: routineN = 'nph_uniaxial_damped', & 1887 routineP = moduleN//':'//routineN 1888 REAL(dp), PARAMETER :: e2 = 1._dp/6._dp, e4 = e2/20._dp, & 1889 e6 = e4/42._dp, e8 = e6/72._dp 1890 1891 INTEGER :: iroll, nparticle, nparticle_kind, nshell 1892 INTEGER, POINTER :: itimes 1893 LOGICAL :: first, first_time, shell_adiabatic, & 1894 shell_present 1895 REAL(KIND=dp) :: aa, aax, dt, gamma1, infree, kin, & 1896 roll_tol, roll_tol_thrs 1897 REAL(KIND=dp), DIMENSION(3) :: vector_r, vector_v 1898 REAL(KIND=dp), DIMENSION(3, 3) :: pv_kin 1899 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds 1900 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 1901 TYPE(cell_type), POINTER :: cell 1902 TYPE(cp_para_env_type), POINTER :: para_env 1903 TYPE(cp_subsys_type), POINTER :: subsys 1904 TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles 1905 TYPE(force_env_type), POINTER :: force_env 1906 TYPE(global_constraint_type), POINTER :: gci 1907 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds 1908 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set 1909 TYPE(molecule_list_type), POINTER :: molecules 1910 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 1911 TYPE(npt_info_type), POINTER :: npt(:, :) 1912 TYPE(old_variables_type), POINTER :: old 1913 TYPE(particle_list_type), POINTER :: core_particles, particles, & 1914 shell_particles 1915 TYPE(particle_type), DIMENSION(:), POINTER :: core_particle_set, particle_set, & 1916 shell_particle_set 1917 TYPE(simpar_type), POINTER :: simpar 1918 TYPE(tmp_variables_type), POINTER :: tmp 1919 TYPE(virial_type), POINTER :: virial 1920 1921 NULLIFY (gci, force_env) 1922 NULLIFY (atomic_kinds, cell, para_env, subsys, local_molecules, local_particles) 1923 NULLIFY (molecule_kinds, molecules, molecule_kind_set, npt) 1924 NULLIFY (core_particles, particles, shell_particles, tmp, old) 1925 NULLIFY (core_particle_set, particle_set, shell_particle_set) 1926 NULLIFY (simpar, virial, itimes) 1927 1928 CALL get_md_env(md_env=md_env, simpar=simpar, force_env=force_env, npt=npt, & 1929 first_time=first_time, para_env=para_env, itimes=itimes) 1930 dt = simpar%dt 1931 infree = 1.0_dp/REAL(simpar%nfree, dp) 1932 gamma1 = simpar%gamma_nph 1933 1934 CALL force_env_get(force_env, subsys=subsys, cell=cell) 1935 1936 CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, & 1937 particles=particles, local_molecules=local_molecules, molecules=molecules, gci=gci, & 1938 molecule_kinds=molecule_kinds, virial=virial) 1939 1940 nparticle_kind = atomic_kinds%n_els 1941 atomic_kind_set => atomic_kinds%els 1942 molecule_kind_set => molecule_kinds%els 1943 1944 nparticle = particles%n_els 1945 particle_set => particles%els 1946 molecule_set => molecules%els 1947 1948 IF (first_time) THEN 1949 CALL virial_evaluate(atomic_kind_set, particle_set, & 1950 local_particles, virial, para_env%group) 1951 END IF 1952 1953 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, & 1954 shell_present=shell_present, shell_adiabatic=shell_adiabatic) 1955 1956 ! Allocate work storage for positions and velocities 1957 CALL allocate_old(old, particle_set, npt) 1958 1959 IF (shell_present) THEN 1960 CALL cp_subsys_get(subsys=subsys, & 1961 shell_particles=shell_particles, core_particles=core_particles) 1962 shell_particle_set => shell_particles%els 1963 nshell = SIZE(shell_particles%els) 1964 IF (shell_adiabatic) THEN 1965 core_particle_set => core_particles%els 1966 END IF 1967 END IF 1968 1969 CALL allocate_tmp(md_env, tmp, nparticle, nshell, shell_adiabatic) 1970 1971 ! perform damping on velocities 1972 CALL damp_v(molecule_kind_set, molecule_set, particle_set, local_molecules, & 1973 gamma1, npt(1, 1), dt, para_env%group) 1974 1975 IF (simpar%constraint) THEN 1976 ! Possibly update the target values 1977 CALL shake_update_targets(gci, local_molecules, molecule_set, & 1978 molecule_kind_set, dt, force_env%root_section) 1979 END IF 1980 1981 ! setting up for ROLL: saving old variables 1982 IF (simpar%constraint) THEN 1983 roll_tol_thrs = simpar%roll_tol 1984 iroll = 1 1985 CALL set(old, atomic_kind_set, particle_set, local_particles, cell, npt, 'F') 1986 CALL getold(gci, local_molecules, molecule_set, & 1987 molecule_kind_set, particle_set, cell) 1988 ELSE 1989 roll_tol_thrs = EPSILON(0.0_dp) 1990 ENDIF 1991 roll_tol = -roll_tol_thrs 1992 1993 SR: DO WHILE (ABS(roll_tol) >= roll_tol_thrs) ! SHAKE-ROLL LOOP 1994 1995 ! perform damping on the barostat momentum 1996 CALL damp_veps(npt(1, 1), gamma1, dt) 1997 1998 IF (simpar%constraint) THEN 1999 CALL set(old, atomic_kind_set, particle_set, local_particles, cell, npt, 'B') 2000 END IF 2001 CALL update_pv(gci, simpar, atomic_kind_set, particle_set, & 2002 local_molecules, molecule_set, molecule_kind_set, & 2003 local_particles, kin, pv_kin, virial, para_env%group) 2004 CALL update_veps(cell, npt, simpar, pv_kin, kin, virial, infree) 2005 2006 ! perform damping on the barostat momentum 2007 CALL damp_veps(npt(1, 1), gamma1, dt) 2008 2009 tmp%arg_r(1) = (0.5_dp*npt(1, 1)%v*dt)* & 2010 (0.5_dp*npt(1, 1)%v*dt) 2011 tmp%poly_r(1) = 1._dp + e2*tmp%arg_r(1) + e4*tmp%arg_r(1)*tmp%arg_r(1) + & 2012 e6*tmp%arg_r(1)**3 + e8*tmp%arg_r(1)**4 2013 2014 aax = npt(1, 1)%v*(1.0_dp + infree) 2015 tmp%arg_v(1) = (0.25_dp*dt*aax)*(0.25_dp*dt*aax) 2016 tmp%poly_v(1) = 1._dp + e2*tmp%arg_v(1) + e4*tmp%arg_v(1)*tmp%arg_v(1) + & 2017 e6*tmp%arg_v(1)**3 + e8*tmp%arg_v(1)**4 2018 2019 aa = npt(1, 1)%v*infree 2020 tmp%arg_v(2) = (0.25_dp*dt*aa)*(0.25_dp*dt*aa) 2021 tmp%poly_v(2) = 1._dp + e2*tmp%arg_v(2) + e4*tmp%arg_v(2)*tmp%arg_v(2) + & 2022 e6*tmp%arg_v(2)**3 + e8*tmp%arg_v(2)**4 2023 tmp%poly_v(3) = 1._dp + e2*tmp%arg_v(2) + e4*tmp%arg_v(2)*tmp%arg_v(2) + & 2024 e6*tmp%arg_v(2)**3 + e8*tmp%arg_v(2)**4 2025 2026 tmp%scale_r(1) = EXP(0.5_dp*dt*npt(1, 1)%v) 2027 tmp%scale_v(1) = EXP(-0.25_dp*dt*aax) 2028 tmp%scale_v(2) = EXP(-0.25_dp*dt*aa) 2029 tmp%scale_v(3) = EXP(-0.25_dp*dt*aa) 2030 2031 ! first half of velocity verlet 2032 CALL vv_first(tmp, atomic_kind_set, local_particles, particle_set, & 2033 core_particle_set, shell_particle_set, nparticle_kind, & 2034 shell_adiabatic, dt) 2035 2036 IF (simpar%variable_dt) CALL variable_timestep(md_env, tmp, dt, simpar, para_env, & 2037 atomic_kind_set, local_particles, particle_set, core_particle_set, & 2038 shell_particle_set, nparticle_kind, shell_adiabatic, npt=npt) 2039 2040 roll_tol = 0._dp 2041 vector_r(:) = 0._dp 2042 vector_v(:) = tmp%scale_v(:)*tmp%poly_v(:) 2043 vector_r(1) = tmp%scale_r(1)*tmp%poly_r(1) 2044 2045 IF (simpar%constraint) CALL shake_roll_control(gci, local_molecules, & 2046 molecule_set, molecule_kind_set, particle_set, tmp%pos, tmp%vel, dt, simpar, & 2047 roll_tol, iroll, vector_r, vector_v, para_env%group, cell=cell, & 2048 local_particles=local_particles) 2049 END DO SR 2050 2051 ! Update h_mat 2052 cell%hmat(1, 1) = cell%hmat(1, 1)*tmp%scale_r(1)*tmp%scale_r(1) 2053 2054 ! Update the inverse 2055 CALL init_cell(cell) 2056 2057 ! Broadcast the new particle positions and deallocate the pos components of temporary 2058 CALL update_dealloc_tmp(tmp, particle_set, shell_particle_set, & 2059 core_particle_set, para_env, shell_adiabatic, pos=.TRUE.) 2060 2061 ! Update forces 2062 CALL force_env_calc_energy_force(force_env) 2063 2064 ! Metadynamics 2065 CALL metadyn_integrator(force_env, itimes, tmp%vel) 2066 2067 ! Velocity Verlet (second part) 2068 CALL vv_second(tmp, atomic_kind_set, local_particles, particle_set, & 2069 core_particle_set, shell_particle_set, nparticle_kind, & 2070 shell_adiabatic, dt) 2071 2072 IF (simpar%constraint) THEN 2073 roll_tol_thrs = simpar%roll_tol 2074 first = .TRUE. 2075 iroll = 1 2076 CALL set(old, atomic_kind_set, particle_set, tmp%vel, local_particles, cell, npt, 'F') 2077 ELSE 2078 roll_tol_thrs = EPSILON(0.0_dp) 2079 ENDIF 2080 roll_tol = -roll_tol_thrs 2081 2082 RR: DO WHILE (ABS(roll_tol) >= roll_tol_thrs) ! RATTLE-ROLL LOOP 2083 roll_tol = 0._dp 2084 IF (simpar%constraint) CALL rattle_roll_setup(old, gci, atomic_kind_set, & 2085 particle_set, local_particles, molecule_kind_set, molecule_set, local_molecules, & 2086 tmp%vel, dt, cell, npt, simpar, virial, vector_v, roll_tol, iroll, infree, first, & 2087 para_env) 2088 ! perform damping on the barostat momentum 2089 CALL damp_veps(npt(1, 1), gamma1, dt) 2090 2091 CALL update_pv(gci, simpar, atomic_kind_set, tmp%vel, particle_set, & 2092 local_molecules, molecule_set, molecule_kind_set, & 2093 local_particles, kin, pv_kin, virial, para_env%group) 2094 CALL update_veps(cell, npt, simpar, pv_kin, kin, virial, infree) 2095 2096 ! perform damping on the barostat momentum 2097 CALL damp_veps(npt(1, 1), gamma1, dt) 2098 2099 END DO RR 2100 2101 ! perform damping on velocities 2102 CALL damp_v(molecule_kind_set, molecule_set, particle_set, local_molecules, & 2103 tmp%vel, gamma1, npt(1, 1), dt, para_env%group) 2104 2105 IF (simpar%annealing) tmp%vel(:, :) = tmp%vel(:, :)*simpar%f_annealing 2106 2107 ! Broadcast the new particle velocities and deallocate temporary 2108 CALL update_dealloc_tmp(tmp, particle_set, shell_particle_set, & 2109 core_particle_set, para_env, shell_adiabatic, vel=.TRUE.) 2110 2111 ! Update constraint virial 2112 IF (simpar%constraint) CALL pv_constraint(gci, local_molecules, & 2113 molecule_set, molecule_kind_set, particle_set, virial, para_env%group) 2114 2115 CALL virial_evaluate(atomic_kind_set, particle_set, & 2116 local_particles, virial, para_env%group) 2117 2118 ! Deallocate old variables 2119 CALL deallocate_old(old) 2120 2121 IF (first_time) THEN 2122 first_time = .FALSE. 2123 CALL set_md_env(md_env, first_time=first_time) 2124 END IF 2125 2126 END SUBROUTINE nph_uniaxial_damped 2127 2128! ************************************************************************************************** 2129!> \brief Velocity Verlet integrator for the NPT ensemble with fully flexible cell 2130!> \param md_env ... 2131!> \param globenv ... 2132!> \par History 2133!> none 2134!> \author CJM 2135! ************************************************************************************************** 2136 SUBROUTINE npt_f(md_env, globenv) 2137 2138 TYPE(md_environment_type), POINTER :: md_env 2139 TYPE(global_environment_type), POINTER :: globenv 2140 2141 CHARACTER(len=*), PARAMETER :: routineN = 'npt_f', routineP = moduleN//':'//routineN 2142 REAL(KIND=dp), PARAMETER :: e2 = 1.0_dp/6.0_dp, e4 = e2/20.0_dp, & 2143 e6 = e4/42.0_dp, e8 = e6/72.0_dp 2144 2145 INTEGER :: i, iroll, j, nparticle, nparticle_kind, & 2146 nshell 2147 INTEGER, POINTER :: itimes 2148 LOGICAL :: first, first_time, shell_adiabatic, & 2149 shell_check_distance, shell_present 2150 REAL(KIND=dp) :: dt, infree, kin, roll_tol, & 2151 roll_tol_thrs, trvg 2152 REAL(KIND=dp), DIMENSION(3) :: vector_r, vector_v 2153 REAL(KIND=dp), DIMENSION(3, 3) :: pv_kin, uh 2154 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds 2155 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 2156 TYPE(barostat_type), POINTER :: barostat 2157 TYPE(cell_type), POINTER :: cell 2158 TYPE(cp_para_env_type), POINTER :: para_env 2159 TYPE(cp_subsys_type), POINTER :: subsys 2160 TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles 2161 TYPE(force_env_type), POINTER :: force_env 2162 TYPE(global_constraint_type), POINTER :: gci 2163 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds 2164 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set 2165 TYPE(molecule_list_type), POINTER :: molecules 2166 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 2167 TYPE(npt_info_type), POINTER :: npt(:, :) 2168 TYPE(old_variables_type), POINTER :: old 2169 TYPE(particle_list_type), POINTER :: core_particles, particles, & 2170 shell_particles 2171 TYPE(particle_type), DIMENSION(:), POINTER :: core_particle_set, particle_set, & 2172 shell_particle_set 2173 TYPE(simpar_type), POINTER :: simpar 2174 TYPE(thermostat_type), POINTER :: thermostat_baro, thermostat_part, & 2175 thermostat_shell 2176 TYPE(tmp_variables_type), POINTER :: tmp 2177 TYPE(virial_type), POINTER :: virial 2178 2179 NULLIFY (gci, thermostat_baro, thermostat_part, thermostat_shell, force_env) 2180 NULLIFY (atomic_kinds, cell, para_env, subsys, local_molecules, local_particles) 2181 NULLIFY (molecule_kinds, molecules, molecule_kind_set, npt, barostat) 2182 NULLIFY (core_particles, particles, shell_particles, tmp, old) 2183 NULLIFY (core_particle_set, particle_set, shell_particle_set) 2184 NULLIFY (simpar, virial, itimes) 2185 2186 CALL get_md_env(md_env=md_env, simpar=simpar, force_env=force_env, & 2187 thermostat_part=thermostat_part, thermostat_baro=thermostat_baro, & 2188 thermostat_shell=thermostat_shell, npt=npt, first_time=first_time, & 2189 para_env=para_env, barostat=barostat, itimes=itimes) 2190 dt = simpar%dt 2191 infree = 1.0_dp/REAL(simpar%nfree, KIND=dp) 2192 2193 CALL force_env_get(force_env, subsys=subsys, cell=cell) 2194 2195 ! Do some checks on coordinates and box 2196 CALL apply_qmmm_walls_reflective(force_env) 2197 2198 CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, & 2199 particles=particles, local_molecules=local_molecules, molecules=molecules, & 2200 gci=gci, molecule_kinds=molecule_kinds, virial=virial) 2201 2202 nparticle_kind = atomic_kinds%n_els 2203 atomic_kind_set => atomic_kinds%els 2204 molecule_kind_set => molecule_kinds%els 2205 2206 nparticle = particles%n_els 2207 particle_set => particles%els 2208 molecule_set => molecules%els 2209 2210 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, & 2211 shell_present=shell_present, shell_adiabatic=shell_adiabatic, & 2212 shell_check_distance=shell_check_distance) 2213 2214 IF (first_time) THEN 2215 CALL virial_evaluate(atomic_kind_set, particle_set, & 2216 local_particles, virial, para_env%group) 2217 END IF 2218 2219 ! Allocate work storage for positions and velocities 2220 CALL allocate_old(old, particle_set, npt) 2221 2222 IF (shell_present) THEN 2223 CALL cp_subsys_get(subsys=subsys, & 2224 shell_particles=shell_particles, core_particles=core_particles) 2225 shell_particle_set => shell_particles%els 2226 nshell = SIZE(shell_particles%els) 2227 IF (shell_adiabatic) THEN 2228 core_particle_set => core_particles%els 2229 END IF 2230 END IF 2231 2232 CALL allocate_tmp(md_env, tmp, nparticle, nshell, shell_adiabatic) 2233 2234 ! Apply Thermostat to Barostat 2235 CALL apply_thermostat_baro(thermostat_baro, npt, para_env%group) 2236 2237 ! Apply Thermostat over the full set of particles 2238 IF (simpar%ensemble /= npe_f_ensemble) THEN 2239 IF (shell_adiabatic) THEN 2240 CALL apply_thermostat_particles(thermostat_part, force_env, molecule_kind_set, molecule_set, & 2241 particle_set, local_molecules, local_particles, para_env%group, shell_adiabatic=shell_adiabatic, & 2242 shell_particle_set=shell_particle_set, core_particle_set=core_particle_set) 2243 ELSE 2244 CALL apply_thermostat_particles(thermostat_part, force_env, molecule_kind_set, molecule_set, & 2245 particle_set, local_molecules, local_particles, para_env%group) 2246 END IF 2247 END IF 2248 2249 ! Apply Thermostat over the core-shell motion 2250 CALL apply_thermostat_shells(thermostat_shell, atomic_kind_set, particle_set, & 2251 local_particles, para_env%group, shell_particle_set=shell_particle_set, & 2252 core_particle_set=core_particle_set) 2253 2254 IF (simpar%constraint) THEN 2255 ! Possibly update the target values 2256 CALL shake_update_targets(gci, local_molecules, molecule_set, & 2257 molecule_kind_set, dt, force_env%root_section) 2258 END IF 2259 2260 ! setting up for ROLL: saving old variables 2261 IF (simpar%constraint) THEN 2262 roll_tol_thrs = simpar%roll_tol 2263 iroll = 1 2264 CALL set(old, atomic_kind_set, particle_set, local_particles, cell, npt, 'F') 2265 CALL getold(gci, local_molecules, molecule_set, & 2266 molecule_kind_set, particle_set, cell) 2267 ELSE 2268 roll_tol_thrs = EPSILON(0.0_dp) 2269 ENDIF 2270 roll_tol = -roll_tol_thrs 2271 2272 SR: DO WHILE (ABS(roll_tol) >= roll_tol_thrs) ! SHAKE-ROLL LOOP 2273 2274 IF (simpar%constraint) THEN 2275 CALL set(old, atomic_kind_set, particle_set, local_particles, cell, npt, 'B') 2276 END IF 2277 CALL update_pv(gci, simpar, atomic_kind_set, particle_set, & 2278 local_molecules, molecule_set, molecule_kind_set, & 2279 local_particles, kin, pv_kin, virial, para_env%group) 2280 CALL update_veps(cell, npt, simpar, pv_kin, kin, virial, infree, & 2281 virial_components=barostat%virial_components) 2282 2283 trvg = npt(1, 1)%v + npt(2, 2)%v + npt(3, 3)%v 2284 ! 2285 ! find eigenvalues and eigenvectors of npt ( :, : )%v 2286 ! 2287 2288 CALL diagonalise(matrix=npt(:, :)%v, mysize=3, & 2289 storageform="UPPER", eigenvalues=tmp%e_val, eigenvectors=tmp%u) 2290 2291 tmp%arg_r(:) = 0.5_dp*tmp%e_val(:)*dt* & 2292 0.5_dp*tmp%e_val(:)*dt 2293 tmp%poly_r = 1.0_dp + e2*tmp%arg_r + e4*tmp%arg_r*tmp%arg_r + & 2294 e6*tmp%arg_r**3 + e8*tmp%arg_r**4 2295 tmp%scale_r(:) = EXP(0.5_dp*dt*tmp%e_val(:)) 2296 2297 tmp%arg_v(:) = 0.25_dp*dt*(tmp%e_val(:) + trvg*infree)* & 2298 0.25_dp*dt*(tmp%e_val(:) + trvg*infree) 2299 tmp%poly_v = 1.0_dp + e2*tmp%arg_v + e4*tmp%arg_v*tmp%arg_v + & 2300 e6*tmp%arg_v**3 + e8*tmp%arg_v**4 2301 tmp%scale_v(:) = EXP(-0.25_dp*dt*(tmp%e_val(:) + trvg*infree)) 2302 2303 CALL vv_first(tmp, atomic_kind_set, local_particles, particle_set, & 2304 core_particle_set, shell_particle_set, nparticle_kind, & 2305 shell_adiabatic, dt, u=tmp%u) 2306 2307 IF (simpar%variable_dt) CALL variable_timestep(md_env, tmp, dt, simpar, para_env, & 2308 atomic_kind_set, local_particles, particle_set, core_particle_set, & 2309 shell_particle_set, nparticle_kind, shell_adiabatic, npt=npt) 2310 2311 roll_tol = 0.0_dp 2312 vector_r = tmp%scale_r*tmp%poly_r 2313 vector_v = tmp%scale_v*tmp%poly_v 2314 2315 IF (simpar%constraint) CALL shake_roll_control(gci, local_molecules, & 2316 molecule_set, molecule_kind_set, particle_set, tmp%pos, tmp%vel, dt, & 2317 simpar, roll_tol, iroll, vector_r, vector_v, & 2318 para_env%group, u=tmp%u, cell=cell, & 2319 local_particles=local_particles) 2320 END DO SR 2321 2322 ! Update h_mat 2323 uh = MATMUL_3X3(TRANSPOSE_3D(tmp%u), cell%hmat) 2324 2325 DO i = 1, 3 2326 DO j = 1, 3 2327 uh(i, j) = uh(i, j)*tmp%scale_r(i)*tmp%scale_r(i) 2328 END DO 2329 END DO 2330 2331 cell%hmat = MATMUL_3x3(tmp%u, uh) 2332 ! Update the inverse 2333 CALL init_cell(cell) 2334 2335 ! Broadcast the new particle positions and deallocate the pos components of temporary 2336 CALL update_dealloc_tmp(tmp, particle_set, shell_particle_set, & 2337 core_particle_set, para_env, shell_adiabatic, pos=.TRUE.) 2338 2339 IF (shell_adiabatic .AND. shell_check_distance) THEN 2340 CALL optimize_shell_core(force_env, particle_set, & 2341 shell_particle_set, core_particle_set, globenv, tmp=tmp, check=.TRUE.) 2342 END IF 2343 2344 ! Update forces 2345 CALL force_env_calc_energy_force(force_env) 2346 2347 ! Metadynamics 2348 CALL metadyn_integrator(force_env, itimes, tmp%vel) 2349 2350 ! Velocity Verlet (second part) 2351 CALL vv_second(tmp, atomic_kind_set, local_particles, particle_set, & 2352 core_particle_set, shell_particle_set, nparticle_kind, & 2353 shell_adiabatic, dt, tmp%u) 2354 2355 IF (simpar%constraint) THEN 2356 roll_tol_thrs = simpar%roll_tol 2357 first = .TRUE. 2358 iroll = 1 2359 CALL set(old, atomic_kind_set, particle_set, tmp%vel, local_particles, cell, npt, 'F') 2360 ELSE 2361 roll_tol_thrs = EPSILON(0.0_dp) 2362 ENDIF 2363 roll_tol = -roll_tol_thrs 2364 2365 RR: DO WHILE (ABS(roll_tol) >= roll_tol_thrs) ! RATTLE-ROLL LOOP 2366 roll_tol = 0.0_dp 2367 IF (simpar%constraint) CALL rattle_roll_setup(old, gci, atomic_kind_set, & 2368 particle_set, local_particles, molecule_kind_set, molecule_set, & 2369 local_molecules, tmp%vel, dt, cell, npt, simpar, virial, vector_v, & 2370 roll_tol, iroll, infree, first, para_env, u=tmp%u) 2371 2372 CALL update_pv(gci, simpar, atomic_kind_set, tmp%vel, particle_set, & 2373 local_molecules, molecule_set, molecule_kind_set, & 2374 local_particles, kin, pv_kin, virial, para_env%group) 2375 CALL update_veps(cell, npt, simpar, pv_kin, kin, virial, infree, & 2376 virial_components=barostat%virial_components) 2377 END DO RR 2378 2379 ! Apply Thermostat over the full set of particles 2380 IF (simpar%ensemble /= npe_f_ensemble) THEN 2381 IF (shell_adiabatic) THEN 2382 CALL apply_thermostat_particles(thermostat_part, force_env, molecule_kind_set, molecule_set, & 2383 particle_set, local_molecules, local_particles, para_env%group, shell_adiabatic=shell_adiabatic, & 2384 vel=tmp%vel, shell_vel=tmp%shell_vel, core_vel=tmp%core_vel) 2385 2386 ELSE 2387 CALL apply_thermostat_particles(thermostat_part, force_env, molecule_kind_set, molecule_set, & 2388 particle_set, local_molecules, local_particles, para_env%group, vel=tmp%vel) 2389 END IF 2390 END IF 2391 2392 ! Apply Thermostat over the core-shell motion 2393 IF (ASSOCIATED(thermostat_shell)) THEN 2394 CALL apply_thermostat_shells(thermostat_shell, atomic_kind_set, particle_set, & 2395 local_particles, para_env%group, vel=tmp%vel, shell_vel=tmp%shell_vel, & 2396 core_vel=tmp%core_vel) 2397 END IF 2398 2399 ! Apply Thermostat to Barostat 2400 CALL apply_thermostat_baro(thermostat_baro, npt, para_env%group) 2401 2402 ! Annealing of particle velocities is only possible when no thermostat is active 2403 IF (simpar%ensemble == npe_f_ensemble .AND. simpar%annealing) THEN 2404 tmp%vel(:, :) = tmp%vel(:, :)*simpar%f_annealing 2405 IF (shell_adiabatic) THEN 2406 CALL shell_scale_comv(atomic_kind_set, local_particles, particle_set, & 2407 tmp%vel, tmp%shell_vel, tmp%core_vel) 2408 END IF 2409 END IF 2410 ! Annealing of CELL velocities is only possible when no thermostat is active 2411 IF (simpar%ensemble == npe_f_ensemble .AND. simpar%annealing_cell) THEN 2412 npt(:, :)%v = npt(:, :)%v*simpar%f_annealing_cell 2413 END IF 2414 2415 ! Broadcast the new particle velocities and deallocate temporary 2416 CALL update_dealloc_tmp(tmp, particle_set, shell_particle_set, & 2417 core_particle_set, para_env, shell_adiabatic, vel=.TRUE.) 2418 2419 ! Update constraint virial 2420 IF (simpar%constraint) & 2421 CALL pv_constraint(gci, local_molecules, molecule_set, & 2422 molecule_kind_set, particle_set, virial, para_env%group) 2423 2424 CALL virial_evaluate(atomic_kind_set, particle_set, & 2425 local_particles, virial, para_env%group) 2426 2427 ! Deallocate old variables 2428 CALL deallocate_old(old) 2429 2430 IF (first_time) THEN 2431 first_time = .FALSE. 2432 CALL set_md_env(md_env, first_time=first_time) 2433 END IF 2434 2435 END SUBROUTINE npt_f 2436 2437! ************************************************************************************************** 2438!> \brief RESPA integrator for nve ensemble for particle positions & momenta 2439!> \param md_env ... 2440!> \author FS 2441! ************************************************************************************************** 2442 SUBROUTINE nve_respa(md_env) 2443 2444 TYPE(md_environment_type), POINTER :: md_env 2445 2446 CHARACTER(len=*), PARAMETER :: routineN = 'nve_respa', routineP = moduleN//':'//routineN 2447 2448 INTEGER :: i_step, iparticle, iparticle_kind, & 2449 iparticle_local, n_time_steps, & 2450 nparticle, nparticle_kind, & 2451 nparticle_local 2452 INTEGER, POINTER :: itimes 2453 REAL(KIND=dp) :: dm, dt, mass 2454 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: pos, vel 2455 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds 2456 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 2457 TYPE(atomic_kind_type), POINTER :: atomic_kind 2458 TYPE(cell_type), POINTER :: cell 2459 TYPE(cp_para_env_type), POINTER :: para_env 2460 TYPE(cp_subsys_type), POINTER :: subsys, subsys_respa 2461 TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles 2462 TYPE(force_env_type), POINTER :: force_env 2463 TYPE(global_constraint_type), POINTER :: gci 2464 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds 2465 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set 2466 TYPE(molecule_list_type), POINTER :: molecules 2467 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 2468 TYPE(particle_list_type), POINTER :: particles, particles_respa 2469 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set, particle_set_respa 2470 TYPE(simpar_type), POINTER :: simpar 2471 2472 NULLIFY (para_env, cell, subsys_respa, particles_respa, particle_set_respa, gci, force_env, atomic_kinds) 2473 NULLIFY (atomic_kind_set, simpar, subsys, particles, particle_set) 2474 NULLIFY (local_molecules, molecule_kinds, molecules, molecule_kind_set, local_particles, itimes) 2475 CALL get_md_env(md_env=md_env, simpar=simpar, force_env=force_env, & 2476 para_env=para_env, itimes=itimes) 2477 dt = simpar%dt 2478 2479 n_time_steps = simpar%n_time_steps 2480 2481 CALL force_env_get(force_env, subsys=subsys, cell=cell) 2482 CALL force_env_get(force_env%sub_force_env(1)%force_env, subsys=subsys_respa) 2483 2484 ! Do some checks on coordinates and box 2485 CALL apply_qmmm_walls_reflective(force_env) 2486 2487 CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, & 2488 particles=particles, local_molecules=local_molecules, molecules=molecules, & 2489 gci=gci, molecule_kinds=molecule_kinds) 2490 2491 CALL cp_subsys_get(subsys=subsys_respa, particles=particles_respa) 2492 particle_set_respa => particles_respa%els 2493 2494 nparticle_kind = atomic_kinds%n_els 2495 atomic_kind_set => atomic_kinds%els 2496 molecule_kind_set => molecule_kinds%els 2497 2498 nparticle = particles%n_els 2499 particle_set => particles%els 2500 molecule_set => molecules%els 2501 2502 ! Allocate work storage for positions and velocities 2503 ALLOCATE (pos(3, nparticle)) 2504 ALLOCATE (vel(3, nparticle)) 2505 vel(:, :) = 0.0_dp 2506 2507 IF (simpar%constraint) CALL getold(gci, local_molecules, molecule_set, & 2508 molecule_kind_set, particle_set, cell) 2509 2510 ! Multiple time step (first part) 2511 DO iparticle_kind = 1, nparticle_kind 2512 atomic_kind => atomic_kind_set(iparticle_kind) 2513 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass) 2514 dm = 0.5_dp*dt/mass 2515 nparticle_local = local_particles%n_el(iparticle_kind) 2516 DO iparticle_local = 1, nparticle_local 2517 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) 2518 vel(:, iparticle) = particle_set(iparticle)%v(:) + & 2519 dm*(particle_set(iparticle)%f(:) - & 2520 particle_set_respa(iparticle)%f(:)) 2521 END DO 2522 END DO 2523 2524 ! Velocity Verlet (first part) 2525 DO i_step = 1, n_time_steps 2526 pos(:, :) = 0.0_dp 2527 DO iparticle_kind = 1, nparticle_kind 2528 atomic_kind => atomic_kind_set(iparticle_kind) 2529 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass) 2530 dm = 0.5_dp*dt/(n_time_steps*mass) 2531 nparticle_local = local_particles%n_el(iparticle_kind) 2532 DO iparticle_local = 1, nparticle_local 2533 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) 2534 vel(:, iparticle) = vel(:, iparticle) + & 2535 dm*particle_set_respa(iparticle)%f(:) 2536 pos(:, iparticle) = particle_set(iparticle)%r(:) + & 2537 (dt/n_time_steps)*vel(:, iparticle) 2538 END DO 2539 END DO 2540 2541 IF (simpar%constraint) THEN 2542 ! Possibly update the target values 2543 CALL shake_update_targets(gci, local_molecules, molecule_set, & 2544 molecule_kind_set, dt, force_env%root_section) 2545 2546 CALL shake_control(gci, local_molecules, molecule_set, & 2547 molecule_kind_set, particle_set, pos, vel, dt, simpar%shake_tol, & 2548 simpar%info_constraint, simpar%lagrange_multipliers, simpar%dump_lm, cell, & 2549 para_env%group, local_particles) 2550 END IF 2551 2552 ! Broadcast the new particle positions 2553 CALL update_particle_set(particle_set, para_env%group, pos=pos) 2554 DO iparticle = 1, SIZE(particle_set) 2555 particle_set_respa(iparticle)%r = particle_set(iparticle)%r 2556 END DO 2557 2558 ! Update forces 2559 CALL force_env_calc_energy_force(force_env%sub_force_env(1)%force_env) 2560 2561 ! Metadynamics 2562 CALL metadyn_integrator(force_env, itimes, vel) 2563 2564 ! Velocity Verlet (second part) 2565 DO iparticle_kind = 1, nparticle_kind 2566 atomic_kind => atomic_kind_set(iparticle_kind) 2567 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass) 2568 dm = 0.5_dp*dt/(n_time_steps*mass) 2569 nparticle_local = local_particles%n_el(iparticle_kind) 2570 DO iparticle_local = 1, nparticle_local 2571 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) 2572 vel(1, iparticle) = vel(1, iparticle) + dm*particle_set_respa(iparticle)%f(1) 2573 vel(2, iparticle) = vel(2, iparticle) + dm*particle_set_respa(iparticle)%f(2) 2574 vel(3, iparticle) = vel(3, iparticle) + dm*particle_set_respa(iparticle)%f(3) 2575 END DO 2576 END DO 2577 2578 IF (simpar%constraint) CALL rattle_control(gci, local_molecules, molecule_set, & 2579 molecule_kind_set, particle_set, vel, dt, simpar%shake_tol, & 2580 simpar%info_constraint, simpar%lagrange_multipliers, & 2581 simpar%dump_lm, cell, para_env%group, local_particles) 2582 2583 IF (simpar%annealing) vel(:, :) = vel(:, :)*simpar%f_annealing 2584 END DO 2585 DEALLOCATE (pos) 2586 2587 ! Multiple time step (second part) 2588 ! Compute forces for respa force_env 2589 CALL force_env_calc_energy_force(force_env) 2590 2591 ! Metadynamics 2592 CALL metadyn_integrator(force_env, itimes, vel) 2593 2594 DO iparticle_kind = 1, nparticle_kind 2595 atomic_kind => atomic_kind_set(iparticle_kind) 2596 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass) 2597 dm = 0.5_dp*dt/mass 2598 nparticle_local = local_particles%n_el(iparticle_kind) 2599 DO iparticle_local = 1, nparticle_local 2600 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) 2601 vel(1, iparticle) = vel(1, iparticle) + dm*(particle_set(iparticle)%f(1) - particle_set_respa(iparticle)%f(1)) 2602 vel(2, iparticle) = vel(2, iparticle) + dm*(particle_set(iparticle)%f(2) - particle_set_respa(iparticle)%f(2)) 2603 vel(3, iparticle) = vel(3, iparticle) + dm*(particle_set(iparticle)%f(3) - particle_set_respa(iparticle)%f(3)) 2604 END DO 2605 END DO 2606 2607 ! Broadcast the new particle velocities 2608 CALL update_particle_set(particle_set, para_env%group, vel=vel) 2609 2610 DEALLOCATE (vel) 2611 2612 END SUBROUTINE nve_respa 2613 2614END MODULE integrator 2615