comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-29 06:35:28 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Motion: Simplify NPT_IA ensemble

Browse source
This commit is contained in:
Ole Schütt 2022-03-13 22:01:19 +01:00 • committed by Ole Schütt
parent a564af00f6
commit ff401dcee0
15 changed files with 46 additions and 1470 deletions
Download patch file Download diff file Expand all files Collapse all files

360
src/motion/integrator.F
View file

@ -59,11 +59,12 @@ MODULE integrator
USE global_types, ONLY: global_environment_type
USE input_constants, ONLY: ehrenfest,&
npe_f_ensemble,&
npe_i_ensemble
npe_i_ensemble,&
npt_ia_ensemble
USE integrator_utils, ONLY: &
allocate_old, allocate_tmp, damp_v, damp_veps, deallocate_old, get_s_ds, &
old_variables_type, rattle_roll_setup, set, tmp_variables_type, update_dealloc_tmp, &
update_pv, update_veps, variable_timestep, vv_first, vv_first_frozen, vv_second
update_pv, update_veps, variable_timestep, vv_first, vv_second
USE kinds, ONLY: dp
USE mathlib, ONLY: matmul_3x3,&
transpose_3d
@ -109,7 +110,7 @@ MODULE integrator
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'integrator'
PUBLIC :: isokin, langevin, nve, nvt, npt_i, npt_ia, npt_f, nve_respa
PUBLIC :: isokin, langevin, nve, nvt, npt_i, npt_f, nve_respa
PUBLIC :: nph_uniaxial_damped, nph_uniaxial, nvt_adiabatic, reftraj
CONTAINS
@ -1152,7 +1153,7 @@ CONTAINS
END SUBROUTINE nvt
! **************************************************************************************************
!> \brief npt_ia integrator for particle positions & momenta
!> \brief npt_i integrator for particle positions & momenta
!> isotropic box changes
!> \param md_env ...
!> \param globenv ...
@ -1160,7 +1161,7 @@ CONTAINS
!> none
!> \author CJM
! **************************************************************************************************
SUBROUTINE npt_ia(md_env, globenv)
SUBROUTINE npt_i(md_env, globenv)
TYPE(md_environment_type), POINTER :: md_env
TYPE(global_environment_type), POINTER :: globenv
@ -1186,7 +1187,8 @@ CONTAINS
TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles
TYPE(force_env_type), POINTER :: force_env
TYPE(global_constraint_type), POINTER :: gci
TYPE(local_fixd_constraint_type), POINTER :: lfixd_list(:)
TYPE(local_fixd_constraint_type), DIMENSION(:), &
POINTER :: lfixd_list
TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
TYPE(molecule_list_type), POINTER :: molecules
@ -1208,7 +1210,7 @@ CONTAINS
NULLIFY (molecule_kinds, molecules, molecule_kind_set, npt)
NULLIFY (core_particles, particles, shell_particles, tmp, old)
NULLIFY (core_particle_set, particle_set, shell_particle_set)
NULLIFY (simpar, virial, rand, itimes)
NULLIFY (simpar, virial, rand, itimes, lfixd_list)
CALL get_md_env(md_env=md_env, simpar=simpar, force_env=force_env, &
thermostat_part=thermostat_part, thermostat_baro=thermostat_baro, &
@ -1247,336 +1249,6 @@ CONTAINS
! Allocate work storage for positions and velocities
CALL allocate_old(old, particle_set, npt)
IF (ASSOCIATED(force_env%meta_env)) THEN
! Allocate random number for Langevin Thermostat acting on COLVARS
IF (force_env%meta_env%langevin) THEN
ALLOCATE (rand(force_env%meta_env%n_colvar))
rand(:) = 0.0_dp
ENDIF
ENDIF
IF (shell_present) THEN
CALL cp_subsys_get(subsys=subsys, &
shell_particles=shell_particles, core_particles=core_particles)
shell_particle_set => shell_particles%els
nshell = SIZE(shell_particles%els)
IF (shell_adiabatic) THEN
core_particle_set => core_particles%els
END IF
END IF
CALL allocate_tmp(md_env, tmp, nparticle, nshell, shell_adiabatic)
! Initialize eps_0 the first time through
IF (first_time) eps_0 = npt(1, 1)%eps
! Apply thermostat to barostat
CALL apply_thermostat_baro(thermostat_baro, npt, para_env%group)
! Apply Thermostat over the full set of particles
IF (simpar%ensemble /= npe_i_ensemble) THEN
IF (shell_adiabatic) THEN
CALL apply_thermostat_particles(thermostat_part, force_env, molecule_kind_set, molecule_set, &
particle_set, local_molecules, local_particles, para_env%group, shell_adiabatic=shell_adiabatic, &
shell_particle_set=shell_particle_set, core_particle_set=core_particle_set)
ELSE
CALL apply_thermostat_particles(thermostat_part, force_env, molecule_kind_set, molecule_set, &
particle_set, local_molecules, local_particles, para_env%group)
END IF
END IF
! Apply Thermostat over the core-shell motion
CALL apply_thermostat_shells(thermostat_shell, atomic_kind_set, particle_set, &
local_particles, para_env%group, shell_particle_set=shell_particle_set, &
core_particle_set=core_particle_set)
IF (simpar%constraint) THEN
! Possibly update the target values
CALL shake_update_targets(gci, local_molecules, molecule_set, &
molecule_kind_set, dt, force_env%root_section)
END IF
! setting up for ROLL: saving old variables
IF (simpar%constraint) THEN
roll_tol_thrs = simpar%roll_tol
iroll = 1
CALL set(old, atomic_kind_set, particle_set, local_particles, cell, npt, 'F')
CALL getold(gci, local_molecules, molecule_set, &
molecule_kind_set, particle_set, cell)
ELSE
roll_tol_thrs = EPSILON(0.0_dp)
ENDIF
roll_tol = -roll_tol_thrs
! *** Velocity Verlet for Langeving *** v(t)--> v(t+1/2)
IF (ASSOCIATED(force_env%meta_env)) THEN
IF (force_env%meta_env%langevin) THEN
DO ivar = 1, force_env%meta_env%n_colvar
rand(ivar) = force_env%meta_env%rng(ivar)%next()
ENDDO
CALL metadyn_velocities_colvar(force_env, rand)
ENDIF
ENDIF
SR: DO WHILE (ABS(roll_tol) >= roll_tol_thrs) ! SHAKE-ROLL LOOP
IF (simpar%constraint) THEN
CALL set(old, atomic_kind_set, particle_set, local_particles, cell, npt, 'B')
END IF
CALL update_pv(gci, simpar, atomic_kind_set, particle_set, &
local_molecules, molecule_set, molecule_kind_set, &
local_particles, kin, pv_kin, virial, para_env%group)
CALL update_veps(cell, npt, simpar, pv_kin, kin, virial, infree)
tmp%arg_r(1) = (0.5_dp*npt(1, 1)%v*dt)* &
(0.5_dp*npt(1, 1)%v*dt)
tmp%poly_r(1:3) = 1.0_dp + e2*tmp%arg_r(1) + e4*tmp%arg_r(1)*tmp%arg_r(1) + &
e6*tmp%arg_r(1)**3 + e8*tmp%arg_r(1)**4
tmp%arg_v(1) = (0.25_dp*npt(1, 1)%v*dt* &
(1.0_dp + 3.0_dp*infree))*(0.25_dp*npt(1, 1)%v* &
dt*(1.0_dp + 3.0_dp*infree))
tmp%poly_v(1:3) = 1.0_dp + e2*tmp%arg_v(1) + e4*tmp%arg_v(1)*tmp%arg_v(1) + &
e6*tmp%arg_v(1)**3 + e8*tmp%arg_v(1)**4
tmp%scale_r(1:3) = EXP(0.5_dp*dt*npt(1, 1)%v)
tmp%scale_v(1:3) = EXP(-0.25_dp*dt*npt(1, 1)%v* &
(1.0_dp + 3.0_dp*infree))
NULLIFY (lfixd_list)
CALL create_local_fixd_list(lfixd_list, nkind, molecule_kind_set, local_particles)
! first half of velocity verlet
CALL vv_first_frozen(tmp, atomic_kind_set, local_particles, particle_set, &
core_particle_set, shell_particle_set, nparticle_kind, &
shell_adiabatic, dt, lfixd_list)
CALL release_local_fixd_list(lfixd_list)
IF (simpar%variable_dt) CALL variable_timestep(md_env, tmp, dt, simpar, para_env, &
atomic_kind_set, local_particles, particle_set, core_particle_set, &
shell_particle_set, nparticle_kind, shell_adiabatic, npt=npt)
roll_tol = 0.0_dp
vector_r(:) = tmp%scale_r(:)*tmp%poly_r(:)
vector_v(:) = tmp%scale_v(:)*tmp%poly_v(:)
IF (simpar%constraint) CALL shake_roll_control(gci, local_molecules, &
molecule_set, molecule_kind_set, particle_set, tmp%pos, tmp%vel, dt, simpar, &
roll_tol, iroll, vector_r, vector_v, para_env%group, cell=cell, &
local_particles=local_particles)
END DO SR
! Update eps:
npt(:, :)%eps = npt(:, :)%eps + dt*npt(:, :)%v
! Update h_mat
cell%hmat(:, :) = cell%hmat(:, :)*EXP(npt(1, 1)%eps - eps_0)
eps_0 = npt(1, 1)%eps
! Update the inverse
CALL init_cell(cell)
! Broadcast the new particle positions and deallocate the pos components of temporary
CALL update_dealloc_tmp(tmp, particle_set, shell_particle_set, &
core_particle_set, para_env, shell_adiabatic, pos=.TRUE.)
IF (shell_adiabatic .AND. shell_check_distance) THEN
CALL optimize_shell_core(force_env, particle_set, &
shell_particle_set, core_particle_set, globenv, tmp=tmp, check=.TRUE.)
END IF
! Update forces
CALL force_env_calc_energy_force(force_env)
! Metadynamics
CALL metadyn_integrator(force_env, itimes, tmp%vel, rand=rand)
! Velocity Verlet (second part)
CALL vv_second(tmp, atomic_kind_set, local_particles, particle_set, &
core_particle_set, shell_particle_set, nparticle_kind, &
shell_adiabatic, dt)
IF (simpar%constraint) THEN
roll_tol_thrs = simpar%roll_tol
first = .TRUE.
iroll = 1
CALL set(old, atomic_kind_set, particle_set, tmp%vel, local_particles, cell, npt, 'F')
ELSE
roll_tol_thrs = EPSILON(0.0_dp)
ENDIF
roll_tol = -roll_tol_thrs
RR: DO WHILE (ABS(roll_tol) >= roll_tol_thrs) ! RATTLE-ROLL LOOP
roll_tol = 0.0_dp
IF (simpar%constraint) CALL rattle_roll_setup(old, gci, atomic_kind_set, &
particle_set, local_particles, molecule_kind_set, molecule_set, &
local_molecules, tmp%vel, dt, cell, npt, simpar, virial, vector_v, &
roll_tol, iroll, infree, first, para_env)
CALL update_pv(gci, simpar, atomic_kind_set, tmp%vel, particle_set, &
local_molecules, molecule_set, molecule_kind_set, &
local_particles, kin, pv_kin, virial, para_env%group)
CALL update_veps(cell, npt, simpar, pv_kin, kin, virial, infree)
END DO RR
! Apply Thermostat over the full set of particles
IF (simpar%ensemble /= npe_i_ensemble) THEN
IF (shell_adiabatic) THEN
CALL apply_thermostat_particles(thermostat_part, force_env, molecule_kind_set, molecule_set, &
particle_set, local_molecules, local_particles, para_env%group, shell_adiabatic=shell_adiabatic, &
vel=tmp%vel, shell_vel=tmp%shell_vel, core_vel=tmp%core_vel)
ELSE
CALL apply_thermostat_particles(thermostat_part, force_env, molecule_kind_set, molecule_set, &
particle_set, local_molecules, local_particles, para_env%group, vel=tmp%vel)
END IF
END IF
! Apply Thermostat over the core-shell motion
IF (ASSOCIATED(thermostat_shell)) THEN
CALL apply_thermostat_shells(thermostat_shell, atomic_kind_set, particle_set, &
local_particles, para_env%group, vel=tmp%vel, shell_vel=tmp%shell_vel, &
core_vel=tmp%core_vel)
END IF
! Apply Thermostat to Barostat
CALL apply_thermostat_baro(thermostat_baro, npt, para_env%group)
! Annealing of particle velocities is only possible when no thermostat is active
IF (simpar%ensemble == npe_i_ensemble .AND. simpar%annealing) THEN
tmp%vel(:, :) = tmp%vel(:, :)*simpar%f_annealing
IF (shell_adiabatic) THEN
CALL shell_scale_comv(atomic_kind_set, local_particles, particle_set, &
tmp%vel, tmp%shell_vel, tmp%core_vel)
END IF
END IF
! Annealing of CELL velocities is only possible when no thermostat is active
IF (simpar%ensemble == npe_i_ensemble .AND. simpar%annealing_cell) THEN
npt(1, 1)%v = npt(1, 1)%v*simpar%f_annealing_cell
END IF
! Broadcast the new particle velocities and deallocate temporary
CALL update_dealloc_tmp(tmp, particle_set, shell_particle_set, &
core_particle_set, para_env, shell_adiabatic, vel=.TRUE.)
! Update constraint virial
IF (simpar%constraint) CALL pv_constraint(gci, local_molecules, &
molecule_set, molecule_kind_set, particle_set, virial, para_env%group)
CALL virial_evaluate(atomic_kind_set, particle_set, &
local_particles, virial, para_env%group)
! Deallocate old variables
CALL deallocate_old(old)
IF (ASSOCIATED(force_env%meta_env)) THEN
IF (force_env%meta_env%langevin) THEN
DEALLOCATE (rand)
ENDIF
ENDIF
IF (first_time) THEN
first_time = .FALSE.
CALL set_md_env(md_env, first_time=first_time)
END IF
END SUBROUTINE npt_ia
! **************************************************************************************************
!> \brief npt_i integrator for particle positions & momenta
!> isotropic box changes
!> \param md_env ...
!> \param globenv ...
!> \par History
!> none
!> \author CJM
! **************************************************************************************************
SUBROUTINE npt_i(md_env, globenv)
TYPE(md_environment_type), POINTER :: md_env
TYPE(global_environment_type), POINTER :: globenv
REAL(KIND=dp), PARAMETER :: e2 = 1.0_dp/6.0_dp, e4 = e2/20.0_dp, &
e6 = e4/42.0_dp, e8 = e6/72.0_dp
INTEGER :: iroll, ivar, nparticle, nparticle_kind, &
nshell
INTEGER, POINTER :: itimes
LOGICAL :: first, first_time, shell_adiabatic, &
shell_check_distance, shell_present
REAL(KIND=dp) :: dt, infree, kin, roll_tol, roll_tol_thrs
REAL(KIND=dp), DIMENSION(3) :: vector_r, vector_v
REAL(KIND=dp), DIMENSION(3, 3) :: pv_kin
REAL(KIND=dp), DIMENSION(:), POINTER :: rand
REAL(KIND=dp), SAVE :: eps_0
TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cell_type), POINTER :: cell
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(cp_subsys_type), POINTER :: subsys
TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles
TYPE(force_env_type), POINTER :: force_env
TYPE(global_constraint_type), POINTER :: gci
TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
TYPE(molecule_list_type), POINTER :: molecules
TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
TYPE(npt_info_type), POINTER :: npt(:, :)
TYPE(old_variables_type), POINTER :: old
TYPE(particle_list_type), POINTER :: core_particles, particles, &
shell_particles
TYPE(particle_type), DIMENSION(:), POINTER :: core_particle_set, particle_set, &
shell_particle_set
TYPE(simpar_type), POINTER :: simpar
TYPE(thermostat_type), POINTER :: thermostat_baro, thermostat_part, &
thermostat_shell
TYPE(tmp_variables_type), POINTER :: tmp
TYPE(virial_type), POINTER :: virial
NULLIFY (gci, thermostat_baro, thermostat_part, thermostat_shell, force_env)
NULLIFY (atomic_kinds, cell, para_env, subsys, local_molecules, local_particles)
NULLIFY (molecule_kinds, molecules, molecule_kind_set, npt)
NULLIFY (core_particles, particles, shell_particles, tmp, old)
NULLIFY (core_particle_set, particle_set, shell_particle_set)
NULLIFY (simpar, virial, rand, itimes)
CALL get_md_env(md_env=md_env, simpar=simpar, force_env=force_env, &
thermostat_part=thermostat_part, thermostat_baro=thermostat_baro, &
thermostat_shell=thermostat_shell, npt=npt, first_time=first_time, &
para_env=para_env, itimes=itimes)
dt = simpar%dt
infree = 1.0_dp/REAL(simpar%nfree, KIND=dp)
CALL force_env_get(force_env=force_env, subsys=subsys, cell=cell)
! Do some checks on coordinates and box
CALL apply_qmmm_walls_reflective(force_env)
CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, &
particles=particles, local_molecules=local_molecules, molecules=molecules, &
gci=gci, molecule_kinds=molecule_kinds, virial=virial)
nparticle_kind = atomic_kinds%n_els
atomic_kind_set => atomic_kinds%els
molecule_kind_set => molecule_kinds%els
nparticle = particles%n_els
particle_set => particles%els
molecule_set => molecules%els
CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
shell_present=shell_present, shell_adiabatic=shell_adiabatic, &
shell_check_distance=shell_check_distance)
IF (first_time) THEN
CALL virial_evaluate(atomic_kind_set, particle_set, &
local_particles, virial, para_env%group)
END IF
! Allocate work storage for positions and velocities
CALL allocate_old(old, particle_set, npt)
IF (ASSOCIATED(force_env%meta_env)) THEN
! Allocate random number for Langevin Thermostat acting on COLVARS
IF (force_env%meta_env%langevin) THEN
@ -1676,9 +1348,17 @@ CONTAINS
(1.0_dp + 3.0_dp*infree))
! first half of velocity verlet
CALL vv_first(tmp, atomic_kind_set, local_particles, particle_set, &
core_particle_set, shell_particle_set, nparticle_kind, &
shell_adiabatic, dt)
IF (simpar%ensemble == npt_ia_ensemble) THEN
CALL create_local_fixd_list(lfixd_list, nkind, molecule_kind_set, local_particles)
CALL vv_first(tmp, atomic_kind_set, local_particles, particle_set, &
core_particle_set, shell_particle_set, nparticle_kind, &
shell_adiabatic, dt, lfixd_list=lfixd_list)
CALL release_local_fixd_list(lfixd_list)
ELSE
CALL vv_first(tmp, atomic_kind_set, local_particles, particle_set, &
core_particle_set, shell_particle_set, nparticle_kind, &
shell_adiabatic, dt)
END IF
IF (simpar%variable_dt) CALL variable_timestep(md_env, tmp, dt, simpar, para_env, &
atomic_kind_set, local_particles, particle_set, core_particle_set, &

245
src/motion/integrator_utils.F
View file

@ -101,7 +101,7 @@ MODULE integrator_utils
PUBLIC :: old_variables_type, tmp_variables_type, allocate_old, deallocate_old, &
allocate_tmp, update_dealloc_tmp, damp_v, set, update_pv, get_s_ds, &
rattle_roll_setup, damp_veps, update_veps, vv_first, vv_first_frozen, &
rattle_roll_setup, damp_veps, update_veps, vv_first, &
vv_second, variable_timestep
CONTAINS
@ -1114,13 +1114,14 @@ CONTAINS
!> \param shell_adiabatic ...
!> \param dt ...
!> \param u ...
!> \param lfixd_list ...
!> \par History
!> none
!> \author MI (February 2008)
! **************************************************************************************************
SUBROUTINE vv_first(tmp, atomic_kind_set, local_particles, particle_set, &
core_particle_set, shell_particle_set, nparticle_kind, &
shell_adiabatic, dt, u)
shell_adiabatic, dt, u, lfixd_list)
TYPE(tmp_variables_type), POINTER :: tmp
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
@ -1131,6 +1132,8 @@ CONTAINS
LOGICAL, INTENT(IN) :: shell_adiabatic
REAL(KIND=dp) :: dt
REAL(KIND=dp), DIMENSION(3, 3), OPTIONAL :: u
TYPE(local_fixd_constraint_type), DIMENSION(:), &
OPTIONAL :: lfixd_list
INTEGER :: iparticle, iparticle_kind, &
iparticle_local, nparticle_local, &
@ -1288,221 +1291,13 @@ CONTAINS
particle_set(iparticle)%r(3)*tmp%scale_r(3)*tmp%scale_r(3) + &
tmp%scale_r(3)*tmp%poly_r(3)*dt*tmp%vel(3, iparticle)
tmp%max_vel = MAX(tmp%max_vel, &
ABS(tmp%vel(1, iparticle)), ABS(tmp%vel(2, iparticle)), ABS(tmp%vel(3, iparticle)))
tmp%max_dr = MAX(tmp%max_dr, &
ABS(particle_set(iparticle)%r(1) - tmp%pos(1, iparticle)), &
ABS(particle_set(iparticle)%r(2) - tmp%pos(2, iparticle)), &
ABS(particle_set(iparticle)%r(3) - tmp%pos(3, iparticle)))
tmp%max_dvel = MAX(tmp%max_dvel, &
ABS(particle_set(iparticle)%v(1) - tmp%vel(1, iparticle)), &
ABS(particle_set(iparticle)%v(2) - tmp%vel(2, iparticle)), &
ABS(particle_set(iparticle)%v(3) - tmp%vel(3, iparticle)))
END DO
END IF
END IF
END IF
END DO
END SUBROUTINE vv_first
! **************************************************************************************************
!> \brief First half of the velocity-verlet algorithm : update velocity by half
!> step and positions by full step, modified for frozen atoms
!> \param tmp ...
!> \param atomic_kind_set ...
!> \param local_particles ...
!> \param particle_set ...
!> \param core_particle_set ...
!> \param shell_particle_set ...
!> \param nparticle_kind ...
!> \param shell_adiabatic ...
!> \param dt ...
!> \param lfixd_list ...
!> \param u ...
!> \par History
!> none
!> \author Paul Clabaut (Jully 2021)
! **************************************************************************************************
SUBROUTINE vv_first_frozen(tmp, atomic_kind_set, local_particles, particle_set, &
core_particle_set, shell_particle_set, nparticle_kind, &
shell_adiabatic, dt, lfixd_list, u)
TYPE(tmp_variables_type), POINTER :: tmp
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(distribution_1d_type), POINTER :: local_particles
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set, core_particle_set, &
shell_particle_set
INTEGER, INTENT(IN) :: nparticle_kind
LOGICAL, INTENT(IN) :: shell_adiabatic
REAL(KIND=dp) :: dt
TYPE(local_fixd_constraint_type), OPTIONAL, &
POINTER :: lfixd_list(:)
REAL(KIND=dp), DIMENSION(3, 3), OPTIONAL :: u
INTEGER :: iparticle, iparticle_kind, &
iparticle_local, nparticle_local, &
shell_index
LOGICAL :: is_shell
REAL(KIND=dp) :: dm, dmc, dms, dsc(3), dvsc(3), &
fac_massc, fac_masss, mass
TYPE(atomic_kind_type), POINTER :: atomic_kind
TYPE(shell_kind_type), POINTER :: shell
NULLIFY (atomic_kind, shell)
tmp%max_vel = 0.0_dp
tmp%max_vel_sc = 0.0_dp
tmp%max_dvel = 0.0_dp
tmp%max_dvel_sc = 0.0_dp
tmp%max_dr = 0.0_dp
tmp%max_dsc = 0.0_dp
dsc = 0.0_dp
dvsc = 0.0_dp
! *** Velocity Verlet (first part) ***
DO iparticle_kind = 1, nparticle_kind
atomic_kind => atomic_kind_set(iparticle_kind)
CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass, shell_active=is_shell)
IF (mass /= 0.0_dp) THEN
dm = 0.5_dp*dt/mass
IF (is_shell .AND. shell_adiabatic) THEN
CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
dms = 0.5_dp*dt/shell%mass_shell
dmc = 0.5_dp*dt/shell%mass_core
fac_masss = shell%mass_shell/mass
fac_massc = shell%mass_core/mass
nparticle_local = local_particles%n_el(iparticle_kind)
IF (PRESENT(u)) THEN
DO iparticle_local = 1, nparticle_local
iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
shell_index = particle_set(iparticle)%shell_index
! Transform positions and velocities and forces of the shell
CALL transform_first(shell_particle_set, tmp%shell_pos, tmp%shell_vel, &
shell_index, u, dms, dt, tmp%poly_v, tmp%poly_r, tmp%scale_v, tmp%scale_r)
! Transform positions and velocities and forces of the core
CALL transform_first(core_particle_set, tmp%core_pos, tmp%core_vel, &
shell_index, u, dmc, dt, tmp%poly_v, tmp%poly_r, tmp%scale_v, tmp%scale_r)
! Derive velocities and positions of the COM
tmp%vel(:, iparticle) = fac_masss*tmp%shell_vel(:, shell_index) + &
fac_massc*tmp%core_vel(:, shell_index)
tmp%pos(:, iparticle) = fac_masss*tmp%shell_pos(:, shell_index) + &
fac_massc*tmp%core_pos(:, shell_index)
tmp%max_vel = MAX(tmp%max_vel, ABS(tmp%vel(1, iparticle)), &
ABS(tmp%vel(2, iparticle)), ABS(tmp%vel(3, iparticle)))
tmp%max_vel_sc = MAX(tmp%max_vel_sc, &
ABS(tmp%shell_vel(1, shell_index) - tmp%core_vel(1, shell_index)), &
ABS(tmp%shell_vel(2, shell_index) - tmp%core_vel(2, shell_index)), &
ABS(tmp%shell_vel(3, shell_index) - tmp%core_vel(3, shell_index)))
tmp%max_dr = MAX(tmp%max_dr, &
ABS(particle_set(iparticle)%r(1) - tmp%pos(1, iparticle)), &
ABS(particle_set(iparticle)%r(2) - tmp%pos(2, iparticle)), &
ABS(particle_set(iparticle)%r(3) - tmp%pos(3, iparticle)))
tmp%max_dvel = MAX(tmp%max_dvel, &
ABS(particle_set(iparticle)%v(1) - tmp%vel(1, iparticle)), &
ABS(particle_set(iparticle)%v(2) - tmp%vel(2, iparticle)), &
ABS(particle_set(iparticle)%v(3) - tmp%vel(3, iparticle)))
dsc(:) = tmp%shell_pos(:, shell_index) - tmp%core_pos(:, shell_index) - &
shell_particle_set(shell_index)%r(:) + core_particle_set(shell_index)%r(:)
tmp%max_dsc = MAX(tmp%max_dsc, ABS(dsc(1)), ABS(dsc(2)), ABS(dsc(3)))
dvsc(:) = tmp%shell_vel(:, shell_index) - tmp%core_vel(:, shell_index) - &
shell_particle_set(shell_index)%v(:) + core_particle_set(shell_index)%v(:)
tmp%max_dvel_sc = MAX(tmp%max_dvel_sc, ABS(dvsc(1)), ABS(dvsc(2)), ABS(dvsc(3)))
END DO ! iparticle_local
ELSE
DO iparticle_local = 1, nparticle_local
iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
shell_index = particle_set(iparticle)%shell_index
tmp%shell_vel(:, shell_index) = &
shell_particle_set(shell_index)%v(:)*tmp%scale_v(:)*tmp%scale_v(:) + &
tmp%scale_v(:)*tmp%poly_v(:)*dms*shell_particle_set(shell_index)%f(:)
tmp%shell_pos(:, shell_index) = &
shell_particle_set(shell_index)%r(:)*tmp%scale_r(:)*tmp%scale_r(:) + &
tmp%scale_r(:)*tmp%poly_r(:)*dt*tmp%shell_vel(:, shell_index)
tmp%core_vel(:, shell_index) = &
core_particle_set(shell_index)%v(:)*tmp%scale_v(:)*tmp%scale_v(:) + &
tmp%scale_v(:)*tmp%poly_v(:)*dmc*core_particle_set(shell_index)%f(:)
tmp%core_pos(:, shell_index) = &
core_particle_set(shell_index)%r(:)*tmp%scale_r(:)*tmp%scale_r(:) + &
tmp%scale_r(:)*tmp%poly_r(:)*dt*tmp%core_vel(:, shell_index)
tmp%vel(:, iparticle) = fac_masss*tmp%shell_vel(:, shell_index) + &
fac_massc*tmp%core_vel(:, shell_index)
tmp%pos(:, iparticle) = fac_masss*tmp%shell_pos(:, shell_index) + &
fac_massc*tmp%core_pos(:, shell_index)
tmp%max_vel = MAX(tmp%max_vel, &
ABS(tmp%vel(1, iparticle)), ABS(tmp%vel(2, iparticle)), ABS(tmp%vel(3, iparticle)))
tmp%max_vel_sc = MAX(tmp%max_vel_sc, &
ABS(tmp%shell_vel(1, shell_index) - tmp%core_vel(1, shell_index)), &
ABS(tmp%shell_vel(2, shell_index) - tmp%core_vel(2, shell_index)), &
ABS(tmp%shell_vel(3, shell_index) - tmp%core_vel(3, shell_index)))
tmp%max_dr = MAX(tmp%max_dr, &
ABS(particle_set(iparticle)%r(1) - tmp%pos(1, iparticle)), &
ABS(particle_set(iparticle)%r(2) - tmp%pos(2, iparticle)), &
ABS(particle_set(iparticle)%r(3) - tmp%pos(3, iparticle)))
tmp%max_dvel = MAX(tmp%max_dvel, &
ABS(particle_set(iparticle)%v(1) - tmp%vel(1, iparticle)), &
ABS(particle_set(iparticle)%v(2) - tmp%vel(2, iparticle)), &
ABS(particle_set(iparticle)%v(3) - tmp%vel(3, iparticle)))
dsc(:) = tmp%shell_pos(:, shell_index) - tmp%core_pos(:, shell_index) - &
shell_particle_set(shell_index)%r(:) + core_particle_set(shell_index)%r(:)
tmp%max_dsc = MAX(tmp%max_dsc, ABS(dsc(1)), ABS(dsc(2)), ABS(dsc(3)))
dvsc(:) = tmp%shell_vel(:, shell_index) - tmp%core_vel(:, shell_index) - &
shell_particle_set(shell_index)%v(:) + core_particle_set(shell_index)%v(:)
tmp%max_dvel_sc = MAX(tmp%max_dvel_sc, ABS(dvsc(1)), ABS(dvsc(2)), ABS(dvsc(3)))
END DO ! iparticle_local
END IF
ELSE
nparticle_local = local_particles%n_el(iparticle_kind)
IF (PRESENT(u)) THEN
DO iparticle_local = 1, nparticle_local
iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
! Transform positions and velocities and forces
CALL transform_first(particle_set, tmp%pos, tmp%vel, &
iparticle, u, dm, dt, tmp%poly_v, tmp%poly_r, tmp%scale_v, tmp%scale_r)
tmp%max_vel = MAX(tmp%max_vel, &
ABS(tmp%vel(1, iparticle)), ABS(tmp%vel(2, iparticle)), ABS(tmp%vel(3, iparticle)))
tmp%max_dr = MAX(tmp%max_dr, ABS(particle_set(iparticle)%r(1) - tmp%pos(1, iparticle)), &
ABS(particle_set(iparticle)%r(2) - tmp%pos(2, iparticle)), &
ABS(particle_set(iparticle)%r(3) - tmp%pos(3, iparticle)))
tmp%max_dvel = MAX(tmp%max_dvel, &
ABS(particle_set(iparticle)%v(1) - tmp%vel(1, iparticle)), &
ABS(particle_set(iparticle)%v(2) - tmp%vel(2, iparticle)), &
ABS(particle_set(iparticle)%v(3) - tmp%vel(3, iparticle)))
END DO ! iparticle_local
ELSE
DO iparticle_local = 1, nparticle_local
iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
tmp%vel(1, iparticle) = &
particle_set(iparticle)%v(1)*tmp%scale_v(1)*tmp%scale_v(1) + &
tmp%scale_v(1)*tmp%poly_v(1)*dm*particle_set(iparticle)%f(1)
tmp%vel(2, iparticle) = &
particle_set(iparticle)%v(2)*tmp%scale_v(2)*tmp%scale_v(2) + &
tmp%scale_v(2)*tmp%poly_v(2)*dm*particle_set(iparticle)%f(2)
tmp%vel(3, iparticle) = &
particle_set(iparticle)%v(3)*tmp%scale_v(3)*tmp%scale_v(3) + &
tmp%scale_v(3)*tmp%poly_v(3)*dm*particle_set(iparticle)%f(3)
IF (ANY(lfixd_list(:)%ifixd_index == iparticle)) THEN
tmp%pos(1, iparticle) = &
particle_set(iparticle)%r(1)
tmp%pos(2, iparticle) = &
particle_set(iparticle)%r(2)
tmp%pos(3, iparticle) = &
particle_set(iparticle)%r(3)
ELSE
tmp%pos(1, iparticle) = &
particle_set(iparticle)%r(1)*tmp%scale_r(1)*tmp%scale_r(1) + &
tmp%scale_r(1)*tmp%poly_r(1)*dt*tmp%vel(1, iparticle)
tmp%pos(2, iparticle) = &
particle_set(iparticle)%r(2)*tmp%scale_r(2)*tmp%scale_r(2) + &
tmp%scale_r(2)*tmp%poly_r(2)*dt*tmp%vel(2, iparticle)
tmp%pos(3, iparticle) = &
particle_set(iparticle)%r(3)*tmp%scale_r(3)*tmp%scale_r(3) + &
tmp%scale_r(3)*tmp%poly_r(3)*dt*tmp%vel(3, iparticle)
! overwrite positions of frozen particles
IF (PRESENT(lfixd_list)) THEN
IF (ANY(lfixd_list(:)%ifixd_index == iparticle)) THEN
tmp%pos(1, iparticle) = particle_set(iparticle)%r(1)
tmp%pos(2, iparticle) = particle_set(iparticle)%r(2)
tmp%pos(3, iparticle) = particle_set(iparticle)%r(3)
END IF
END IF
tmp%max_vel = MAX(tmp%max_vel, &
@ -1520,7 +1315,7 @@ CONTAINS
END IF
END IF
END DO
END SUBROUTINE vv_first_frozen
END SUBROUTINE vv_first
! **************************************************************************************************
!> \brief Second half of the velocity-verlet algorithm : update velocity by half
@ -1903,7 +1698,7 @@ CONTAINS
core_particle_set, shell_particle_set, nparticle_kind, &
shell_adiabatic, dt)
CASE (npt_i_ensemble, npe_i_ensemble)
CASE (npt_i_ensemble, npt_ia_ensemble, npe_i_ensemble)
arg_r = arg_r*simpar%dt_fact*simpar%dt_fact
tmp%poly_r(1:3) = 1.0_dp + e2*arg_r(1) + e4*arg_r(1)*arg_r(1) + e6*arg_r(1)**3 + e8*arg_r(1)**4
arg_v = arg_v*simpar%dt_fact*simpar%dt_fact
@ -1915,18 +1710,6 @@ CONTAINS
core_particle_set, shell_particle_set, nparticle_kind, &
shell_adiabatic, dt)
CASE (npt_ia_ensemble)
arg_r = arg_r*simpar%dt_fact*simpar%dt_fact
tmp%poly_r(1:3) = 1.0_dp + e2*arg_r(1) + e4*arg_r(1)*arg_r(1) + e6*arg_r(1)**3 + e8*arg_r(1)**4
arg_v = arg_v*simpar%dt_fact*simpar%dt_fact
tmp%poly_v(1:3) = 1.0_dp + e2*arg_v(1) + e4*arg_v(1)*arg_v(1) + e6*arg_v(1)**3 + e8*arg_v(1)**4
tmp%scale_r(1:3) = EXP(0.5_dp*dt*npt(1, 1)%v)
tmp%scale_v(1:3) = EXP(-0.25_dp*dt*npt(1, 1)%v* &
(1.0_dp + 3.0_dp*infree))
CALL vv_first_frozen(tmp, atomic_kind_set, local_particles, particle_set, &
core_particle_set, shell_particle_set, nparticle_kind, &
shell_adiabatic, dt)
CASE (npt_f_ensemble, npe_f_ensemble)
trvg = npt(1, 1)%v + npt(2, 2)%v + npt(3, 3)%v
arg_r(:) = arg_r(:)*simpar%dt_fact*simpar%dt_fact

