Skip to content

GitLab

Commit 0a214970 authored by peastman's avatar peastman
Browse files

Support all values of NBXMod in CHARMM parameter files

parent b71e92d9
Hide whitespace changes
Inline Side-by-side
Showing with 117 additions and 36 deletions
wrappers/python/simtk/openmm/app/charmmparameterset.py
View file @ 0a214970
...@@ -116,6 +116,7 @@ class CharmmParameterSet(object): ...@@ -116,6 +116,7 @@ class CharmmParameterSet(object):
self.nbfix_types = dict() self.nbfix_types = dict()
self.nbthole_types = dict() self.nbthole_types = dict()
self.parametersets = [] self.parametersets = []
self.nbxmod = 5
# Load all of the files # Load all of the files
tops, pars, strs = [], [], [] tops, pars, strs = [], [], []
...@@ -228,8 +229,17 @@ class CharmmParameterSet(object): ...@@ -228,8 +229,17 @@ class CharmmParameterSet(object):
nonbonded_types = dict() # Holder nonbonded_types = dict() # Holder
parameterset = None parameterset = None
read_first_nonbonded = False read_first_nonbonded = False
previous = ''
for line in f: for line in f:
line = line.strip() line = previous+line.strip()
previous = ''
if line.endswith('-'):
# This will be continued on the next line.
previous = line[:-1]
continue
if line.startswith('!'):
# This is a comment.
continue
if not line: if not line:
# This is a blank line # This is a blank line
continue continue
...@@ -258,6 +268,12 @@ class CharmmParameterSet(object): ...@@ -258,6 +268,12 @@ class CharmmParameterSet(object):
if line.startswith('NONBONDED'): if line.startswith('NONBONDED'):
read_first_nonbonded = False read_first_nonbonded = False
section = 'NONBONDED' section = 'NONBONDED'
fields = line.upper().split()
if 'NBXMOD' in fields:
nbxmod = int(fields[fields.index('NBXMOD')+1])
if nbxmod not in list(range(-5, 6)):
raise CharmmFileError('Unsupported value for NBXMOD: %d' % nbxmod)
self.nbxmod = nbxmod
continue continue
if line.startswith('NBFIX'): if line.startswith('NBFIX'):
section = 'NBFIX' section = 'NBFIX'
......
wrappers/python/simtk/openmm/app/charmmpsffile.py
View file @ 0a214970
...@@ -1323,26 +1323,28 @@ class CharmmPsfFile(object): ...@@ -1323,26 +1323,28 @@ class CharmmPsfFile(object):
# Add 1-4 interactions # Add 1-4 interactions
excluded_atom_pairs = set() # save these pairs so we don't zero them out excluded_atom_pairs = set() # save these pairs so we don't zero them out
sigma_scale = 2**(-1/6) sigma_scale = 2**(-1/6)
for tor in self.dihedral_parameter_list: nbxmod = abs(params.nbxmod)
# First check to see if atoms 1 and 4 are already excluded because if nbxmod == 5:
# they are 1-2 or 1-3 pairs (would happen in 6-member rings or for tor in self.dihedral_parameter_list:
# fewer). Then check that they're not already added as exclusions # First check to see if atoms 1 and 4 are already excluded because
if tor.atom1 in tor.atom4.bond_partners: continue # they are 1-2 or 1-3 pairs (would happen in 6-member rings or
if tor.atom1 in tor.atom4.angle_partners: continue # fewer). Then check that they're not already added as exclusions
key = min((tor.atom1.idx, tor.atom4.idx), if tor.atom1 in tor.atom4.bond_partners: continue
(tor.atom4.idx, tor.atom1.idx)) if tor.atom1 in tor.atom4.angle_partners: continue
if key in excluded_atom_pairs: continue # multiterm... key = min((tor.atom1.idx, tor.atom4.idx),
charge_prod = (tor.atom1.charge * tor.atom4.charge) (tor.atom4.idx, tor.atom1.idx))
epsilon = (sqrt(abs(tor.atom1.type.epsilon_14) * ene_conv * if key in excluded_atom_pairs: continue # multiterm...
abs(tor.atom4.type.epsilon_14) * ene_conv)) charge_prod = (tor.atom1.charge * tor.atom4.charge)
sigma = (tor.atom1.type.rmin_14 + tor.atom4.type.rmin_14) * ( epsilon = (sqrt(abs(tor.atom1.type.epsilon_14) * ene_conv *
length_conv * sigma_scale) abs(tor.atom4.type.epsilon_14) * ene_conv))
force.addException(tor.atom1.idx, tor.atom4.idx, sigma = (tor.atom1.type.rmin_14 + tor.atom4.type.rmin_14) * (
charge_prod, sigma, epsilon) length_conv * sigma_scale)
excluded_atom_pairs.add( force.addException(tor.atom1.idx, tor.atom4.idx,
min((tor.atom1.idx, tor.atom4.idx), charge_prod, sigma, epsilon)
(tor.atom4.idx, tor.atom1.idx)) excluded_atom_pairs.add(
) min((tor.atom1.idx, tor.atom4.idx),
(tor.atom4.idx, tor.atom1.idx))
)
# Add excluded atoms # Add excluded atoms
# Drude and lonepairs will be excluded based on their parent atoms # Drude and lonepairs will be excluded based on their parent atoms
...@@ -1368,21 +1370,24 @@ class CharmmPsfFile(object): ...@@ -1368,21 +1370,24 @@ class CharmmPsfFile(object):
force.addException(excludeterm[j], excludeterm[i], 0.0, 0.1, 0.0) force.addException(excludeterm[j], excludeterm[i], 0.0, 0.1, 0.0)
# Exclude all bonds and angles, as well as the lonepair/Drude attached onto them # Exclude all bonds and angles, as well as the lonepair/Drude attached onto them
for atom in self.atom_list: for atom in self.atom_list:
for atom2 in atom.bond_partners: if nbxmod > 1:
if atom2.idx > atom.idx: for atom2 in atom.bond_partners:
for excludeatom in [atom.idx]+parent_exclude_list[atom.idx]: if atom2.idx > atom.idx:
for excludeatom2 in [atom2.idx]+parent_exclude_list[atom2.idx]: for excludeatom in [atom.idx]+parent_exclude_list[atom.idx]:
force.addException(excludeatom, excludeatom2, 0.0, 0.1, 0.0) for excludeatom2 in [atom2.idx]+parent_exclude_list[atom2.idx]:
for atom2 in atom.angle_partners: force.addException(excludeatom, excludeatom2, 0.0, 0.1, 0.0)
if atom2.idx > atom.idx: if nbxmod > 2:
for excludeatom in [atom.idx]+parent_exclude_list[atom.idx]: for atom2 in atom.angle_partners:
for excludeatom2 in [atom2.idx]+parent_exclude_list[atom2.idx]: if atom2.idx > atom.idx:
force.addException(excludeatom, excludeatom2, 0.0, 0.1, 0.0) for excludeatom in [atom.idx]+parent_exclude_list[atom.idx]:
for atom2 in atom.dihedral_partners: for excludeatom2 in [atom2.idx]+parent_exclude_list[atom2.idx]:
if atom2.idx <= atom.idx: continue force.addException(excludeatom, excludeatom2, 0.0, 0.1, 0.0)
if ((atom.idx, atom2.idx) in excluded_atom_pairs): if nbxmod > 3:
continue for atom2 in atom.dihedral_partners:
force.addException(atom.idx, atom2.idx, 0.0, 0.1, 0.0) if atom2.idx <= atom.idx: continue
if ((atom.idx, atom2.idx) in excluded_atom_pairs):
continue
force.addException(atom.idx, atom2.idx, 0.0, 0.1, 0.0)
system.addForce(force) system.addForce(force)
# Add Drude particles (Drude force) # Add Drude particles (Drude force)
......
wrappers/python/tests/TestCharmmFiles.py
View file @ 0a214970
...@@ -5,7 +5,12 @@ from simtk.openmm import * ...@@ -5,7 +5,12 @@ from simtk.openmm import *
from simtk.unit import * from simtk.unit import *
import simtk.openmm.app.element as elem import simtk.openmm.app.element as elem
import math import math
import tempfile
import warnings import warnings
if sys.version_info >= (3,0):
from io import StringIO
else:
from cStringIO import StringIO
class TestCharmmFiles(unittest.TestCase): class TestCharmmFiles(unittest.TestCase):
...@@ -324,6 +329,24 @@ class TestCharmmFiles(unittest.TestCase): ...@@ -324,6 +329,24 @@ class TestCharmmFiles(unittest.TestCase):
psf.setBox(3.0584*nanometers,3.0584*nanometers,3.0584*nanometers) psf.setBox(3.0584*nanometers,3.0584*nanometers,3.0584*nanometers)
psf.createSystem(parameters, nonbondedMethod=PME) psf.createSystem(parameters, nonbondedMethod=PME)
def test_NBXMod(self):
"""Test that all values of NBXMod are interpreted correctly."""
crd = CharmmCrdFile('systems/ala_ala_ala.crd')
with open('systems/charmm22.par') as parfile:
par = parfile.read()
# The following values were computed with CHARMM.
modeEnergy = {0: 754318.20507, 1: 754318.20507, 2: 908.35224, 3: 59.65279, 4: -241.12856, 5: 39.13169}
for nbxmod in range(-5, 6):
with tempfile.NamedTemporaryFile(suffix='.par', mode='w') as parfile:
parfile.write(par.replace('nbxmod 5', 'nbxmod %d' % nbxmod))
parfile.flush()
params = CharmmParameterSet('systems/charmm22.rtf', parfile.name)
system = self.psf_c.createSystem(params, nonbondedMethod=NoCutoff)
context = Context(system, VerletIntegrator(1*femtoseconds), Platform.getPlatformByName('Reference'))
context.setPositions(crd.positions)
energy = context.getState(getEnergy=True).getPotentialEnergy().value_in_unit(kilocalories_per_mole)
self.assertAlmostEqual(energy, modeEnergy[abs(nbxmod)], delta=1e-3*abs(energy))
if __name__ == '__main__': if __name__ == '__main__':
unittest.main() unittest.main()
......
wrappers/python/tests/systems/ala_ala_ala.crd 0 → 100644
View file @ 0a214970
* NONE *
* DATE: 9/ 5/19 10:32: 5 CREATED BY USER: peastman
*
33
1 1 ALA N 0.02400 -0.10300 -0.10100 AAL 1 0.00000
2 1 ALA HT1 0.02700 -1.13200 -0.23900 AAL 1 0.00000
3 1 ALA HT2 -0.80500 0.16300 0.47100 AAL 1 0.00000
4 1 ALA HT3 -0.05900 0.38400 -1.01900 AAL 1 0.00000
5 1 ALA CA 1.24700 0.37500 0.63600 AAL 1 0.00000
6 1 ALA HA 0.81400 0.86100 1.49500 AAL 1 0.00000
7 1 ALA CB 2.05700 -0.77200 1.28900 AAL 1 0.00000
8 1 ALA HB1 3.13600 -0.75200 1.03200 AAL 1 0.00000
9 1 ALA HB2 1.99000 -0.64100 2.39500 AAL 1 0.00000
10 1 ALA HB3 1.65600 -1.78200 1.06300 AAL 1 0.00000
11 1 ALA C 1.95600 1.57900 0.03600 AAL 1 0.00000
12 1 ALA O 1.21900 2.52500 -0.20100 AAL 1 0.00000
13 2 ALA N 3.28900 1.63100 -0.20200 AAL 2 0.00000
14 2 ALA HN 3.93900 0.86800 -0.17400 AAL 2 0.00000
15 2 ALA CA 3.99000 2.90900 -0.21500 AAL 2 0.00000
16 2 ALA HA 3.74200 3.44000 0.69500 AAL 2 0.00000
17 2 ALA CB 3.66200 3.80200 -1.43400 AAL 2 0.00000
18 2 ALA HB1 4.19200 4.77800 -1.35800 AAL 2 0.00000
19 2 ALA HB2 3.95600 3.31100 -2.38200 AAL 2 0.00000
20 2 ALA HB3 2.57700 4.02700 -1.46700 AAL 2 0.00000
21 2 ALA C 5.48700 2.65400 -0.12800 AAL 2 0.00000
22 2 ALA O 5.88900 1.48900 -0.13700 AAL 2 0.00000
23 3 ALA N 6.27500 3.73300 -0.03700 AAL 3 0.00000
24 3 ALA HN 5.96300 4.69100 -0.02800 AAL 3 0.00000
25 3 ALA CA 7.70700 3.80200 0.06800 AAL 3 0.00000
26 3 ALA HA 8.16000 3.41800 -0.83300 AAL 3 0.00000
27 3 ALA CB 8.23300 3.09300 1.33300 AAL 3 0.00000
28 3 ALA HB1 9.34200 3.14900 1.35600 AAL 3 0.00000
29 3 ALA HB2 7.83500 3.59300 2.24000 AAL 3 0.00000
30 3 ALA HB3 7.92300 2.03000 1.33200 AAL 3 0.00000
31 3 ALA C 8.01800 5.32300 0.13600 AAL 3 0.00000
32 3 ALA OT1 7.03200 6.11900 0.12700 AAL 3 0.00000
33 3 ALA OT2 9.21900 5.69200 0.18800 AAL 3 0.00000
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment