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

add CMD to PINT from richard beckman (#1957)

Browse source 
add CMD to PINT from richard beckmann
This commit is contained in:
hforbert 2022-02-23 14:59:15 +01:00 • committed by GitHub
parent f2b1148c02
commit d40d240e9e
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
8 changed files with 264 additions and 11 deletions
Download patch file Download diff file Expand all files Collapse all files

3
src/input_constants.F
View file

@ -296,7 +296,8 @@ MODULE input_constants
! pint propagator mode
INTEGER, PARAMETER, PUBLIC :: propagator_pimd = 1, &
propagator_rpmd = 2
propagator_rpmd = 2, &
propagator_cmd = 3
! pint propagator mode
INTEGER, PARAMETER, PUBLIC :: integrate_numeric = 1, &
integrate_exact = 2

87
src/motion/pint_methods.F
View file

@ -68,6 +68,7 @@ MODULE pint_methods
USE helium_types, ONLY: helium_solvent_p_type
USE input_constants, ONLY: integrate_exact,&
integrate_numeric,&
propagator_cmd,&
propagator_rpmd,&
transformation_normal,&
transformation_stage
@ -299,6 +300,11 @@ CONTAINS
CALL section_vals_val_get(pint_section, "transformation", &
i_val=pint_env%transform)
IF ((pint_env%propagator%prop_kind == propagator_cmd) .AND. &
(pint_env%transform /= transformation_normal)) THEN
CPABORT("CMD Propagator without normal modes not implemented!")
END IF
NULLIFY (pint_env%tx, pint_env%tv, pint_env%tv_t, pint_env%tv_old, pint_env%tv_new, pint_env%tf)
pint_env%nnos = 0
@ -1619,6 +1625,9 @@ CONTAINS
END DO
END DO
END DO
IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
pint_env%tv(:, 1, :) = 0.0_dp
END IF
NULLIFY (input_section)
input_section => section_vals_get_subs_vals(pint_env%input, &
@ -1644,6 +1653,9 @@ CONTAINS
END DO
END IF
END IF
IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
pint_env%tx(:, 1, :) = 0.0_dp
END IF
NULLIFY (input_section)
input_section => section_vals_get_subs_vals(pint_env%input, &
@ -1668,6 +1680,9 @@ CONTAINS
END DO
END DO
END IF
IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
pint_env%tv(:, 1, :) = 0.0_dp
END IF
END IF
ELSEIF (pint_env%pimd_thermostat == thermostat_gle) THEN
@ -2072,12 +2087,25 @@ CONTAINS
END DO
END IF
! Exempt centroid from thermostat for CMD
IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
pint_env%tx(:, 1, :) = 0.0_dp
pint_env%tv(:, 1, :) = 0.0_dp
pint_env%tf(:, 1, :) = 0.0_dp
END IF
DO i = pint_env%first_propagated_mode, pint_env%p
pint_env%ux(i, :) = pint_env%ux(i, :) - &
dti22*pint_env%uv(i, :)*pint_env%tv(1, i, :)
END DO
pint_env%tx = pint_env%tx + dti*pint_env%tv + dti22*pint_env%tf
IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
pint_env%tx(:, 1, :) = 0.0_dp
pint_env%tv(:, 1, :) = 0.0_dp
pint_env%tf(:, 1, :) = 0.0_dp
END IF
END IF
!Integrate position in harmonic springs (uf_h) and physical potential
!(uf)
@ -2091,6 +2119,13 @@ CONTAINS
! apply thermostats to velocities
SELECT CASE (pint_env%pimd_thermostat)
CASE (thermostat_nose)
IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
pint_env%tx(:, 1, :) = 0.0_dp
pint_env%tv(:, 1, :) = 0.0_dp
pint_env%tf(:, 1, :) = 0.0_dp
END IF
pint_env%uv_t = pint_env%uv - dti2* &
pint_env%uv*pint_env%tv(1, :, :)
tmp => pint_env%tv_t
@ -2114,6 +2149,13 @@ CONTAINS
END DO
END IF
! Exempt centroid from thermostat for CMD
IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
pint_env%tx(:, 1, :) = 0.0_dp
pint_env%tv(:, 1, :) = 0.0_dp
pint_env%tf(:, 1, :) = 0.0_dp
END IF
!Integrate harmonic velocities and physical velocities
pint_env%uv_t = pint_env%uv_t + dti2*(pint_env%uf_h + pint_env%uf)
@ -2160,6 +2202,12 @@ CONTAINS
END DO
END IF
! Exempt centroid from thermostat for CMD
IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
pint_env%tx(:, 1, :) = 0.0_dp
pint_env%tv(:, 1, :) = 0.0_dp
pint_env%tf(:, 1, :) = 0.0_dp
END IF
CALL pint_calc_uf_h(pint_env=pint_env, e_h=e_h)
pint_env%uv_t = pint_env%uv_t + dti2*(pint_env%uf_h + pint_env%uf)
@ -2196,6 +2244,12 @@ CONTAINS
! Apply second half of thermostats
SELECT CASE (pint_env%pimd_thermostat)
CASE (thermostat_nose)
! Exempt centroid from thermostat for CMD
IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
pint_env%tx(:, 1, :) = 0.0_dp
pint_env%tv(:, 1, :) = 0.0_dp
pint_env%tf(:, 1, :) = 0.0_dp
END IF
DO i = 1, 6
tol = 0._dp
pint_env%uv_new = pint_env%uv_t/(1.+dti2*pint_env%tv(1, :, :))
@ -2217,6 +2271,13 @@ CONTAINS
END DO
END IF
! Exempt centroid from thermostat for CMD
IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
pint_env%tx(:, 1, :) = 0.0_dp
pint_env%tv(:, 1, :) = 0.0_dp
pint_env%tf(:, 1, :) = 0.0_dp
END IF
DO idim = 1, pint_env%ndim
DO ib = 1, pint_env%p
DO inos = 1, pint_env%nnos - 1
@ -2240,6 +2301,13 @@ CONTAINS
END DO
END IF
! Exempt centroid from thermostat for CMD
IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
pint_env%tx(:, 1, :) = 0.0_dp
pint_env%tv(:, 1, :) = 0.0_dp
pint_env%tf(:, 1, :) = 0.0_dp
END IF
pint_env%tv_new(pint_env%nnos, ib, idim) = &
pint_env%tv_t(pint_env%nnos, ib, idim) + &
dti2*pint_env%tf(pint_env%nnos, ib, idim)
@ -2256,6 +2324,12 @@ CONTAINS
END DO
END DO
END IF
! Exempt centroid from thermostat for CMD
IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
pint_env%tx(:, 1, :) = 0.0_dp
pint_env%tv(:, 1, :) = 0.0_dp
pint_env%tf(:, 1, :) = 0.0_dp
END IF
END DO
END DO
@ -2263,6 +2337,12 @@ CONTAINS
pint_env%uv = pint_env%uv_new
pint_env%tv = pint_env%tv_new
IF (tol <= pint_env%v_tol) EXIT
! Exempt centroid from thermostat for CMD
IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
pint_env%tx(:, 1, :) = 0.0_dp
pint_env%tv(:, 1, :) = 0.0_dp
pint_env%tf(:, 1, :) = 0.0_dp
END IF
END DO
! Apply centroid constraints (RATTLE)
@ -2309,6 +2389,13 @@ CONTAINS
- pint_env%tv(inos, :, :)*pint_env%tv(inos + 1, :, :)
END DO
! Exempt centroid from thermostat for CMD
IF (pint_env%propagator%prop_kind == propagator_cmd) THEN
pint_env%tx(:, 1, :) = 0.0_dp
pint_env%tv(:, 1, :) = 0.0_dp
pint_env%tf(:, 1, :) = 0.0_dp
END IF
CASE (thermostat_gle)
CALL pint_gle_step(pint_env)
pint_env%uv = pint_env%uv_t

38
src/motion/pint_normalmode.F
View file

@ -16,7 +16,8 @@
!> [Felix Uhl
! **************************************************************************************************
MODULE pint_normalmode
USE input_constants, ONLY: propagator_pimd,&
USE input_constants, ONLY: propagator_cmd,&
propagator_pimd,&
propagator_rpmd
USE input_section_types, ONLY: section_vals_type,&
section_vals_val_get
@ -61,7 +62,9 @@ CONTAINS
INTEGER, INTENT(in) :: propagator
INTEGER :: i, j, k, li
REAL(kind=dp) :: invsqrtp, pip, sqrt2p, twopip
LOGICAL :: explicit_gamma, explicit_modefactor
REAL(kind=dp) :: gamma_parameter, invsqrtp, pip, sqrt2p, &
twopip
CPASSERT(.NOT. ASSOCIATED(normalmode_env))
ALLOCATE (normalmode_env)
@ -79,7 +82,27 @@ CONTAINS
CALL section_vals_val_get(normalmode_section, "Q_BEAD", &
r_val=normalmode_env%Q_bead)
CALL section_vals_val_get(normalmode_section, "MODEFACTOR", &
explicit=explicit_modefactor, &
r_val=normalmode_env%modefactor)
CALL section_vals_val_get(normalmode_section, "GAMMA", &
r_val=gamma_parameter, &
explicit=explicit_gamma)
IF (explicit_modefactor .AND. explicit_gamma) THEN
CPABORT("Both GAMMA and MODEFACTOR have been declared. Please use only one.")
END IF
IF (explicit_gamma) THEN
normalmode_env%modefactor = 1.0_dp/gamma_parameter**2
END IF
IF (propagator == propagator_cmd) THEN
IF (.NOT. explicit_gamma) THEN
CPABORT("GAMMA needs to be specified with CMD PROPAGATOR")
END IF
IF (gamma_parameter <= 1.0_dp) THEN
CPWARN("GAMMA should be larger than 1.0 for CMD PROPAGATOR")
END IF
END IF
IF (normalmode_env%Q_centroid < 0.0_dp) THEN
normalmode_env%Q_centroid = -normalmode_env%Q_centroid/(kT*p)
@ -89,8 +112,14 @@ CONTAINS
END IF
!Use different normal mode transformations depending on the propagator
IF (propagator == propagator_pimd) THEN
normalmode_env%harm = p*kT*kT/normalmode_env%modefactor
IF (propagator == propagator_pimd .OR. propagator == propagator_cmd) THEN
IF (propagator == propagator_pimd) THEN
normalmode_env%harm = p*kT*kT/normalmode_env%modefactor
ELSE IF (propagator == propagator_cmd) THEN
normalmode_env%harm = p*kT*kT*gamma_parameter*gamma_parameter
normalmode_env%modefactor = 1.0_dp/(gamma_parameter*gamma_parameter)
END IF
! set up the transformation matrices
DO i = 1, p
@ -118,6 +147,7 @@ CONTAINS
END DO
END DO
normalmode_env%lambda(:) = normalmode_env%harm
ELSE IF (propagator == propagator_rpmd) THEN
normalmode_env%harm = kT/normalmode_env%modefactor
sqrt2p = SQRT(2.0_dp/REAL(p, dp))

11
src/motion/pint_pile.F
View file

@ -13,6 +13,7 @@
!> 10.2014 created [Felix Uhl]
! **************************************************************************************************
MODULE pint_pile
USE input_constants, ONLY: propagator_cmd
USE input_section_types, ONLY: section_vals_get,&
section_vals_get_subs_vals,&
section_vals_type,&
@ -58,11 +59,11 @@ CONTAINS
TYPE(section_vals_type), POINTER :: section
CHARACTER(LEN=rng_record_length) :: rng_record
INTEGER :: i, j, p
INTEGER :: i, i_propagator, j, p
LOGICAL :: explicit
REAL(KIND=dp) :: dti2, ex
REAL(KIND=dp), DIMENSION(3, 2) :: initial_seed
TYPE(section_vals_type), POINTER :: rng_section
TYPE(section_vals_type), POINTER :: pint_section, rng_section
pint_env%e_pile = 0.0_dp
ALLOCATE (pile_therm)
@ -72,6 +73,12 @@ CONTAINS
CALL section_vals_val_get(section, "TAU", r_val=pile_therm%tau)
CALL section_vals_val_get(section, "LAMBDA", r_val=pile_therm%lamb)
CALL section_vals_val_get(section, "THERMOSTAT_ENERGY", r_val=pile_therm%thermostat_energy)
pint_section => section_vals_get_subs_vals(pint_env%input, "MOTION%PINT")
CALL section_vals_val_get(pint_section, "PROPAGATOR", i_val=i_propagator)
IF (i_propagator == propagator_cmd) THEN
pile_therm%tau = 0.0_dp
END IF
p = pint_env%p
dti2 = 0.5_dp*pint_env%dt

15
src/start/input_cp2k_motion.F
View file

@ -39,7 +39,7 @@ MODULE input_cp2k_motion
helium_sampling_ceperley, helium_sampling_worm, helium_solute_intpot_mwater, &
helium_solute_intpot_none, integrate_exact, integrate_numeric, ls_2pnt, ls_3pnt, ls_fit, &
ls_gold, ls_none, matrix_init_cholesky, matrix_init_diagonal, numerical, perm_cycle, &
perm_plain, propagator_pimd, propagator_rpmd, transformation_normal, transformation_stage
perm_plain, propagator_pimd, propagator_rpmd, propagator_cmd, transformation_normal, transformation_stage
USE input_cp2k_constraints, ONLY: create_constraint_section
USE input_cp2k_free_energy, ONLY: create_fe_section
USE input_cp2k_md, ONLY: create_md_section
@ -1721,10 +1721,10 @@ CONTAINS
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="propagator", &
description="Specifies the real time propagator to use", &
usage="PROPAGATOR (PIMD|RPMD)", &
usage="PROPAGATOR (PIMD|RPMD|CMD)", &
default_i_val=propagator_pimd, &
enum_c_vals=s2a("PIMD", "RPMD"), &
enum_i_vals=(/propagator_pimd, propagator_rpmd/))
enum_c_vals=s2a("PIMD", "RPMD", "CMD"), &
enum_i_vals=(/propagator_pimd, propagator_rpmd, propagator_cmd/))
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="FIX_CENTROID_POS", &
@ -1755,6 +1755,13 @@ CONTAINS
repeats=.FALSE., default_r_val=1.0_dp)
CALL section_add_keyword(subsection, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="GAMMA", &
description="mass scale factor for non-centroid degrees of freedom, &
& naming convention according to Witt, 2008, doi.org/10.1063/1.3125009", &
repeats=.FALSE., default_r_val=8.0_dp)
CALL section_add_keyword(subsection, keyword)
CALL keyword_release(keyword)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)

2
tests/Pimd/regtest-1/TEST_FILES
View file

@ -16,4 +16,6 @@ he32_only_restart.inp 40 2e-13 -1.437
water-in-helium.inp 9 7e-14 1.0717546824070511E-003
water-in-helium-restart.inp 9 1e-13 1.0716120146125518E-003
water-in-helium-lastforce.inp 9 1e-13 1.0708166830830485E-003
h2o_cmd_nose.inp 9 1e-13 2.5049926484361770E-002
h2o_cmd_pile.inp 9 1e-13 6.7246877655376420E-003
#EOF

59
tests/Pimd/regtest-1/h2o_cmd_nose.inp Normal file
View file

@ -0,0 +1,59 @@
&GLOBAL
PROJECT_NAME h2o_cmd_nose
RUN_TYPE PINT
PRINT_LEVEL LOW
&END GLOBAL
&MOTION
&PINT
P 10
PROC_PER_REPLICA 1
NUM_STEPS 40
DT 0.5
NRESPA 2
TEMP 50.0
TRANSFORMATION NORMAL
PROPAGATOR CMD
&NORMALMODE
GAMMA 8
&END NORMALMODE
&NOSE
NNOS 3
&END NOSE
&END PINT
&END MOTION
&FORCE_EVAL
METHOD FIST
&MM
&FORCEFIELD
PARM_FILE_NAME ../../Fist/sample_pot/water.pot
PARMTYPE CHM
&CHARGE
ATOM OT
CHARGE -0.8476
&END CHARGE
&CHARGE
ATOM HT
CHARGE 0.4238
&END CHARGE
&END FORCEFIELD
&POISSON
&EWALD
EWALD_TYPE SPME
ALPHA 0.44
GMAX 24
O_SPLINE 6
&END EWALD
&END POISSON
&END MM
&SUBSYS
&CELL
ABC 24.955 24.955 24.955
&END CELL
&TOPOLOGY
COORD_FILE_NAME ../../Fist/sample_pdb/water_1.pdb
COORD_FILE_FORMAT PDB
&END TOPOLOGY
&END SUBSYS
&END FORCE_EVAL

60
tests/Pimd/regtest-1/h2o_cmd_pile.inp Normal file
View file

@ -0,0 +1,60 @@
&GLOBAL
PROJECT_NAME h2o_cmd_pile
RUN_TYPE PINT
PRINT_LEVEL LOW
&END GLOBAL
&MOTION
&PINT
P 10
PROC_PER_REPLICA 1
NUM_STEPS 40
DT 0.5
NRESPA 2
TEMP 50.0
TRANSFORMATION NORMAL
PROPAGATOR CMD
&NORMALMODE
GAMMA 8
&END NORMALMODE
&PILE
TAU 1000
LAMBDA 0.5
&END PILE
&END PINT
&END MOTION
&FORCE_EVAL
METHOD FIST
&MM
&FORCEFIELD
PARM_FILE_NAME ../../Fist/sample_pot/water.pot
PARMTYPE CHM
&CHARGE
ATOM OT
CHARGE -0.8476
&END CHARGE
&CHARGE
ATOM HT
CHARGE 0.4238
&END CHARGE
&END FORCEFIELD
&POISSON
&EWALD
EWALD_TYPE SPME
ALPHA 0.44
GMAX 24
O_SPLINE 6
&END EWALD
&END POISSON
&END MM
&SUBSYS
&CELL
ABC 24.955 24.955 24.955
&END CELL
&TOPOLOGY
COORD_FILE_NAME ../../Fist/sample_pdb/water_1.pdb
COORD_FILE_FORMAT PDB
&END TOPOLOGY
&END SUBSYS
&END FORCE_EVAL