Skip to content

GitLab

Commit 5a81b23d authored by peastman's avatar peastman
Browse files

Fixes to GromacsTopFile with different combining rules

parent 1665962a
Hide whitespace changes
Inline Side-by-side
Showing with 36 additions and 36 deletions
wrappers/python/simtk/openmm/app/gromacstopfile.py
View file @ 5a81b23d
...@@ -596,15 +596,13 @@ class GromacsTopFile(object): ...@@ -596,15 +596,13 @@ class GromacsTopFile(object):
sys.setDefaultPeriodicBoxVectors(*boxVectors) sys.setDefaultPeriodicBoxVectors(*boxVectors)
elif nonbondedMethod in (ff.CutoffPeriodic, ff.Ewald, ff.PME, ff.LJPME): elif nonbondedMethod in (ff.CutoffPeriodic, ff.Ewald, ff.PME, ff.LJPME):
raise ValueError('Illegal nonbonded method for a non-periodic system') raise ValueError('Illegal nonbonded method for a non-periodic system')
if self._defaults[1] == '1': nb = mm.NonbondedForce()
# factor of 138.935456 for coulomb taken from tests/TestCustomNonbondedForce.h:574
nb = mm.CustomNonbondedForce('138.935456*q1*q2/r-sqrt(C1*C2)/r^6+sqrt(A1*A2)/r^12')
nb.addPerParticleParameter('q')
nb.addPerParticleParameter('C')
nb.addPerParticleParameter('A')
elif self._defaults[1] == '2':
nb = mm.NonbondedForce()
sys.addForce(nb) sys.addForce(nb)
if self._defaults[1] == '1':
lj = mm.CustomNonbondedForce('sqrt(A1*A2)/r^12-sqrt(C1*C2)/r^6')
lj.addPerParticleParameter('C')
lj.addPerParticleParameter('A')
sys.addForce(lj)
if implicitSolvent is OBC2: if implicitSolvent is OBC2:
gb = mm.GBSAOBCForce() gb = mm.GBSAOBCForce()
gb.setSoluteDielectric(soluteDielectric) gb.setSoluteDielectric(soluteDielectric)
...@@ -622,7 +620,7 @@ class GromacsTopFile(object): ...@@ -622,7 +620,7 @@ class GromacsTopFile(object):
mapIndices = {} mapIndices = {}
bondIndices = [] bondIndices = []
topologyAtoms = list(self.topology.atoms()) topologyAtoms = list(self.topology.atoms())
exceptions = [] exclusions = []
pairs = [] pairs = []
fudgeQQ = float(self._defaults[4]) fudgeQQ = float(self._defaults[4])
fudgeLJ = float(self._defaults[3]) fudgeLJ = float(self._defaults[3])
...@@ -793,7 +791,7 @@ class GromacsTopFile(object): ...@@ -793,7 +791,7 @@ class GromacsTopFile(object):
phi0 = float(params[5]) phi0 = float(params[5])
if k != 0: if k != 0:
if harmonicTorsion is None: if harmonicTorsion is None:
harmonicTorsion = mm.CustomTorsionForce('0.5*k*(thetap-theta0)^2; thetap = step(-(theta-theta0+pi))*2*pi+theta+step(theta-theta0-pi)*(-2*pi); pi = 3.14159') harmonicTorsion = mm.CustomTorsionForce('0.5*k*(thetap-theta0)^2; thetap = step(-(theta-theta0+pi))*2*pi+theta+step(theta-theta0-pi)*(-2*pi); pi = %.15g' % math.pi)
harmonicTorsion.addPerTorsionParameter('theta0') harmonicTorsion.addPerTorsionParameter('theta0')
harmonicTorsion.addPerTorsionParameter('k') harmonicTorsion.addPerTorsionParameter('k')
sys.addForce(harmonicTorsion) sys.addForce(harmonicTorsion)
...@@ -849,7 +847,8 @@ class GromacsTopFile(object): ...@@ -849,7 +847,8 @@ class GromacsTopFile(object):
q = float(params[4]) q = float(params[4])
if self._defaults[1] == '1': if self._defaults[1] == '1':
nb.addParticle([q, float(params[6]), float(params[7])]) nb.addParticle(q, 1.0, 0.0)
lj.addParticle([float(params[6]), float(params[7])])
elif self._defaults[1] == '2': elif self._defaults[1] == '2':
nb.addParticle(q, float(params[6]), float(params[7])) nb.addParticle(q, float(params[6]), float(params[7]))
...@@ -880,44 +879,37 @@ class GromacsTopFile(object): ...@@ -880,44 +879,37 @@ class GromacsTopFile(object):
continue # We'll use the automatically generated parameters continue # We'll use the automatically generated parameters
atom1params = nb.getParticleParameters(baseAtomIndex+atoms[0]) atom1params = nb.getParticleParameters(baseAtomIndex+atoms[0])
atom2params = nb.getParticleParameters(baseAtomIndex+atoms[1]) atom2params = nb.getParticleParameters(baseAtomIndex+atoms[1])
# enable tracking pairs and exceptions separately pairs.append((baseAtomIndex+atoms[0], baseAtomIndex+atoms[1], atom1params[0]*atom2params[0]*fudgeQQ, params[0], params[1]))
if self._defaults[1] == '1':
pairs.append((baseAtomIndex+atoms[0], baseAtomIndex+atoms[1], atom1params[0]*atom2params[0]*fudgeQQ, params[0], params[1]))
elif self._defaults[1] == '2':
exceptions.append((baseAtomIndex+atoms[0], baseAtomIndex+atoms[1], atom1params[0]*atom2params[0]*fudgeQQ, params[0], params[1]))
for fields in moleculeType.exclusions: for fields in moleculeType.exclusions:
atoms = [int(x)-1 for x in fields] atoms = [int(x)-1 for x in fields]
for atom in atoms[1:]: for atom in atoms[1:]:
if atom > atoms[0]: if atom > atoms[0]:
if self._defaults[1] == '1': exclusions.append((baseAtomIndex+atoms[0], baseAtomIndex+atom))
exceptions.append((baseAtomIndex+atoms[0], baseAtomIndex+atom))
elif self._defaults[1] == '2':
exceptions.append((baseAtomIndex+atoms[0], baseAtomIndex+atom, 0, 0, 0))
# Create nonbonded exceptions. # Create nonbonded exceptions.
nb.createExceptionsFromBonds(bondIndices, fudgeQQ, fudgeLJ)
for exclusion in exclusions:
nb.addException(exclusion[0], exclusion[1], 0.0, 1.0, 0.0, True)
if self._defaults[1] == '1': if self._defaults[1] == '1':
# implement named exceptions (generally used for 1-4 interactions) # We're using a CustomNonbondedForce for LJ interactions, so also create a CustomBondForce
# to handle the exceptions.
pair_bond = mm.CustomBondForce('138.935456*q/r-C/r^6+A/r^12') pair_bond = mm.CustomBondForce('138.935456*q/r-C/r^6+A/r^12')
pair_bond.addPerBondParameter('q') pair_bond.addPerBondParameter('q')
pair_bond.addPerBondParameter('C') pair_bond.addPerBondParameter('C')
pair_bond.addPerBondParameter('A') pair_bond.addPerBondParameter('A')
sys.addForce(pair_bond) sys.addForce(pair_bond)
lj.createExclusionsFromBonds(bondIndices, 3)
# generate 1-2, 1-3, and 1-4 exclusions
nb.createExclusionsFromBonds(bondIndices, 3)
for pair in pairs: for pair in pairs:
pair_bond.addBond(pair[0], pair[1], [pair[2],float(pair[3])*fudgeLJ, float(pair[4])*fudgeLJ]) nb.addException(pair[0], pair[1], pair[2], 1.0, 0.0, True)
pair_bond.addBond(pair[0], pair[1], [pair[2],float(pair[3]), float(pair[4])])
# create requested exclusions for exclusion in exclusions:
for exception in exceptions: lj.addExclusion(exclusion[0], exclusion[1])
nb.addExclusion(exception[0], exception[1])
elif self._defaults[1] == '2': elif self._defaults[1] == '2':
nb.createExceptionsFromBonds(bondIndices, fudgeQQ, fudgeLJ) for pair in pairs:
for exception in exceptions: nb.addException(pair[0], pair[1], pair[2], float(pair[3]), float(pair[4]), True)
nb.addException(exception[0], exception[1], exception[2], float(exception[3])*fudgeLJ, float(exception[4]), True)
# Finish configuring the NonbondedForce. # Finish configuring the NonbondedForce.
...@@ -929,8 +921,16 @@ class GromacsTopFile(object): ...@@ -929,8 +921,16 @@ class GromacsTopFile(object):
ff.LJPME:mm.NonbondedForce.LJPME} ff.LJPME:mm.NonbondedForce.LJPME}
nb.setNonbondedMethod(methodMap[nonbondedMethod]) nb.setNonbondedMethod(methodMap[nonbondedMethod])
nb.setCutoffDistance(nonbondedCutoff) nb.setCutoffDistance(nonbondedCutoff)
if self._defaults[1] == '2': nb.setEwaldErrorTolerance(ewaldErrorTolerance)
nb.setEwaldErrorTolerance(ewaldErrorTolerance) if self._defaults[1] == '1':
methodMap = {ff.NoCutoff:mm.CustomNonbondedForce.NoCutoff,
ff.CutoffNonPeriodic:mm.CustomNonbondedForce.CutoffNonPeriodic,
ff.CutoffPeriodic:mm.CustomNonbondedForce.CutoffPeriodic,
ff.Ewald:mm.CustomNonbondedForce.CutoffPeriodic,
ff.PME:mm.CustomNonbondedForce.CutoffPeriodic,
ff.LJPME:mm.CustomNonbondedForce.CutoffPeriodic}
lj.setNonbondedMethod(methodMap[nonbondedMethod])
lj.setCutoffDistance(nonbondedCutoff)
# Adjust masses. # Adjust masses.
......
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