Skip to content

add minus sign to all "charge penetration" terms #79

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 4 commits into from
Jan 5, 2023
Merged

add minus sign to all "charge penetration" terms #79

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
2690192
Fix mbar example (#76)
WangXinyan940 Dec 3, 2022
db60c61
add minus sign to all "charge penetration" terms
junminchen Dec 6, 2022
2e180d1
Adding ADMP Slater-type Decompositions Test
junminchen Dec 29, 2022
29e01cc
Adding slater type function decomposition test
junminchen Dec 29, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
45 changes: 41 additions & 4 deletions dmff/generators/admp.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -599,26 +599,63 @@ def getMetaData(self):


# Here are all the short range "charge penetration" terms
# They all have the exchange form
# They all have the exchange form with minus sign
class SlaterSrEsGenerator(SlaterExGenerator):
def __init__(self, ff):
super().__init__(ff)
self.name = "SlaterSrEsForce"

def createForce(self, system, data, nonbondedMethod, nonbondedCutoff,
args):

n_atoms = len(data.atoms)
# build index map
map_atomtype = np.zeros(n_atoms, dtype=int)
for i in range(n_atoms):
atype = data.atomType[data.atoms[i]]
map_atomtype[i] = np.where(self.atomTypes == atype)[0][0]
self.map_atomtype = map_atomtype
# build covalent map
self.covalent_map = build_covalent_map(data, 6)

self._meta["cov_map"] = self.covalent_map
self._meta[self.name+"_map_atomtype"] = self.map_atomtype

pot_fn_sr = generate_pairwise_interaction(slater_sr_kernel,
static_args={})

def potential_fn(positions, box, pairs, params):
params = params[self.name]
mScales = params["mScales"]
a_list = params["A"][map_atomtype]
b_list = params["B"][map_atomtype] / 10 # nm^-1 to A^-1

# add minus sign
return - pot_fn_sr(positions, box, pairs, mScales, a_list, b_list)

self._jaxPotential = potential_fn
# self._top_data = data

def getJaxPotential(self):
return self._jaxPotential

def getMetaData(self):
return self._meta


class SlaterSrPolGenerator(SlaterExGenerator):
class SlaterSrPolGenerator(SlaterSrEsGenerator):
def __init__(self, ff):
super().__init__(ff)
self.name = "SlaterSrPolForce"


class SlaterSrDispGenerator(SlaterExGenerator):
class SlaterSrDispGenerator(SlaterSrEsGenerator):
def __init__(self, ff):
super().__init__(ff)
self.name = "SlaterSrDispForce"


class SlaterDhfGenerator(SlaterExGenerator):
class SlaterDhfGenerator(SlaterSrEsGenerator):
def __init__(self, ff):
super().__init__(ff)
self.name = "SlaterDhfForce"
Expand Down
6 changes: 6 additions & 0 deletions docs/user_guide/installation.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,12 @@ pip install jax==0.3.17
```bash
pip install jax-md==0.2.0
```
+ Install [mdtraj](https://github.com/mdtraj/mdtraj), [optax](https://github.com/deepmind/optax) and [pymbar](https://github.com/choderalab/pymbar):
```bash
conda install -c conda-forge mdtraj==1.9.7
pip install optax==0.1.3
pip install pymbar==4.0.1
```
+ Install [OpenMM](https://openmm.org/):
```bash
conda install -c conda-forge openmm==7.7.0
Expand Down
5 changes: 0 additions & 5 deletions examples/mbar/ben-prm.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,11 +18,6 @@
<HarmonicBondForce>
<Bond type1="ca" type2="ca" length="0.24075506225207396" k="30000.0"/>
</HarmonicBondForce>
<NonbondedForce coulomb14scale="0.8333333333333334" lj14scale="0.5">
<UseAttributeFromResidue name="charge"/>
<Atom type="ca" sigma="1.0" epsilon="0.0"/>
<Atom type="ha" sigma="1.0" epsilon="0.0"/>
</NonbondedForce>
<LennardJonesForce lj14scale="0.5" useDispersionCorrection="False">
<Atom type="ca" sigma="0.45" epsilon="1.0"/>
</LennardJonesForce>
Expand Down
58 changes: 58 additions & 0 deletions examples/peg_slater_isa/check_calc.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,58 @@
#!/usr/bin/env python
import sys
import numpy as np
import openmm
from openmm import *
from openmm.app import *
from openmm.unit import *
import jax
import jax_md
import jax.numpy as jnp
import dmff
from dmff.api import Hamiltonian
from dmff.common import nblist
import pickle
import time
from jax import value_and_grad, jit

if __name__ == '__main__':
ff = 'peg.xml'
pdb_AB = PDBFile('peg2.pdb')
H_AB = Hamiltonian(ff)
rc = 15
# get potential functions
pots_AB = H_AB.createPotential(pdb_AB.topology, nonbondedCutoff=rc*angstrom, nonbondedMethod=CutoffPeriodic, ethresh=1e-4)
pot_pme_AB = pots_AB.dmff_potentials['ADMPPmeForce']
pot_disp_AB = pots_AB.dmff_potentials['ADMPDispPmeForce']
pot_ex_AB = pots_AB.dmff_potentials['SlaterExForce']
pot_sr_es_AB = pots_AB.dmff_potentials['SlaterSrEsForce']
pot_sr_pol_AB = pots_AB.dmff_potentials['SlaterSrPolForce']
pot_sr_disp_AB = pots_AB.dmff_potentials['SlaterSrDispForce']
pot_dhf_AB = pots_AB.dmff_potentials['SlaterDhfForce']
pot_dmp_es_AB = pots_AB.dmff_potentials['QqTtDampingForce']
pot_dmp_disp_AB = pots_AB.dmff_potentials['SlaterDampingForce']

paramtree = H_AB.getParameters()

# init positions used to set up neighbor list
pos_AB0 = jnp.array(pdb_AB.positions._value) * 10
n_atoms = len(pos_AB0)
box = jnp.array(pdb_AB.topology.getPeriodicBoxVectors()._value) * 10

# nn list initial allocation
nbl_AB = nblist.NeighborList(box, rc, H_AB.getGenerators()[0].covalent_map)
nbl_AB.allocate(pos_AB0)
pairs_AB = nbl_AB.pairs
pairs_AB = pairs_AB[pairs_AB[:, 0] < pairs_AB[:, 1]]

pos_AB = jnp.array(pos_AB0)
E_es = pot_pme_AB(pos_AB, box, pairs_AB, paramtree)
E_disp = pot_disp_AB(pos_AB, box, pairs_AB, paramtree)
E_ex_AB = pot_ex_AB(pos_AB, box, pairs_AB, paramtree)
E_sr_es = pot_sr_es_AB(pos_AB, box, pairs_AB, paramtree)
E_sr_pol = pot_sr_pol_AB(pos_AB, box, pairs_AB, paramtree)
E_sr_disp = pot_sr_disp_AB(pos_AB, box, pairs_AB, paramtree)
E_dhf = pot_dhf_AB(pos_AB, box, pairs_AB, paramtree)
E_dmp_es = pot_dmp_es_AB(pos_AB, box, pairs_AB, paramtree)
E_dmp_disp = pot_dmp_disp_AB(pos_AB, box, pairs_AB, paramtree)
print(E_es, E_disp, E_ex_AB, E_sr_es, E_sr_pol, E_sr_disp, E_dhf, E_dmp_es, E_dmp_disp)
145 changes: 145 additions & 0 deletions examples/peg_slater_isa/peg.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,145 @@
<ForceField>
<AtomTypes>
<Type class="CT" element="C" mass="12.0107" name="1" />
<Type class="HC" element="H" mass="1.00784" name="2" />
<Type class="OS" element="O" mass="15.999" name="3" />
<Type class="CT" element="C" mass="12.0107" name="4" />
<Type class="HC" element="H" mass="1.00784" name="5" />
</AtomTypes>
<Residues>
<Residue name="TER">
<Atom name="C00" type="1" />
<Atom name="H01" type="2" />
<Atom name="H02" type="2" />
<Atom name="O03" type="3" />
<Atom name="C04" type="4" />
<Atom name="H05" type="5" />
<Atom name="H06" type="5" />
<Atom name="H07" type="5" />
<Bond from="0" to="1" />
<Bond from="0" to="2" />
<Bond from="0" to="3" />
<Bond from="3" to="4" />
<Bond from="4" to="5" />
<Bond from="4" to="6" />
<Bond from="4" to="7" />
<ExternalBond atomName="C00" />
</Residue>
<Residue name="INT">
<Atom name="C00" type="1" />
<Atom name="H01" type="2" />
<Atom name="H02" type="2" />
<Atom name="O03" type="3" />
<Atom name="C04" type="1" />
<Atom name="H05" type="2" />
<Atom name="H06" type="2" />
<Bond from="0" to="1" />
<Bond from="0" to="2" />
<Bond from="0" to="3" />
<Bond from="3" to="4" />
<Bond from="4" to="5" />
<Bond from="4" to="6" />
<ExternalBond atomName="C00" />
<ExternalBond atomName="C04" />
</Residue>
</Residues>
<ADMPPmeForce lmax="2"
mScale12="0.00" mScale13="0.00" mScale14="0.00" mScale15="0.00" mScale16="0.00"
pScale12="0.00" pScale13="0.00" pScale14="0.00" pScale15="0.00" pScale16="0.00"
dScale12="1.00" dScale13="1.00" dScale14="1.00" dScale15="1.00" dScale16="1.00">
<Atom type="1" kx="1" kz="3"
c0="0.14574378"
dX="0.00158195" dY="0.00011841" dZ="0.00693851"
qXX="-0.00005399" qXY="0.00000367" qYY="-0.00005356" qXZ="-0.00002917" qYZ="0.00000185" qZZ="0.00010755" />
<Atom type="2" kx="2" kz="1"
c0="0.02178918"
dX="-0.00043629" dY="0.00004082" dZ="-0.00219400"
qXX="-0.00002509" qXY="0.00000025" qYY="0.00000855" qXZ="0.00000249" qYZ="-0.00000070" qZZ="0.00001654" />
<Atom type="3" kx="4" kz="-1"
c0="-0.34690552"
dX="0.00074502" dY="-0.00008518" dZ="0.01025277"
qXX="0.00033647" qXY="-0.00000009" qYY="-0.00023636" qXZ="0.00000021" qYZ="0.00000353" qZZ="-0.00010012" />
<Atom type="3" kx="1" kz="-1"
c0="-0.34690552"
dX="0.00074502" dY="-0.00008518" dZ="0.01025277"
qXX="0.00033647" qXY="-0.00000009" qYY="-0.00023636" qXZ="0.00000021" qYZ="0.00000353" qZZ="-0.00010012" />
<Atom type="4" kx="5" kz="3"
c0="0.00942772"
dX="0.00001267" dY="0.00003435" dZ="0.00815618"
qXX="-0.00006481" qXY="0.00000355" qYY="-0.00006754" qXZ="0.00000077" qYZ="0.00000115" qZZ="0.00013235" />
<Atom type="5" kx="3" kz="4"
c0="0.04964274"
dX="0.00046553" dY="-0.00001268" dZ="-0.00206894"
qXX="0.00001716" qXY="-0.00000005" qYY="-0.00003286" qXZ="0.00000060" qYZ="0.00000033" qZZ="0.00001570" />
<Polarize type="1" polarizabilityXX="1.072970e-03" polarizabilityYY="1.072970e-03" polarizabilityZZ="1.072970e-03" thole="0.33" />
<Polarize type="2" polarizabilityXX="3.680091e-04" polarizabilityYY="3.680091e-04" polarizabilityZZ="3.680091e-04" thole="0.33" />
<Polarize type="3" polarizabilityXX="6.192140e-04" polarizabilityYY="6.192140e-04" polarizabilityZZ="6.192140e-04" thole="0.33" />
<Polarize type="4" polarizabilityXX="1.099428e-03" polarizabilityYY="1.099428e-03" polarizabilityZZ="1.099428e-03" thole="0.33" />
<Polarize type="5" polarizabilityXX="3.365314e-04" polarizabilityYY="3.365314e-04" polarizabilityZZ="3.365314e-04" thole="0.33" />
</ADMPPmeForce>
<ADMPDispPmeForce
mScale12="0.00" mScale13="0.00" mScale14="0.00" mScale15="0.00" mScale16="0.00" >
<Atom type="1" C6="1.507393e-03" C8="1.241793e-04" C10="4.890285e-06" />
<Atom type="2" C6="1.212045e-04" C8="4.577425e-06" C10="8.708729e-08" />
<Atom type="3" C6="7.159800e-04" C8="5.846871e-05" C10="2.282115e-06" />
<Atom type="4" C6="1.523005e-03" C8="1.127912e-04" C10="4.005600e-06" />
<Atom type="5" C6="1.136931e-04" C8="4.123377e-06" C10="7.495037e-08" />
</ADMPDispPmeForce>
<SlaterExForce
mScale12="0.00" mScale13="0.00" mScale14="0.00" mScale15="0.00" mScale16="0.00" >
<Atom type="1" A="1" B="3.977508e+01" />
<Atom type="2" A="1" B="4.596271e+01" />
<Atom type="3" A="1" B="4.637414e+01" />
<Atom type="4" A="1" B="3.831504e+01" />
<Atom type="5" A="1" B="4.632228e+01" />
</SlaterExForce>
<SlaterSrEsForce
mScale12="0.00" mScale13="0.00" mScale14="0.00" mScale15="0.00" mScale16="0.00" >
<Atom type="1" A="1" B="3.977508e+01" Q="0.14574378"/>
<Atom type="2" A="1" B="4.596271e+01" Q="0.02178918"/>
<Atom type="3" A="1" B="4.637414e+01" Q="-0.34690552"/>
<Atom type="4" A="1" B="3.831504e+01" Q="0.00942772"/>
<Atom type="5" A="1" B="4.632228e+01" Q="0.04964274"/>
</SlaterSrEsForce>
<SlaterSrPolForce
mScale12="0.00" mScale13="0.00" mScale14="0.00" mScale15="0.00" mScale16="0.00" >
<Atom type="1" A="1" B="3.977508e+01" />
<Atom type="2" A="1" B="4.596271e+01" />
<Atom type="3" A="1" B="4.637414e+01" />
<Atom type="4" A="1" B="3.831504e+01" />
<Atom type="5" A="1" B="4.632228e+01" />
</SlaterSrPolForce>
<SlaterSrDispForce
mScale12="0.00" mScale13="0.00" mScale14="0.00" mScale15="0.00" mScale16="0.00" >
<Atom type="1" A="1" B="3.977508e+01" />
<Atom type="2" A="1" B="4.596271e+01" />
<Atom type="3" A="1" B="4.637414e+01" />
<Atom type="4" A="1" B="3.831504e+01" />
<Atom type="5" A="1" B="4.632228e+01" />
</SlaterSrDispForce>
<SlaterDhfForce
mScale12="0.00" mScale13="0.00" mScale14="0.00" mScale15="0.00" mScale16="0.00" >
<Atom type="1" A="1" B="3.977508e+01" />
<Atom type="2" A="1" B="4.596271e+01" />
<Atom type="3" A="1" B="4.637414e+01" />
<Atom type="4" A="1" B="3.831504e+01" />
<Atom type="5" A="1" B="4.632228e+01" />
</SlaterDhfForce>
<QqTtDampingForce
mScale12="0.00" mScale13="0.00" mScale14="0.00" mScale15="0.00" mScale16="0.00" >
<Atom type="1" B="3.977508e+01" Q="0.14574378"/>
<Atom type="2" B="4.596271e+01" Q="0.02178918"/>
<Atom type="3" B="4.637414e+01" Q="-0.34690552"/>
<Atom type="4" B="3.831504e+01" Q="0.00942772"/>
<Atom type="5" B="4.632228e+01" Q="0.04964274"/>
</QqTtDampingForce>
<SlaterDampingForce
mScale12="0.00" mScale13="0.00" mScale14="0.00" mScale15="0.00" mScale16="0.00" >
<Atom type="1" B="3.977508e+01" C6="1.507393e-03" C8="1.241793e-04" C10="4.890285e-06" />
<Atom type="2" B="4.596271e+01" C6="1.212045e-04" C8="4.577425e-06" C10="8.708729e-08" />
<Atom type="3" B="4.637414e+01" C6="7.159800e-04" C8="5.846871e-05" C10="2.282115e-06" />
<Atom type="4" B="3.831504e+01" C6="1.523005e-03" C8="1.127912e-04" C10="4.005600e-06" />
<Atom type="5" B="4.632228e+01" C6="1.136931e-04" C8="4.123377e-06" C10="7.495037e-08" />
</SlaterDampingForce>
</ForceField>

Loading