10
src/motion/velocity_verlet_control.F
View file

@ -19,8 +19,8 @@ MODULE velocity_verlet_control
nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, npt_f_ensemble, npt_i_ensemble, &
npt_ia_ensemble, nve_ensemble, nvt_adiabatic_ensemble, nvt_ensemble, reftraj_ensemble
USE integrator, ONLY: &
isokin, langevin, nph_uniaxial, nph_uniaxial_damped, npt_f, npt_i, npt_ia, nve, nve_respa, &
nvt, nvt_adiabatic, reftraj
isokin, langevin, nph_uniaxial, nph_uniaxial_damped, npt_f, npt_i, nve, nve_respa, nvt, &
nvt_adiabatic, reftraj
USE md_environment_types, ONLY: get_md_env,&
md_environment_type
USE simpar_types, ONLY: simpar_type
@ -79,10 +79,8 @@ CONTAINS
CALL nvt_adiabatic(md_env, globenv)
CASE (isokin_ensemble)
CALL isokin(md_env)
CASE (npt_i_ensemble)
CASE (npt_i_ensemble, npt_ia_ensemble, npe_i_ensemble)
CALL npt_i(md_env, globenv)
CASE (npt_ia_ensemble)
CALL npt_ia(md_env, globenv)
CASE (npt_f_ensemble)
CALL npt_f(md_env, globenv)
CASE (nph_uniaxial_ensemble)
@ -95,8 +93,6 @@ CONTAINS
CALL langevin(md_env)
CASE (npe_f_ensemble)
CALL npt_f(md_env, globenv)
CASE (npe_i_ensemble)
CALL npt_i(md_env, globenv)
END SELECT
CALL timestop(handle)

2
tests/Fist/regtest-gal21/TEST_FILES → tests/Fist/regtest-16/TEST_FILES
View file

@ -3,5 +3,7 @@
# e.g. 0 means do not compare anything, running is enough
# 1 compares the last total energy in the file
# for details see cp2k/tools/do_regtest
test_gal19.inp 11 8e-12 -2.717713309171477
test_gal21.inp 11 8e-12 -2.729695004946153
test_npt_ia.inp 97 1e-6 3.350000
#EOF

0
tests/Fist/regtest-gal19/coord.inc → tests/Fist/regtest-16/coord-gal19.inc
View file

0
tests/Fist/regtest-gal21/coord.inc → tests/Fist/regtest-16/coord-gal21.inc
View file

0
tests/Fist/regtest-npt_ia/coord.inc → tests/Fist/regtest-16/coord-npt_ia.inc
View file

0
tests/Fist/regtest-gal19/gcn.inc → tests/Fist/regtest-16/gcn.inc
View file

4
tests/Fist/regtest-gal19/test_gal.inp → tests/Fist/regtest-16/test_gal19.inp
View file

@ -1,5 +1,5 @@
&GLOBAL
PROJECT md_gal
PROJECT md_gal19
RUN_TYPE energy
&END GLOBAL
@ -55,7 +55,7 @@
PERIODIC XYZ
&END CELL
&COORD
@INCLUDE "coord.inc"
@INCLUDE "coord-gal19.inc"
&END COORD
&TOPOLOGY
&GENERATE

2
tests/Fist/regtest-gal21/test_gal21.inp → tests/Fist/regtest-16/test_gal21.inp
View file

@ -55,7 +55,7 @@
PERIODIC XYZ
&END CELL
&COORD
@INCLUDE "coord.inc"
@INCLUDE "coord-gal21.inc"
&END COORD
&TOPOLOGY
&GENERATE

4
tests/Fist/regtest-npt_ia/test_npt_ia.inp → tests/Fist/regtest-16/test_npt_ia.inp
View file

@ -1,5 +1,5 @@
&GLOBAL
PROJECT test_barostat
PROJECT md_npt_ia
RUN_TYPE md
&PRINT DEBUG
&END
@ -97,7 +97,7 @@
PERIODIC XYZ
&END CELL
&COORD
@INCLUDE "coord.inc"
@INCLUDE "coord-npt_ia.inc"
&END COORD
&TOPOLOGY
&GENERATE

7
tests/Fist/regtest-gal19/TEST_FILES
View file

@ -1,7 +0,0 @@
# runs are executed in the same order as in this file
# the second field tells which test should be run in order to compare with the last available output
# e.g. 0 means do not compare anything, running is enough
# 1 compares the last total energy in the file
# for details see cp2k/tools/do_regtest
test_gal.inp 11 8e-12 -2.717713309171477
#EOF

870
tests/Fist/regtest-gal21/gcn.inc
View file

@ -1,870 +0,0 @@
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
11.25
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
7.5
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0

7
tests/Fist/regtest-npt_ia/TEST_FILES
View file

@ -1,7 +0,0 @@
# runs are executed in the same order as in this file
# the second field tells which test should be run in order to compare with the last available output
# e.g. 0 means do not compare anything, running is enough
# 1 compares the last total energy in the file
# for details see cp2k/tools/do_regtest
test_npt_ia.inp 97 1e-6 3.350000
#EOF

5
tests/TEST_DIRS
View file

@ -306,12 +306,11 @@ SE/regtest-3-1
QMMM/QS/regtest-lrigpw
Fist/regtest-field
Fist/regtest-gauss
Fist/regtest-gal19
Fist/regtest-gal21
QS/regtest-fftw fftw3
QS/regtest-elpa elpa mpiranks==1||mpiranks%2==0
ATOM/regtest-libxc libxc
Fist/regtest-spgr spglib
Fist/regtest-plumed2 plumed2
Fist/regtest-quip quip
Fist/regtest-npt_ia
Fist/regtest-16