Fist: Add switch to turn off real space electrostatics

This commit is contained in:
Da Teng 2023-06-30 08:12:06 -05:00 committed by GitHub
parent a20b57f0a3
commit 72e1b4a66a
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
9 changed files with 176 additions and 16 deletions

View file

@ -75,6 +75,7 @@ MODULE fist_nonbond_env_types
INTEGER :: last_update = -1
INTEGER :: num_update = -1
LOGICAL :: do_nonbonded = .FALSE.
LOGICAL :: do_electrostatics = .FALSE.
LOGICAL :: shift_cutoff = .FALSE.
CHARACTER(len=default_string_length) :: unit_type = ""
REAL(KIND=dp) :: lup = 0.0_dp
@ -115,6 +116,7 @@ CONTAINS
!> \param ei_scale14 ...
!> \param vdw_scale14 ...
!> \param shift_cutoff ...
!> \param do_electrostatics ...
!> \param r_last_update ...
!> \param r_last_update_pbc ...
!> \param rshell_last_update_pbc ...
@ -137,7 +139,7 @@ CONTAINS
! **************************************************************************************************
SUBROUTINE fist_nonbond_env_get(fist_nonbond_env, potparm14, potparm, &
nonbonded, rlist_cut, rlist_lowsq, aup, lup, ei_scale14, vdw_scale14, &
shift_cutoff, r_last_update, r_last_update_pbc, rshell_last_update_pbc, &
shift_cutoff, do_electrostatics, r_last_update, r_last_update_pbc, rshell_last_update_pbc, &
rcore_last_update_pbc, cell_last_update, num_update, last_update, &
counter, natom_types, long_range_correction, ij_kind_full_fac, eam_data, &
quip_data, nequip_data, allegro_data, charges)
@ -147,7 +149,7 @@ CONTAINS
TYPE(fist_neighbor_type), OPTIONAL, POINTER :: nonbonded
REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: rlist_cut, rlist_lowsq
REAL(KIND=dp), OPTIONAL :: aup, lup, ei_scale14, vdw_scale14
LOGICAL, INTENT(OUT), OPTIONAL :: shift_cutoff
LOGICAL, INTENT(OUT), OPTIONAL :: shift_cutoff, do_electrostatics
TYPE(pos_type), DIMENSION(:), OPTIONAL, POINTER :: r_last_update, r_last_update_pbc, &
rshell_last_update_pbc, &
rcore_last_update_pbc
@ -189,6 +191,7 @@ CONTAINS
IF (PRESENT(vdw_scale14)) vdw_scale14 = fist_nonbond_env%vdw_scale14
IF (PRESENT(shift_cutoff)) &
shift_cutoff = fist_nonbond_env%shift_cutoff
IF (PRESENT(do_electrostatics)) do_electrostatics = fist_nonbond_env%do_electrostatics
IF (PRESENT(natom_types)) natom_types = fist_nonbond_env%natom_types
IF (PRESENT(counter)) counter = fist_nonbond_env%counter
IF (PRESENT(last_update)) last_update = fist_nonbond_env%last_update
@ -210,6 +213,7 @@ CONTAINS
!> \param ei_scale14 ...
!> \param vdw_scale14 ...
!> \param shift_cutoff ...
!> \param do_electrostatics ...
!> \param r_last_update ...
!> \param r_last_update_pbc ...
!> \param rshell_last_update_pbc ...
@ -231,7 +235,7 @@ CONTAINS
! **************************************************************************************************
SUBROUTINE fist_nonbond_env_set(fist_nonbond_env, potparm14, potparm, &
rlist_cut, rlist_lowsq, nonbonded, aup, lup, ei_scale14, vdw_scale14, &
shift_cutoff, r_last_update, r_last_update_pbc, rshell_last_update_pbc, &
shift_cutoff, do_electrostatics, r_last_update, r_last_update_pbc, rshell_last_update_pbc, &
rcore_last_update_pbc, cell_last_update, num_update, last_update, &
counter, natom_types, long_range_correction, eam_data, quip_data, &
nequip_data, allegro_data, charges)
@ -241,7 +245,7 @@ CONTAINS
REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: rlist_cut, rlist_lowsq
TYPE(fist_neighbor_type), OPTIONAL, POINTER :: nonbonded
REAL(KIND=dp), OPTIONAL :: aup, lup, ei_scale14, vdw_scale14
LOGICAL, INTENT(IN), OPTIONAL :: shift_cutoff
LOGICAL, INTENT(IN), OPTIONAL :: shift_cutoff, do_electrostatics
TYPE(pos_type), DIMENSION(:), OPTIONAL, POINTER :: r_last_update, r_last_update_pbc, &
rshell_last_update_pbc, &
rcore_last_update_pbc
@ -281,6 +285,7 @@ CONTAINS
IF (PRESENT(vdw_scale14)) fist_nonbond_env%vdw_scale14 = vdw_scale14
IF (PRESENT(shift_cutoff)) &
fist_nonbond_env%shift_cutoff = shift_cutoff
IF (PRESENT(do_electrostatics)) fist_nonbond_env%do_electrostatics = do_electrostatics
IF (PRESENT(natom_types)) fist_nonbond_env%natom_types = natom_types
IF (PRESENT(counter)) fist_nonbond_env%counter = counter
IF (PRESENT(last_update)) fist_nonbond_env%last_update = last_update
@ -296,6 +301,7 @@ CONTAINS
!> \param potparm14 ...
!> \param potparm ...
!> \param do_nonbonded ...
!> \param do_electrostatics ...
!> \param verlet_skin ...
!> \param ewald_rcut ...
!> \param ei_scale14 ...
@ -306,12 +312,12 @@ CONTAINS
!> \author Fawzi Mohamed
! **************************************************************************************************
SUBROUTINE fist_nonbond_env_create(fist_nonbond_env, atomic_kind_set, &
potparm14, potparm, do_nonbonded, verlet_skin, ewald_rcut, ei_scale14, &
vdw_scale14, shift_cutoff)
potparm14, potparm, do_nonbonded, do_electrostatics, verlet_skin, ewald_rcut, &
ei_scale14, vdw_scale14, shift_cutoff)
TYPE(fist_nonbond_env_type), INTENT(OUT) :: fist_nonbond_env
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(pair_potential_pp_type), OPTIONAL, POINTER :: potparm14, potparm
LOGICAL, INTENT(IN) :: do_nonbonded
LOGICAL, INTENT(IN) :: do_nonbonded, do_electrostatics
REAL(KIND=dp), INTENT(IN) :: verlet_skin, ewald_rcut, ei_scale14, &
vdw_scale14
LOGICAL, INTENT(IN) :: shift_cutoff
@ -333,8 +339,8 @@ CONTAINS
NULLIFY (fist_nonbond_env%allegro_data)
NULLIFY (fist_nonbond_env%charges)
CALL init_fist_nonbond_env(fist_nonbond_env, atomic_kind_set, potparm14, &
potparm, do_nonbonded, verlet_skin, ewald_rcut, ei_scale14, vdw_scale14, &
shift_cutoff)
potparm, do_nonbonded, do_electrostatics, verlet_skin, ewald_rcut, ei_scale14, &
vdw_scale14, shift_cutoff)
END SUBROUTINE fist_nonbond_env_create
! **************************************************************************************************
@ -344,6 +350,7 @@ CONTAINS
!> \param potparm14 ...
!> \param potparm ...
!> \param do_nonbonded ...
!> \param do_electrostatics ...
!> \param verlet_skin ...
!> \param ewald_rcut ...
!> \param ei_scale14 ...
@ -351,13 +358,13 @@ CONTAINS
!> \param shift_cutoff ...
! **************************************************************************************************
SUBROUTINE init_fist_nonbond_env(fist_nonbond_env, atomic_kind_set, &
potparm14, potparm, do_nonbonded, verlet_skin, ewald_rcut, ei_scale14, &
potparm14, potparm, do_nonbonded, do_electrostatics, verlet_skin, ewald_rcut, ei_scale14, &
vdw_scale14, shift_cutoff)
TYPE(fist_nonbond_env_type), INTENT(INOUT) :: fist_nonbond_env
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(pair_potential_pp_type), OPTIONAL, POINTER :: potparm14, potparm
LOGICAL, INTENT(IN) :: do_nonbonded
LOGICAL, INTENT(IN) :: do_nonbonded, do_electrostatics
REAL(KIND=dp), INTENT(IN) :: verlet_skin, ewald_rcut, ei_scale14, &
vdw_scale14
LOGICAL, INTENT(IN) :: shift_cutoff
@ -381,6 +388,7 @@ CONTAINS
NULLIFY (fist_nonbond_env%ij_kind_full_fac)
fist_nonbond_env%unit_type = "ANGSTROM"
fist_nonbond_env%do_nonbonded = do_nonbonded
fist_nonbond_env%do_electrostatics = do_electrostatics
fist_nonbond_env%lup = 0
fist_nonbond_env%aup = 0
fist_nonbond_env%ei_scale14 = ei_scale14

View file

@ -234,8 +234,8 @@ CONTAINS
fac_ei = 0.5_dp*fac_ei
fac_vdw = 0.5_dp*fac_vdw
END IF
! decide which interactions to compute
IF (do_multipoles) THEN
! decide which interactions to compute\b
IF (do_multipoles .OR. (.NOT. fist_nonbond_env%do_electrostatics)) THEN
fac_ei = 0.0_dp
END IF
IF (ipair <= neighbor_kind_pair%nscale) THEN

View file

@ -185,6 +185,7 @@ MODULE force_field_types
TYPE force_field_type
LOGICAL :: shift_cutoff, &
do_nonbonded, &
do_electrostatics, &
electrostatics, &
multiple_potential, &
ignore_missing_critical

View file

@ -164,6 +164,7 @@ CONTAINS
CALL section_vals_val_get(ff_section, "MULTIPLE_POTENTIAL", l_val=ff_type%multiple_potential)
CALL section_vals_val_get(ff_section, "DO_NONBONDED", l_val=ff_type%do_nonbonded)
CALL section_vals_val_get(ff_section, "DO_ELECTROSTATICS", l_val=ff_type%do_electrostatics)
tmp_section => section_vals_get_subs_vals(ff_section, "NONBONDED")
CALL section_vals_get(tmp_section, explicit=explicit)
IF (explicit .AND. ff_type%do_nonbonded) THEN

