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

add BCMD propagator to PINT (#5247)

Browse source
This commit is contained in:
hforbert 2026-05-19 09:57:35 +02:00 • committed by GitHub
parent f658ef2d91
commit e068a9b9a7
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
8 changed files with 154 additions and 26 deletions
Download patch file Download diff file Expand all files Collapse all files

9
src/common/bibliography.F
View file

@ -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"), &

3
src/input_constants.F
View file

@ -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

63
src/motion/pint_methods.F
View file

@ -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

8
src/motion/pint_normalmode.F
View file

@ -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

33
src/motion/pint_pile.F
View file

@ -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))

9
src/start/input_cp2k_motion.F
View file

@ -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", &

1
tests/Pimd/regtest-1/TEST_FILES.toml
View file

@ -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

54
tests/Pimd/regtest-1/h2o_bcmd.inp Normal file
View file

@ -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