From d40d240e9e932ae2deae73189bc2060db0eb963a Mon Sep 17 00:00:00 2001 From: hforbert <42545632+hforbert@users.noreply.github.com> Date: Wed, 23 Feb 2022 14:59:15 +0100 Subject: [PATCH] add CMD to PINT from richard beckman (#1957) add CMD to PINT from richard beckmann --- src/input_constants.F | 3 +- src/motion/pint_methods.F | 87 +++++++++++++++++++++++++++ src/motion/pint_normalmode.F | 38 ++++++++++-- src/motion/pint_pile.F | 11 +++- src/start/input_cp2k_motion.F | 15 +++-- tests/Pimd/regtest-1/TEST_FILES | 2 + tests/Pimd/regtest-1/h2o_cmd_nose.inp | 59 ++++++++++++++++++ tests/Pimd/regtest-1/h2o_cmd_pile.inp | 60 ++++++++++++++++++ 8 files changed, 264 insertions(+), 11 deletions(-) create mode 100644 tests/Pimd/regtest-1/h2o_cmd_nose.inp create mode 100644 tests/Pimd/regtest-1/h2o_cmd_pile.inp diff --git a/src/input_constants.F b/src/input_constants.F index 7796480ea1..cd9946f3d7 100644 --- a/src/input_constants.F +++ b/src/input_constants.F @@ -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 diff --git a/src/motion/pint_methods.F b/src/motion/pint_methods.F index 297f9cc8a8..16330aba21 100644 --- a/src/motion/pint_methods.F +++ b/src/motion/pint_methods.F @@ -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 diff --git a/src/motion/pint_normalmode.F b/src/motion/pint_normalmode.F index 3b67a75c7c..752567f36a 100644 --- a/src/motion/pint_normalmode.F +++ b/src/motion/pint_normalmode.F @@ -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)) diff --git a/src/motion/pint_pile.F b/src/motion/pint_pile.F index 6163d7d81d..2ea74beebb 100644 --- a/src/motion/pint_pile.F +++ b/src/motion/pint_pile.F @@ -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 diff --git a/src/start/input_cp2k_motion.F b/src/start/input_cp2k_motion.F index 82da8e1768..968763bebb 100644 --- a/src/start/input_cp2k_motion.F +++ b/src/start/input_cp2k_motion.F @@ -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) diff --git a/tests/Pimd/regtest-1/TEST_FILES b/tests/Pimd/regtest-1/TEST_FILES index da26173312..bcc56c9c4a 100644 --- a/tests/Pimd/regtest-1/TEST_FILES +++ b/tests/Pimd/regtest-1/TEST_FILES @@ -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 diff --git a/tests/Pimd/regtest-1/h2o_cmd_nose.inp b/tests/Pimd/regtest-1/h2o_cmd_nose.inp new file mode 100644 index 0000000000..3cdf1abc9c --- /dev/null +++ b/tests/Pimd/regtest-1/h2o_cmd_nose.inp @@ -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 diff --git a/tests/Pimd/regtest-1/h2o_cmd_pile.inp b/tests/Pimd/regtest-1/h2o_cmd_pile.inp new file mode 100644 index 0000000000..3985aa5ccf --- /dev/null +++ b/tests/Pimd/regtest-1/h2o_cmd_pile.inp @@ -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