View file

@ -321,7 +321,7 @@ CONTAINS
ALLOCATE (fist_nonbond_env)
CALL fist_nonbond_env_create(fist_nonbond_env, atomic_kind_set, &
potparm_nonbond14, potparm_nonbond, ff_type%do_nonbonded, &
verlet_skin, ewald_rcut, ff_type%ei_scale14, &
ff_type%do_electrostatics, verlet_skin, ewald_rcut, ff_type%ei_scale14, &
ff_type%vdw_scale14, ff_type%shift_cutoff)
CALL fist_nonbond_env_set(fist_nonbond_env, charges=charges)
! Compute the electrostatic interaction cutoffs.

View file

@ -278,6 +278,14 @@ CONTAINS
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="DO_ELECTROSTATICS", &
description="Controls the computation of all the real-space "// &
"(short-range) electrostatics interactions. This does not "// &
"affect the QM/MM electrostatic coupling when turned off.", &
usage="DO_ELECTROSTATICS T", default_l_val=.TRUE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="IGNORE_MISSING_CRITICAL_PARAMS", &
description="Do not abort when critical force-field parameters "// &
"are missing. CP2K will run as if the terms containing the "// &

View file

@ -1202,8 +1202,8 @@ CONTAINS
CALL section_vals_val_get(se_section, "NEIGHBOR_LISTS%VERLET_SKIN", &
r_val=verlet_skin)
ALLOCATE (se_nonbond_env)
CALL fist_nonbond_env_create(se_nonbond_env, atomic_kind_set, &
do_nonbonded=.TRUE., verlet_skin=verlet_skin, ewald_rcut=ewald_rcut, &
CALL fist_nonbond_env_create(se_nonbond_env, atomic_kind_set, do_nonbonded=.TRUE., &
do_electrostatics=.TRUE., verlet_skin=verlet_skin, ewald_rcut=ewald_rcut, &
ei_scale14=0.0_dp, vdw_scale14=0.0_dp, shift_cutoff=.FALSE.)
! Create and Setup NDDO multipole environment
CALL nddo_mpole_setup(se_nddo_mpole, natom)

View file

@ -0,0 +1,141 @@
&FORCE_EVAL
METHOD QMMM
&DFT
BASIS_SET_FILE_NAME ../GTH_BASIS_SETS
POTENTIAL_FILE_NAME ../POTENTIAL
&MGRID
COMMENSURATE
CUTOFF 50
&END MGRID
&QS
&END QS
&SCF
SCF_GUESS atomic
&END SCF
&XC
&XC_FUNCTIONAL pade
&END XC_FUNCTIONAL
&END XC
&END DFT
&MM
&FORCEFIELD
DO_ELECTROSTATICS F
&BEND
ATOMS H O H
K 0.
THETA0 1.8
&END BEND
&BOND
ATOMS O H
K 0.
R0 1.8
&END BOND
&CHARGE
ATOM O
CHARGE -0.8476
&END CHARGE
&CHARGE
ATOM H
CHARGE 0.4238
&END CHARGE
&NONBONDED
&LENNARD-JONES
atoms O O
EPSILON 78.198
SIGMA 3.166
RCUT 11.4
&END LENNARD-JONES
&LENNARD-JONES
atoms O H
EPSILON 0.0
SIGMA 3.6705
RCUT 11.4
&END LENNARD-JONES
&LENNARD-JONES
atoms H H
EPSILON 0.0
SIGMA 3.30523
RCUT 11.4
&END LENNARD-JONES
&END NONBONDED
&END FORCEFIELD
&POISSON
&EWALD
EWALD_TYPE none
ALPHA .44
GMAX 21
&END EWALD
&END POISSON
&END MM
&QMMM
MM_POTENTIAL_FILE_NAME ../MM_POTENTIAL
&CELL
ABC 6.0 6.0 6.0
&END CELL
ECOUPL GAUSS
NOCOMPATIBILITY
&INTERPOLATOR
EPS_R 1.0e-14
EPS_X 1.0e-14
MAXITER 100
&END INTERPOLATOR
&MM_KIND H
RADIUS 0.44
&END MM_KIND
&MM_KIND O
RADIUS 0.78
&END MM_KIND
&QM_KIND H
MM_INDEX 2 3
&END QM_KIND
&QM_KIND O
MM_INDEX 1
&END QM_KIND
#
# QM_KINDS
#
#
# MM_KINDS
#
#
&END QMMM
&SUBSYS
&CELL
ABC 24.955 24.955 24.955
&END CELL
&COORD
O 0.000000 0.000000 0.000000 H2O1
H 0.000000 0.000000 1.000000 H2O1
H 0.942809 0.000000 -0.333333 H2O1
O -1.617979 -0.948062 -2.341650 H2O2
H -2.529195 -1.296822 -2.122437 H2O2
H -1.534288 -0.833088 -3.331486 H2O2
O -1.447990 2.117783 1.555094 H2O3
H -1.501128 2.645178 2.403050 H2O3
H -2.090603 1.352766 1.597519 H2O3
&END COORD
&KIND H
BASIS_SET SZV-GTH
POTENTIAL GTH-PADE-q1
&END KIND
&KIND O
BASIS_SET SZV-GTH
POTENTIAL GTH-PADE-q6
&END KIND
&TOPOLOGY
&END TOPOLOGY
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PROJECT H2O-qmmm-gauss-2
PRINT_LEVEL MEDIUM
RUN_TYPE MD
&END GLOBAL
&MOTION
&MD
ENSEMBLE NVE
STEPS 2
TIMESTEP 0.5
TEMPERATURE 298
&END MD
&END MOTION

View file

@ -27,4 +27,5 @@ H2O-qmmm-gauss-18.inp 7 6e-12
H2O-qmmm-gauss-19.inp 7 6e-12 -17.4328106395
H2O-qmmm-none-force-mixing-1.inp 1 6e-13 -33.47241414134874
H2O-qmmm-gauss-force-mixing-1.inp 1 7e-13 -50.19552620361094
H2O-qmmm-noelectrostatics.inp 1 2e-13 -16.86180642492571
#EOF