From e068a9b9a79ff801c506cc5bc38f62940e03fd5f Mon Sep 17 00:00:00 2001 From: hforbert <42545632+hforbert@users.noreply.github.com> Date: Tue, 19 May 2026 09:57:35 +0200 Subject: [PATCH] add BCMD propagator to PINT (#5247) --- src/common/bibliography.F | 9 +++- src/input_constants.F | 3 +- src/motion/pint_methods.F | 63 +++++++++++++++++++++++++--- src/motion/pint_normalmode.F | 8 ++-- src/motion/pint_pile.F | 33 ++++++++++----- src/start/input_cp2k_motion.F | 9 ++-- tests/Pimd/regtest-1/TEST_FILES.toml | 1 + tests/Pimd/regtest-1/h2o_bcmd.inp | 54 ++++++++++++++++++++++++ 8 files changed, 154 insertions(+), 26 deletions(-) create mode 100644 tests/Pimd/regtest-1/h2o_bcmd.inp diff --git a/src/common/bibliography.F b/src/common/bibliography.F index 4591077e3e..762ee7f6e7 100644 --- a/src/common/bibliography.F +++ b/src/common/bibliography.F @@ -95,7 +95,8 @@ MODULE bibliography VazdaCruz2021, Chen2025, Hernandez2025, Marek2025, Hehn2022, Hehn2024, Pasquier2025, & Hanasaki2025, Tan2025, Broyden1965, Johnson1988, Kerker1981, & FuHo1983, MethfesselPaxton1989, Marzari1999, dosSantos2023, Mermin1965, & - KuhneHeskeProdan2020, Schreder2021, Schreder2024_1, Schreder2024_2 + KuhneHeskeProdan2020, Schreder2021, Schreder2024_1, Schreder2024_2, & + Shiga2022 CONTAINS @@ -2048,6 +2049,12 @@ CONTAINS "the restricted subspace approximation", & source="Phys. Chem. Chem. Phys.", volume="23", pages="1835-1848", & year=2021, doi="10.1039/d0cp04726k") + CALL add_reference(key=Shiga2022, & + authors=s2a("M. Shiga"), & + title="Path integral Brownian chain molecular dynamics: A simple approximation "// & + "of quantum vibrational dynamics", & + source="J. Comput. Chem.", volume="43", pages="1864-1879", & + year=2022, doi="10.1002/jcc.26989") CALL add_reference(key=Chen2025, & authors=s2a("Z. Chen", "Y. Yang"), & diff --git a/src/input_constants.F b/src/input_constants.F index 7bf4b498d2..12e590a177 100644 --- a/src/input_constants.F +++ b/src/input_constants.F @@ -310,7 +310,8 @@ MODULE input_constants ! pint propagator mode INTEGER, PARAMETER, PUBLIC :: propagator_pimd = 1, & propagator_rpmd = 2, & - propagator_cmd = 3 + propagator_cmd = 3, & + propagator_bcmd = 4 ! 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 ccb4dfdc14..1f59bcbb61 100644 --- a/src/motion/pint_methods.F +++ b/src/motion/pint_methods.F @@ -23,6 +23,7 @@ MODULE pint_methods USE bibliography, ONLY: Brieuc2016,& Ceriotti2010,& Ceriotti2012,& + Shiga2022,& cite_reference USE cell_types, ONLY: cell_type USE constraint, ONLY: rattle_control,& @@ -67,15 +68,16 @@ MODULE pint_methods USE helium_types, ONLY: helium_solvent_p_type USE input_constants, ONLY: integrate_exact,& integrate_numeric,& + propagator_bcmd,& propagator_cmd,& propagator_rpmd,& transformation_normal,& transformation_stage USE input_cp2k_restarts, ONLY: write_restart USE input_section_types, ONLY: & - section_type, section_vals_get, section_vals_get_subs_vals, section_vals_release, & - section_vals_retain, section_vals_type, section_vals_val_get, section_vals_val_set, & - section_vals_val_unset + section_type, section_vals_add_values, section_vals_get, section_vals_get_subs_vals, & + section_vals_release, section_vals_retain, section_vals_type, section_vals_val_get, & + section_vals_val_set, section_vals_val_unset USE kinds, ONLY: default_path_length,& default_string_length,& dp @@ -294,7 +296,12 @@ CONTAINS IF ((pint_env%propagator%prop_kind == propagator_cmd) .AND. & (pint_env%transform /= transformation_normal)) THEN - CPABORT("CMD Propagator without normal modes not implemented!") + CPABORT("CMD propagator without normal modes not implemented!") + END IF + + IF ((pint_env%propagator%prop_kind == propagator_bcmd) .AND. & + (pint_env%transform /= transformation_normal)) THEN + CPABORT("BCMD 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) @@ -305,7 +312,10 @@ CONTAINS CALL section_vals_get(nose_section, explicit=explicit) IF (explicit) THEN IF (pint_env%propagator%prop_kind == propagator_rpmd) THEN - CPABORT("RPMD Propagator with Nose-thermostat not implemented!") + CPABORT("RPMD propagator with Nose-thermostat not implemented!") + END IF + IF (pint_env%propagator%prop_kind == propagator_bcmd) THEN + CPABORT("BCMD propagator with Nose-thermostat not implemented!") END IF CALL section_vals_val_get(nose_section, "nnos", i_val=pint_env%nnos) IF (pint_env%nnos > 0) THEN @@ -395,6 +405,16 @@ CONTAINS ELSE transform_section => section_vals_get_subs_vals(input, & "MOTION%PINT%NORMALMODE") + IF (pint_env%propagator%prop_kind == propagator_bcmd) THEN + CALL section_vals_val_get(transform_section, "GAMMA", explicit=explicit) + IF (explicit) THEN + CALL section_vals_val_set(transform_section, "GAMMA", & + r_val=SQRT(2.0_dp/(pint_env%p*pint_env%dt*pint_env%kT))) + ELSE + CALL section_vals_val_set(transform_section, "MODEFACTOR", & + r_val=0.5_dp*pint_env%p*pint_env%dt*pint_env%kT) + END IF + END IF ALLOCATE (pint_env%normalmode_env) CALL normalmode_env_create(pint_env%normalmode_env, & transform_section, p=pint_env%p, kT=pint_env%kT, propagator=pint_env%propagator%prop_kind) @@ -517,7 +537,7 @@ CONTAINS piglet_section, & dt=pint_env%dt, para_env=para_env) ELSE - CPABORT("PILE thermostat can't be used with another thermostat.") + CPABORT("PIGLET thermostat can't be used with another thermostat.") END IF END IF @@ -541,6 +561,31 @@ CONTAINS END IF END IF + ! force PILE for BCMD + IF (pint_env%propagator%prop_kind == propagator_bcmd) THEN + IF (pint_env%pimd_thermostat == thermostat_none) THEN + !crate section + pile_section => section_vals_get_subs_vals(input, "MOTION%PINT%PILE") + CALL section_vals_get(pile_section, explicit=explicit) + IF (.NOT. explicit) THEN + CALL section_vals_add_values(pile_section) + END IF + CALL section_vals_val_get(pint_env%input, & + "MOTION%PINT%INIT%THERMOSTAT_SEED", & + i_val=pint_env%thermostat_rng_seed) + pint_env%pimd_thermostat = thermostat_pile + ALLOCATE (pint_env%pile_therm) + CALL pint_pile_init(pile_therm=pint_env%pile_therm, & + pint_env=pint_env, & + normalmode_env=pint_env%normalmode_env, & + section=pile_section) + END IF + IF (pint_env%pimd_thermostat /= thermostat_pile) THEN + CPABORT("PILE/no thermostat currently needed for BCMD") + END IF + CALL cite_reference(Shiga2022) + END IF + !Initialize integrator scheme CALL section_vals_val_get(pint_section, "HARM_INT", i_val=pint_env%harm_integrator) IF (pint_env%harm_integrator == integrate_exact) THEN @@ -575,6 +620,12 @@ CONTAINS CPWARN(msg) pint_env%harm_integrator = integrate_exact END IF + IF (pint_env%propagator%prop_kind == propagator_bcmd) THEN + WRITE (UNIT=msg, FMT=*) "PINT WARNING| BCMD needs the exact harmonic "// & + "integrator. Switching to exact harmonic integrator." + CPWARN(msg) + pint_env%harm_integrator = integrate_exact + END IF END IF IF (pint_env%harm_integrator == integrate_exact) THEN diff --git a/src/motion/pint_normalmode.F b/src/motion/pint_normalmode.F index 509d6de2a8..d2a579826b 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_cmd,& + USE input_constants, ONLY: propagator_bcmd,& + propagator_cmd,& propagator_pimd,& propagator_rpmd USE input_section_types, ONLY: section_vals_type,& @@ -105,9 +106,10 @@ CONTAINS END IF !Use different normal mode transformations depending on the propagator - IF (propagator == propagator_pimd .OR. propagator == propagator_cmd) THEN + IF (propagator == propagator_pimd .OR. propagator == propagator_cmd .OR. & + propagator == propagator_bcmd) THEN - IF (propagator == propagator_pimd) THEN + IF (propagator == propagator_pimd .OR. propagator == propagator_bcmd) 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 diff --git a/src/motion/pint_pile.F b/src/motion/pint_pile.F index 7b78b648e7..9892158811 100644 --- a/src/motion/pint_pile.F +++ b/src/motion/pint_pile.F @@ -13,7 +13,8 @@ !> 10.2014 created [Felix Uhl] ! ************************************************************************************************** MODULE pint_pile - USE input_constants, ONLY: propagator_cmd + USE input_constants, ONLY: propagator_bcmd,& + propagator_cmd USE input_section_types, ONLY: section_vals_get,& section_vals_get_subs_vals,& section_vals_type,& @@ -71,7 +72,7 @@ CONTAINS 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 + IF (i_propagator == propagator_cmd .OR. i_propagator == propagator_bcmd) THEN pile_therm%tau = 0.0_dp END IF @@ -89,15 +90,25 @@ CONTAINS ELSE pile_therm%g_fric(1) = 1.0_dp/pile_therm%tau END IF - DO i = 2, p - pile_therm%g_fric(i) = 2.0_dp*pile_therm%lamb*SQRT(normalmode_env%lambda(i)) - END DO - DO i = 1, p - ex = -dti2*pile_therm%g_fric(i) - pile_therm%c1(i) = EXP(ex) - ex = pile_therm%c1(i)*pile_therm%c1(i) - pile_therm%c2(i) = SQRT(1.0_dp - ex) - END DO + IF (i_propagator == propagator_bcmd) THEN + pile_therm%c1(1) = EXP(-dti2*pile_therm%g_fric(1)) + ex = pile_therm%c1(1)*pile_therm%c1(1) + pile_therm%c2(1) = SQRT(1.0_dp - ex) + DO i = 2, p + pile_therm%c1(i) = 0.0_dp + pile_therm%c2(i) = 1.0_dp + END DO + ELSE + DO i = 2, p + pile_therm%g_fric(i) = 2.0_dp*pile_therm%lamb*SQRT(normalmode_env%lambda(i)) + END DO + DO i = 1, p + ex = -dti2*pile_therm%g_fric(i) + pile_therm%c1(i) = EXP(ex) + ex = pile_therm%c1(i)*pile_therm%c1(i) + pile_therm%c2(i) = SQRT(1.0_dp - ex) + END DO + END IF DO j = 1, pint_env%ndim DO i = 1, pint_env%p pile_therm%massfact(i, j) = SQRT(pint_env%kT/pint_env%mass_fict(i, j)) diff --git a/src/start/input_cp2k_motion.F b/src/start/input_cp2k_motion.F index 8c41eec28f..13fea8a074 100644 --- a/src/start/input_cp2k_motion.F +++ b/src/start/input_cp2k_motion.F @@ -39,7 +39,8 @@ MODULE input_cp2k_motion helium_sampling_ceperley, helium_sampling_worm, helium_solute_intpot_mwater, helium_solute_intpot_nnp, & 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, propagator_cmd, transformation_normal, transformation_stage + perm_plain, propagator_pimd, propagator_rpmd, propagator_cmd, propagator_bcmd, 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 @@ -1752,10 +1753,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|CMD)", & + usage="PROPAGATOR (PIMD|RPMD|CMD|BCMD)", & default_i_val=propagator_pimd, & - enum_c_vals=s2a("PIMD", "RPMD", "CMD"), & - enum_i_vals=[propagator_pimd, propagator_rpmd, propagator_cmd]) + enum_c_vals=s2a("PIMD", "RPMD", "CMD", "BCMD"), & + enum_i_vals=[propagator_pimd, propagator_rpmd, propagator_cmd, propagator_bcmd]) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) CALL keyword_create(keyword, __LOCATION__, name="FIX_CENTROID_POS", & diff --git a/tests/Pimd/regtest-1/TEST_FILES.toml b/tests/Pimd/regtest-1/TEST_FILES.toml index ee66807360..d68ace4cee 100644 --- a/tests/Pimd/regtest-1/TEST_FILES.toml +++ b/tests/Pimd/regtest-1/TEST_FILES.toml @@ -18,4 +18,5 @@ "water-in-helium-lastforce.inp" = [{matcher="M009", tol=1e-13, ref=1.0708166830830485E-003}] "h2o_cmd_nose.inp" = [{matcher="M009", tol=1e-13, ref=2.5049926484361770E-002}] "h2o_cmd_pile.inp" = [{matcher="M009", tol=1e-13, ref=6.7246877655376420E-003}] +"h2o_bcmd.inp" = [{matcher="M009", tol=3e-14, ref=1.1854234994934626E-002}] #EOF diff --git a/tests/Pimd/regtest-1/h2o_bcmd.inp b/tests/Pimd/regtest-1/h2o_bcmd.inp new file mode 100644 index 0000000000..f782c0e6a5 --- /dev/null +++ b/tests/Pimd/regtest-1/h2o_bcmd.inp @@ -0,0 +1,54 @@ +&GLOBAL + PRINT_LEVEL LOW + PROJECT_NAME h2o_bcmd + RUN_TYPE PINT +&END GLOBAL + +&MOTION + &PINT + DT 0.5 + HARM_INT EXACT + NRESPA 1 + NUM_STEPS 100 + P 4 + PROC_PER_REPLICA 1 + PROPAGATOR BCMD + TEMP 200.0 + TRANSFORMATION NORMAL + &END PINT +&END MOTION + +&FORCE_EVAL + METHOD FIST + &MM + &FORCEFIELD + PARMTYPE CHM + PARM_FILE_NAME ../../Fist/sample_pot/water.pot + &CHARGE + ATOM OT + CHARGE -0.8476 + &END CHARGE + &CHARGE + ATOM HT + CHARGE 0.4238 + &END CHARGE + &END FORCEFIELD + &POISSON + &EWALD + ALPHA 0.44 + EWALD_TYPE SPME + 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_FORMAT PDB + COORD_FILE_NAME ../../Fist/sample_pdb/water_1.pdb + &END TOPOLOGY + &END SUBSYS +&END FORCE_EVAL