Skip to content

GitLab

Unverified Commit 489e2c46 authored by João Morado's avatar João Morado Committed by GitHub
Browse files

Cleanup TinkerReader, ForceField, and amoebaforces (#5080) 

* Fresh branch refactoring the new AMOEBA code

* Finish cleaning up AmoebaAngleForce and AmoebaInPlaneAngleForce

* Cleanup AmoebaTorsionTorsionForce

* Cleanup AmoebaOutOfPlaneBend

* Cleanup AmoebaMultipoleForce

* Remove unnecessary gkForce

* Simplify usage of atomClasses in ForceField

* Formatting

* Fix type of class on WCA

* Simplify angle forces

* Add parsing of units to PiTorsion and StretchBond, and misc. formatting adjustments

* Update code per review feedback

* Clearly defined API for TorsionTorsion,  and correct matching for UB,

* Unindent break statements

* Raise ValueError if classes and types are mixed in a Urey-Bradley term definition
parent 064cee2f
Hide whitespace changes
Inline Side-by-side
Showing with 998 additions and 1515 deletions
wrappers/python/openmm/app/forcefield.py
View file @ 489e2c46
......@@ -637,6 +637,7 @@ class ForceField(object):
self.isAngleConstrained = []
self.constraints = {}
self.bondedToAtom = [set() for _ in self.atoms]
self.atomClasses = []
# Record which atoms are bonded to each other atom
......@@ -669,6 +670,18 @@ class ForceField(object):
if match == site.index:
self.virtualSites[atom] = (site, [matchAtoms[i].index for i in site.atoms], matchAtoms[site.excludeWith].index)
def setAtomClasses(self, forcefield):
"""
Set atom classes for all atoms in the system.
Parameters
----------
forcefield : ForceField
The ForceField object containing atom type definitions
"""
classNameForType = dict((t.name, str(t.atomClass)) for t in forcefield._atomTypes.values())
self.atomClasses = [classNameForType[self.atomType[atom]] for atom in self.atoms]
class _TemplateData(object):
"""Inner class used to encapsulate data about a residue template definition."""
def __init__(self, name):
......@@ -1260,6 +1273,11 @@ class ForceField(object):
elif rigidWater:
rigidResidue[res.index] = True
# Calculate atom classes once for reuse by all generators
# Important for Amoeba forces
data.setAtomClasses(self)
# Create the System and add atoms
sys = mm.System()
......@@ -3531,7 +3549,6 @@ def getAtomPrint(data, atomIndex):
return returnString
#=============================================================================================
def countConstraint(data):
......@@ -3550,64 +3567,35 @@ def countConstraint(data):
## @private
class AmoebaBondGenerator(object):
#=============================================================================================
"""An AmoebaBondGenerator constructs a AmoebaBondForce."""
#=============================================================================================
def __init__(self, cubic, quartic):
self.cubic = cubic
self.quartic = quartic
self.types1 = []
self.types2 = []
self.length = []
self.k = []
#=============================================================================================
self.builder = amoebaforces.AmoebaBondForceBuilder(cubic, quartic)
@staticmethod
def parseElement(element, forceField):
# <AmoebaBondForce bond-cubic="-25.5" bond-quartic="379.3125">
# <Bond class1="1" class2="2" length="0.1437" k="156900.0"/>
generator = AmoebaBondGenerator(element.attrib['bond-cubic'], element.attrib['bond-quartic'])
forceField._forces.append(generator)
for bond in element.findall('Bond'):
types = forceField._findAtomTypes(bond.attrib, 2)
if None not in types:
generator.types1.append(types[0])
generator.types2.append(types[1])
generator.length.append(float(bond.attrib['length']))
generator.k.append(float(bond.attrib['k']))
else:
outputString = "AmoebaBondGenerator: error getting types: %s %s" % (
bond.attrib['class1'],
bond.attrib['class2'])
try:
generator.builder.registerParams((bond.attrib['class1'], bond.attrib['class2']), float(bond.attrib['length']), float(bond.attrib['k']))
except:
outputString = "AmoebaBondGenerator: error getting classes: %s %s" % (bond.attrib['class1'], bond.attrib['class2'])
raise ValueError(outputString)
#=============================================================================================
def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
builder = amoebaforces.AmoebaBondForceBuilder(self.cubic, self.quartic)
force = builder.getForce(sys)
force = self.builder.getForce(sys)
bondsConstraints = []
for bond in data.bonds:
type1 = data.atomType[data.atoms[bond.atom1]]
type2 = data.atomType[data.atoms[bond.atom2]]
for i in range(len(self.types1)):
types1 = self.types1[i]
types2 = self.types2[i]
if (type1 in types1 and type2 in types2) or (type1 in types2 and type2 in types1):
bond.length = self.length[i]
if bond.isConstrained:
data.addConstraint(sys, bond.atom1, bond.atom2, self.length[i])
if self.k[i] != 0:
if not bond.isConstrained or args.get('flexibleConstraints', False):
force.addBond(bond.atom1, bond.atom2, [self.length[i], self.k[i]])
break
bondsConstraints.append(bond.isConstrained)
bondType = (data.atomClasses[bond.atom1], data.atomClasses[bond.atom2])
params = self.builder._findMatchingParams(self.builder.bondParams, bondType)
bond.length = params[0]
if bond.isConstrained:
data.addConstraint(sys, bond.atom1, bond.atom2, params[0])
self.builder.addBonds(force, data.atomClasses, _findBondsForExclusions(data, sys), bondsConstraints, args.get('flexibleConstraints', False))
parsers["AmoebaBondForce"] = AmoebaBondGenerator.parseElement
......@@ -3639,38 +3627,20 @@ def addAngleConstraint(angle, idealAngle, data, sys):
data.addConstraint(sys, angle[0], angle[2], length)
return
#=============================================================================================
## @private
class AmoebaAngleGenerator(object):
#=============================================================================================
"""An AmoebaAngleGenerator constructs a AmoebaAngleForce."""
#=============================================================================================
def __init__(self, forceField, cubic, quartic, pentic, sextic):
self.angleBuilder = amoebaforces.AmoebaAngleForceBuilder(float(cubic), float(quartic), float(pentic), float(sextic))
self.inPlaneAngleBuilder = amoebaforces.AmoebaInPlaneAngleForceBuilder(float(cubic), float(quartic), float(pentic), float(sextic))
self.forceField = forceField
self.cubic = cubic
self.quartic = quartic
self.pentic = pentic
self.sextic = sextic
self.types1 = []
self.types2 = []
self.types3 = []
self.angle = []
self.k = []
self.inPlane = []
#=============================================================================================
self.angleParams = {}
@staticmethod
def parseElement(element, forceField):
# <AmoebaAngleForce angle-cubic="-0.014" angle-quartic="5.6e-05" angle-pentic="-7e-07" angle-sextic="2.2e-08">
# <Angle class1="2" class2="1" class3="3" k="0.0637259642196" angle1="122.00" />
existing = [f for f in forceField._forces if isinstance(f, AmoebaAngleGenerator)]
if len(existing) == 0:
generator = AmoebaAngleGenerator(forceField, element.attrib['angle-cubic'], element.attrib['angle-quartic'], element.attrib['angle-pentic'], element.attrib['angle-sextic'])
......@@ -3680,29 +3650,19 @@ class AmoebaAngleGenerator(object):
if tuple(element.attrib[x] for x in ('angle-cubic', 'angle-quartic', 'angle-pentic', 'angle-sextic')) != (generator.cubic, generator.quartic, generator.pentic, generator.sextic):
raise ValueError('All <AmoebaAngleForce> tags must use identical scale factors')
for angle in element.findall('Angle'):
types = forceField._findAtomTypes(angle.attrib, 3)
if None not in types:
generator.types1.append(types[0])
generator.types2.append(types[1])
generator.types3.append(types[2])
angleList = []
angleList.append(float(angle.attrib['angle1']))
try:
theta0 = [float(angle.attrib['angle1'])]
if 'angle2' in angle.attrib:
angleList.append(float(angle.attrib['angle2']))
theta0.append(float(angle.attrib['angle2']))
if 'angle3' in angle.attrib:
angleList.append(float(angle.attrib['angle3']))
generator.angle.append(angleList)
generator.k.append(float(angle.attrib['k']))
if 'inPlane' in angle.attrib:
generator.inPlane.append(angle.attrib['inPlane'].lower() == 'true')
else:
generator.inPlane.append(None) # for backward compatibility with older versions of AMOEBA
else:
outputString = "AmoebaAngleGenerator: error getting types: %s %s %s" % (
theta0.append(float(angle.attrib['angle3']))
generator.angleParams[(angle.attrib['class1'],
angle.attrib['class2'],
angle.attrib['class3'])] = {"k": float(angle.attrib['k']),
"theta0": theta0,
"inPlane": angle.attrib.get('inPlane', None)}
except Exception as e:
outputString = "AmoebaAngleGenerator: error getting classes: %s %s %s" % (
angle.attrib['class1'],
angle.attrib['class2'],
angle.attrib['class3'])
......@@ -3722,87 +3682,61 @@ class AmoebaAngleGenerator(object):
# non-in-plane angles
#=============================================================================================
def createForcePostOpBendAngle(self, sys, data, nonbondedMethod, nonbondedCutoff, angleList, args):
builder = amoebaforces.AmoebaAngleForceBuilder(self.cubic, self.quartic, self.pentic, self.sextic)
force = builder.getForce(sys)
DEG_TO_RAD = math.pi / 180
for angleDict in angleList:
angle = angleDict['angle']
isConstrained = angleDict['isConstrained']
inPlane = angleDict['inPlane']
type1 = data.atomType[data.atoms[angle[0]]]
type2 = data.atomType[data.atoms[angle[1]]]
type3 = data.atomType[data.atoms[angle[2]]]
for i in range(len(self.types1)):
# self.inPlane is used for modern force fields. inPlane is used for legacy ones that don't specify it.
if self.inPlane[i] or (self.inPlane[i] is None and inPlane):
continue
types1 = self.types1[i]
types2 = self.types2[i]
types3 = self.types3[i]
if (type1 in types1 and type2 in types2 and type3 in types3) or (type1 in types3 and type2 in types2 and type3 in types1):
if isConstrained and self.k[i] != 0.0:
angleDict['idealAngle'] = self.angle[i][0]
addAngleConstraint(angle, self.angle[i][0]*DEG_TO_RAD, data, sys)
if self.k[i] != 0 and (not isConstrained or args.get('flexibleConstraints', False)):
lenAngle = len(self.angle[i])
if (lenAngle > 1):
# get k-index by counting number of non-angle hydrogens on the central atom
# based on kangle.f
numberOfHydrogens = 0
for bond in data.atomBonds[angle[1]]:
atom1 = data.bonds[bond].atom1
atom2 = data.bonds[bond].atom2
if (atom1 == angle[1] and atom2 != angle[0] and atom2 != angle[2] and (sys.getParticleMass(atom2)/unit.dalton) < 1.90):
numberOfHydrogens += 1
if (atom2 == angle[1] and atom1 != angle[0] and atom1 != angle[2] and (sys.getParticleMass(atom1)/unit.dalton) < 1.90):
numberOfHydrogens += 1
if (numberOfHydrogens < lenAngle):
angleValue = self.angle[i][numberOfHydrogens]
else:
outputString = "AmoebaAngleGenerator angle index=%d is out of range: [0, %5d] " % (numberOfHydrogens, lenAngle)
raise ValueError(outputString)
else:
angleValue = self.angle[i][0]
angleDict['idealAngle'] = angleValue
force.addAngle(angle[0], angle[1], angle[2], [angleValue, self.k[i]])
break
def createForcePostOpBendAngle(self, sys, data, nonbondedMethod, nonbondedCutoff, angles, anglesConstraints, args):
idealAngles = {}
for angle, isConstrained in zip(angles, anglesConstraints):
angleClasses = (data.atomClasses[angle[0]], data.atomClasses[angle[1]], data.atomClasses[angle[2]])
params = self.angleParams.get(angleClasses) or self.angleParams.get(angleClasses[::-1])
if params is None:
raise ValueError(
f"AmoebaAngleGenerator: no parameters found for angle: "
f"{angleClasses[0]}-{angleClasses[2]}-"
f"{getAtomPrint(data, angle[0])}, "
f"{getAtomPrint(data, angle[1])}, "
f"{getAtomPrint(data, angle[2])}"
)
angleTuple = tuple(angle)
if isConstrained and params["k"] != 0.0:
idealAngle = params["theta0"][0]
addAngleConstraint(angle, idealAngle*math.pi/180.0, data, sys)
else:
idealAngle = self.angleBuilder.getIdealAngle(angleTuple, params["theta0"], data.atoms, data.bondedToAtom)
self.angleBuilder.registerParams(angleTuple, (idealAngle, params["k"]))
idealAngles[angleTuple] = idealAngle*math.pi/180.0
force = self.angleBuilder.getForce(sys)
self.angleBuilder.addAngles(force, angles, anglesConstraints, args.get('flexibleConstraints', False))
return idealAngles
#=============================================================================================
# createForcePostOpBendInPlaneAngle is called by AmoebaOutOfPlaneBendForce with the list of
# in-plane angles
#=============================================================================================
def createForcePostOpBendInPlaneAngle(self, sys, data, nonbondedMethod, nonbondedCutoff, angleList, args):
builder = amoebaforces.AmoebaInPlaneAngleForceBuilder(self.cubic, self.quartic, self.pentic, self.sextic)
force = builder.getForce(sys)
for angleDict in angleList:
angle = angleDict['angle']
isConstrained = angleDict['isConstrained']
inPlane = angleDict['inPlane']
type1 = data.atomType[data.atoms[angle[0]]]
type2 = data.atomType[data.atoms[angle[1]]]
type3 = data.atomType[data.atoms[angle[2]]]
for i in range(len(self.types1)):
# self.inPlane is used for modern force fields. inPlane is used for legacy ones that don't specify it.
if self.inPlane[i] == False or (self.inPlane[i] is None and not inPlane):
continue
types1 = self.types1[i]
types2 = self.types2[i]
types3 = self.types3[i]
def createForcePostOpBendInPlaneAngle(self, sys, data, nonbondedMethod, nonbondedCutoff, inPlaneAngles, inPlaneAnglesConstraints, args):
idealAngles = {}
for inPlaneAngle, isConstrained in zip(inPlaneAngles, inPlaneAnglesConstraints):
angleClasses = (data.atomClasses[inPlaneAngle[0]], data.atomClasses[inPlaneAngle[1]], data.atomClasses[inPlaneAngle[2]])
params = self.angleParams.get(angleClasses) or self.angleParams.get(angleClasses[::-1])
if params is None:
raise ValueError(
f"AmoebaAngleGenerator: no parameters found for angle: "
f"{angleClasses[0]}-{angleClasses[2]}-"
f"{getAtomPrint(data, inPlaneAngle[0])}, "
f"{getAtomPrint(data, inPlaneAngle[1])}, "
f"{getAtomPrint(data, inPlaneAngle[2])}"
)
angleTuple = tuple(inPlaneAngle)
if isConstrained and params["k"] != 0.0:
idealAngle = params["theta0"][0]
addAngleConstraint(inPlaneAngle, idealAngle*math.pi/180.0, data, sys)
else:
idealAngle = self.angleBuilder.getIdealAngle(angleTuple, params["theta0"], data.atoms, data.bondedToAtom)
self.inPlaneAngleBuilder.registerParams(angleClasses, (idealAngle, params["k"]))
idealAngles[angleTuple[:3]] = idealAngle*math.pi/180.0
force = self.inPlaneAngleBuilder.getForce(sys)
self.inPlaneAngleBuilder.addInPlaneAngles(force, data.atomClasses, inPlaneAngles, inPlaneAnglesConstraints, args.get('flexibleConstraints', False))
if (type1 in types1 and type2 in types2 and type3 in types3) or (type1 in types3 and type2 in types2 and type3 in types1):
angleDict['idealAngle'] = self.angle[i][0]
if (isConstrained and self.k[i] != 0.0):
addAngleConstraint(angle, self.angle[i][0]*math.pi/180.0, data, sys)
if self.k[i] != 0.0 and (not isConstrained or args.get('flexibleConstraints', False)):
force.addBond((angle[0], angle[1], angle[2], angle[3]), (self.angle[i][0], self.k[i]))
break
return idealAngles
parsers["AmoebaAngleForce"] = AmoebaAngleGenerator.parseElement
......@@ -3813,56 +3747,14 @@ parsers["AmoebaAngleForce"] = AmoebaAngleGenerator.parseElement
#=============================================================================================
## @private
class AmoebaOutOfPlaneBendGenerator(object):
#=============================================================================================
"""An AmoebaOutOfPlaneBendGenerator constructs a AmoebaOutOfPlaneBendForce."""
#=============================================================================================
def __init__(self, forceField, type, cubic, quartic, pentic, sextic):
self.forceField = forceField
self.type = type
self.cubic = cubic
self.quartic = quartic
self.pentic = pentic
self.sextic = sextic
self.types1 = []
self.types2 = []
self.types3 = []
self.types4 = []
self.ks = []
#=============================================================================================
# Local version of findAtomTypes needed since class indices are 0 (i.e., not recognized)
# for types3 and 4
#=============================================================================================
def findAtomTypes(self, forceField, node, num):
"""Parse the attributes on an XML tag to find the set of atom types for each atom it involves."""
types = []
attrib = node.attrib
for i in range(num):
if num == 1:
suffix = ''
else:
suffix = str(i+1)
classAttrib = 'class'+suffix
if classAttrib in attrib:
if attrib[classAttrib] in forceField._atomClasses:
types.append(forceField._atomClasses[attrib[classAttrib]])
else:
types.append(set())
return types
#=============================================================================================
self.builder = amoebaforces.AmoebaOutOfPlaneBendForceBuilder(float(cubic), float(quartic), float(pentic), float(sextic))
@staticmethod
def parseElement(element, forceField):
# <AmoebaOutOfPlaneBendForce type="ALLINGER" opbend-cubic="-0.014" opbend-quartic="5.6e-05" opbend-pentic="-7e-07" opbend-sextic="2.2e-08">
# <Angle class1="2" class2="1" class3="0" class4="0" k="0.0531474541591"/>
# <Angle class1="3" class2="1" class3="0" class4="0" k="0.0898536095496"/>
......@@ -3879,154 +3771,62 @@ class AmoebaOutOfPlaneBendGenerator(object):
raise ValueError('All <AmoebaOutOfPlaneBendForce> tags must use identical scale factors')
for angle in element.findall('Angle'):
if 'class3' in angle.attrib and 'class4' in angle.attrib and angle.attrib['class3'] == '0' and angle.attrib['class4'] == '0':
# This is needed for backward compatibility with old AMOEBA force fields that specified wildcards in a nonstandard way.
angle.attrib['class3'] = ''
angle.attrib['class4'] = ''
types = generator.findAtomTypes(forceField, angle, 4)
if types is not None:
generator.types1.append(types[0])
generator.types2.append(types[1])
generator.types3.append(types[2])
generator.types4.append(types[3])
generator.ks.append(float(angle.attrib['k']))
else:
outputString = "AmoebaOutOfPlaneBendGenerator error getting types: %s %s %s %s." % (
try:
generator.builder.registerParams((angle.attrib['class1'], angle.attrib['class2'], angle.attrib['class3'], angle.attrib['class4']), (float(angle.attrib['k']),))
except:
outputString = "AmoebaOutOfPlaneBendGenerator error getting classes: %s %s %s %s." % (
angle.attrib['class1'], angle.attrib['class2'], angle.attrib['class3'], angle.attrib['class4'])
raise ValueError(outputString)
#=============================================================================================
def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
self._nonbondedMethod = nonbondedMethod
self._nonbondedCutoff = nonbondedCutoff
def postprocessSystem(self, sys, data, args):
# We need to wait until after all bonds have been added so their lengths will be set correctly.
builder = amoebaforces.AmoebaOutOfPlaneBendForceBuilder(self.cubic, self.quartic, self.pentic, self.sextic)
force = builder.getForce(sys)
# this hash is used to ensure the out-of-plane-bend bonds
# are only added once
skipAtoms = dict()
angles = []
def addBond(particles):
types = [data.atomType[data.atoms[p]] for p in particles]
for i in range(len(self.types1)):
if types[1] in self.types2[i] and types[3] in self.types1[i]:
if (types[0] in self.types3[i] and types[2] in self.types4[i]) or (types[2] in self.types3[i] and types[0] in self.types4[i]):
force.addBond(particles, [self.ks[i]])
return
for (angle, isConstrained) in zip(data.angles, data.isAngleConstrained):
middleAtom = angle[1]
middleType = data.atomType[data.atoms[middleAtom]]
middleCovalency = len(data.atomBonds[middleAtom])
# if middle atom has covalency of 3 and
# the types of the middle atom and the partner atom (atom bonded to
# middle atom, but not in angle) match types1 and types2, then
# three out-of-plane bend angles are generated. Three in-plane angle
# are also generated. If the conditions are not satisfied, the angle is marked as 'generic' angle (not a in-plane angle)
if middleCovalency == 3 and middleAtom not in skipAtoms:
partners = []
for bond in data.atomBonds[middleAtom]:
atom1 = data.bonds[bond].atom1
atom2 = data.bonds[bond].atom2
if atom1 != middleAtom:
partner = atom1
else:
partner = atom2
partnerType = data.atomType[data.atoms[partner]]
for i in range(len(self.types1)):
types1 = self.types1[i]
types2 = self.types2[i]
if middleType in types2 and partnerType in types1:
partners.append(partner)
break
if len(partners) == 3:
addBond([partners[0], middleAtom, partners[1], partners[2]])
addBond([partners[2], middleAtom, partners[0], partners[1]])
addBond([partners[1], middleAtom, partners[2], partners[0]])
# skipAtoms is used to ensure angles are only included once
skipAtoms[middleAtom] = set(partners[:3])
# in-plane angle
angleDict = {}
angleList = list(angle[:3])
for atomIndex in partners:
if atomIndex not in angleList:
angleList.append(atomIndex)
angleDict['angle'] = angleList
angleDict['isConstrained'] = 0
angleDict['inPlane'] = True
angles.append(angleDict)
else:
angleDict = {}
angleList = list(angle[:3])
for atomIndex in partners:
if atomIndex not in angleList:
angleList.append(atomIndex)
angleDict['angle'] = angleList
angleDict['isConstrained'] = isConstrained
angleDict['inPlane'] = False
angles.append(angleDict)
elif middleCovalency == 3 and middleAtom in skipAtoms:
angleDict = {}
angleList = list(angle[:3])
for atomIndex in skipAtoms[middleAtom]:
if atomIndex not in angleList:
angleList.append(atomIndex)
angleDict['angle'] = angleList
angleDict['isConstrained'] = isConstrained
angleDict['inPlane'] = True
angles.append(angleDict)
else:
angleDict = {}
angleDict['angle'] = angle
angleDict['isConstrained'] = isConstrained
angleDict['inPlane'] = False
angles.append(angleDict)
# get AmoebaAngleGenerator and add AmoebaAngle and AmoebaInPlaneAngle forces
# Classify angles into in-plane, out-of-plane, and generic angles
opbendTypes = [opbendType for opbendType, _ in self.builder.outOfPlaneBendParams]
inPlaneAngles, outOfPlaneAngles, genericAngles = self.builder.classifyAngles(
data.angles,
data.atomClasses,
data.bondedToAtom,
opbendTypes
)
# Add out-of-plane bend force
force = self.builder.getForce(sys)
self.builder.addOutOfPlaneBends(force, data.atomClasses, outOfPlaneAngles)
# Constraints lists
inPlaneAnglesConstraints = []
genericAnglesConstraints = []
inPlaneAnglesCmp = [tuple(angle)[:3] for angle in inPlaneAngles]
for angle, isConstrained in zip(data.angles, data.isAngleConstrained):
if angle in genericAngles:
genericAnglesConstraints.append(isConstrained)
elif angle in inPlaneAnglesCmp:
inPlaneAnglesConstraints.append(isConstrained)
# Get AmoebaAngleGenerator and add AmoebaAngle and AmoebaInPlaneAngle forces
idealAngles = {}
for force in self.forceField._forces:
if (force.__class__.__name__ == 'AmoebaAngleGenerator'):
force.createForcePostOpBendAngle(sys, data, self._nonbondedMethod, self._nonbondedCutoff, angles, args)
force.createForcePostOpBendInPlaneAngle(sys, data, self._nonbondedMethod, self._nonbondedCutoff, angles, args)
idealAnglesTmp = force.createForcePostOpBendAngle(sys, data, self._nonbondedMethod, self._nonbondedCutoff, genericAngles, genericAnglesConstraints, args)
idealAngles.update(idealAnglesTmp)
idealAnglesTmp = force.createForcePostOpBendInPlaneAngle(sys, data, self._nonbondedMethod, self._nonbondedCutoff, inPlaneAngles, inPlaneAnglesConstraints, args)
idealAngles.update(idealAnglesTmp)
for force in self.forceField._forces:
if (force.__class__.__name__ == 'AmoebaStretchBendGenerator'):
force.createForcePostAmoebaBondForce(sys, data, self._nonbondedMethod, self._nonbondedCutoff, angles, args)
force.createForcePostAmoebaBondForce(sys, data, self._nonbondedMethod, self._nonbondedCutoff, data.angles, idealAngles, args)
parsers["AmoebaOutOfPlaneBendForce"] = AmoebaOutOfPlaneBendGenerator.parseElement
#=============================================================================================
## @private
class AmoebaTorsionGenerator(object):
#=============================================================================================
# TODO: remove once amoeba2013.xml has been updated
"""An AmoebaTorsionGenerator constructs a AmoebaTorsionForce."""
#=============================================================================================
def __init__(self, torsionUnit):
self.torsionUnit = torsionUnit
self.types1 = []
......@@ -4038,11 +3838,8 @@ class AmoebaTorsionGenerator(object):
self.t2 = []
self.t3 = []
#=============================================================================================
@staticmethod
def parseElement(element, forceField):
# <AmoebaTorsionForce torsionUnit="0.5">
# <Torsion class1="3" class2="1" class3="2" class4="3" amp1="0.0" angle1="0.0" amp2="0.0" angle2="3.14159265359" amp3="0.0" angle3="0.0" />
# <Torsion class1="3" class2="1" class3="2" class4="6" amp1="0.0" angle1="0.0" amp2="0.0" angle2="3.14159265359" amp3="-0.263592" angle3="0.0" />
......@@ -4079,15 +3876,13 @@ class AmoebaTorsionGenerator(object):
generator.t3.append(tInfo)
else:
outputString = "AmoebaTorsionGenerator: error getting types: %s %s %s %s" % (
outputString = "AmoebaTorsionGenerator: error getting classes: %s %s %s %s" % (
torsion.attrib['class1'],
torsion.attrib['class2'],
torsion.attrib['class3'],
torsion.attrib['class4'])
raise ValueError(outputString)
#=============================================================================================
def createForce(self, sys, data, nontorsionedMethod, nontorsionedCutoff, args):
builder = amoebaforces.AmoebaTorsionForceBuilder()
force = builder.getForce(sys)
......@@ -4119,126 +3914,61 @@ class AmoebaTorsionGenerator(object):
parsers["AmoebaTorsionForce"] = AmoebaTorsionGenerator.parseElement
#=============================================================================================
## @private
class AmoebaPiTorsionGenerator(object):
#=============================================================================================
"""An AmoebaPiTorsionGenerator constructs a AmoebaPiTorsionForce."""
#=============================================================================================
def __init__(self):
self.types1 = []
self.types2 = []
self.k = []
#=============================================================================================
self.builder = amoebaforces.AmoebaPiTorsionForceBuilder()
@staticmethod
def parseElement(element, forceField):
# <AmoebaPiTorsionForce piTorsionUnit="1.0">
# <PiTorsion class1="1" class2="3" k="28.6604" />
# <PiTorsion class1="3" class2="15" k="28.6604" />
generator = AmoebaPiTorsionGenerator()
forceField._forces.append(generator)
piTorsionUnit = float(element.attrib['piTorsionUnit']) if 'piTorsionUnit' in element.attrib else 1.0
for piTorsion in element.findall('PiTorsion'):
types = forceField._findAtomTypes(piTorsion.attrib, 2)
if None not in types:
generator.types1.append(types[0])
generator.types2.append(types[1])
generator.k.append(float(piTorsion.attrib['k']))
else:
outputString = "AmoebaPiTorsionGenerator: error getting types: %s %s " % (
try:
k=float(piTorsion.attrib['k'])*piTorsionUnit
generator.builder.registerParams((piTorsion.attrib['class1'], piTorsion.attrib['class2']), (k,))
except:
outputString = "AmoebaPiTorsionGenerator: error getting classes: %s %s " % (
piTorsion.attrib['class1'],
piTorsion.attrib['class2'])
raise ValueError(outputString)
#=============================================================================================
def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
builder = amoebaforces.AmoebaPiTorsionForceBuilder()
force = builder.getForce(sys)
for bond1 in data.bonds:
# search for bonds with both atoms in bond having covalency == 3
atom1 = bond1.atom1
atom2 = bond1.atom2
if (len(data.atomBonds[atom1]) == 3 and len(data.atomBonds[atom2]) == 3):
type1 = data.atomType[data.atoms[atom1]]
type2 = data.atomType[data.atoms[atom2]]
for i in range(len(self.types1)):
types1 = self.types1[i]
types2 = self.types2[i]
if (type1 in types1 and type2 in types2) or (type1 in types2 and type2 in types1):
# piTorsionAtom1, piTorsionAtom2 are the atoms bonded to atom1, excluding atom2
# piTorsionAtom5, piTorsionAtom6 are the atoms bonded to atom2, excluding atom1
piTorsionAtom1 = -1
piTorsionAtom2 = -1
piTorsionAtom3 = atom1
piTorsionAtom4 = atom2
piTorsionAtom5 = -1
piTorsionAtom6 = -1
for bond in data.atomBonds[atom1]:
bondedAtom1 = data.bonds[bond].atom1
bondedAtom2 = data.bonds[bond].atom2
if (bondedAtom1 != atom1):
b1 = bondedAtom1
else:
b1 = bondedAtom2
if (b1 != atom2):
if (piTorsionAtom1 == -1):
piTorsionAtom1 = b1
else:
piTorsionAtom2 = b1
for bond in data.atomBonds[atom2]:
bondedAtom1 = data.bonds[bond].atom1
bondedAtom2 = data.bonds[bond].atom2
if (bondedAtom1 != atom2):
b1 = bondedAtom1
else:
b1 = bondedAtom2
if (b1 != atom1):
if (piTorsionAtom5 == -1):
piTorsionAtom5 = b1
else:
piTorsionAtom6 = b1
force.addBond([piTorsionAtom1, piTorsionAtom2, piTorsionAtom3, piTorsionAtom4, piTorsionAtom5, piTorsionAtom6], [self.k[i]])
processedPiTorsions = self.builder.getAllPiTorsions(data.atomClasses, data.bondedToAtom, _findBondsForExclusions(data, sys))
force = self.builder.getForce(sys)
self.builder.addPiTorsions(force, data.atomClasses, processedPiTorsions)
parsers["AmoebaPiTorsionForce"] = AmoebaPiTorsionGenerator.parseElement
#=============================================================================================
## @private
class AmoebaStretchTorsionGenerator(object):
"""An AmoebaStretchTorsionGenerator constructs a AmoebaStretchTorsionForce."""
def __init__(self):
self.torsions = []
self.builder = amoebaforces.AmoebaStretchTorsionForceBuilder()
@staticmethod
def parseElement(element, forceField):
# <Torsion class1="44" class2="46" class3="68" class4="65" v11="0.0" v12="0.0" v13="62.760000000000005" v21="0.0" v22="0.0" v23="-167.36" v31="0.0" v32="0.0" v33="217.568"/>
generator = AmoebaStretchTorsionGenerator()
forceField._forces.append(generator)
params = ('v11', 'v12', 'v13', 'v21', 'v22', 'v23', 'v31', 'v32', 'v33')
for torsion in element.findall('Torsion'):
types = forceField._findAtomTypes(torsion.attrib, 4)
if None not in types:
v = [float(torsion.attrib[param]) for param in params]
generator.torsions.append((types, v))
generator.classNameForType = dict((t.name, int(t.atomClass)) for t in forceField._atomTypes.values())
try:
params = tuple(float(torsion.attrib[p]) for p in ('v11', 'v12', 'v13', 'v21', 'v22', 'v23', 'v31', 'v32', 'v33'))
generator.builder.registerParams((torsion.attrib['class1'], torsion.attrib['class2'], torsion.attrib['class3'], torsion.attrib['class4']), params)
except Exception as e:
outputString = "AmoebaStretchTorsionGenerator: error getting classes: %s %s %s %s" % (
torsion.attrib['class1'],
torsion.attrib['class2'],
torsion.attrib['class3'],
torsion.attrib['class4'])
raise ValueError(outputString)
def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
pass
......@@ -4246,67 +3976,35 @@ class AmoebaStretchTorsionGenerator(object):
def postprocessSystem(self, sys, data, args):
# We need to wait until after all bonds and torsions have been added before adding the stretch-torsions,
# since it needs parameters from them.
builder = amoebaforces.AmoebaStretchTorsionForceBuilder()
force = builder.getForce(sys)
# Record parameters for bonds and torsions so we can look them up quickly.
bondForce = [f for f in sys.getForces() if type(f) == mm.CustomBondForce and f.getName() == 'AmoebaBond'][0]
torsionForce = [f for f in sys.getForces() if type(f) == mm.PeriodicTorsionForce][0]
bondLength = {}
torsionPhase = defaultdict(lambda: [0.0, math.pi, 0.0])
for i in range(bondForce.getNumBonds()):
p1, p2, params = bondForce.getBondParameters(i)
bondLength[(p1, p2)] = params[0]
bondLength[(p2, p1)] = params[0]
for i in range(torsionForce.getNumTorsions()):
p1, p2, p3, p4, periodicity, phase, k = torsionForce.getTorsionParameters(i)
if periodicity < 4:
phase = phase.value_in_unit(unit.radian)
torsionPhase[(p1, p2, p3, p4)][periodicity-1] = phase
torsionPhase[(p4, p3, p2, p1)][periodicity-1] = phase
# Add stretch-torsions.
for torsion in data.propers:
type1 = data.atomType[data.atoms[torsion[0]]]
type2 = data.atomType[data.atoms[torsion[1]]]
type3 = data.atomType[data.atoms[torsion[2]]]
type4 = data.atomType[data.atoms[torsion[3]]]
for types, v in self.torsions:
if (type1 in types[3] and type2 in types[2] and type3 in types[1] and type4 in types[0]):
type1, type2, type3, type4 = type4, type3, type2, type1
torsion = tuple(reversed(torsion))
if (type1 in types[0] and type2 in types[1] and type3 in types[2] and type4 in types[3]):
params = list(v)
params.append(bondLength[(torsion[0], torsion[1])])
params.append(bondLength[(torsion[1], torsion[2])])
params.append(bondLength[(torsion[2], torsion[3])])
params += torsionPhase[torsion]
force.addBond(torsion, params)
break
force = self.builder.getForce(sys)
self.builder.addStretchTorsions(sys, force, data.atomClasses, data.propers)
parsers["AmoebaStretchTorsionForce"] = AmoebaStretchTorsionGenerator.parseElement
#=============================================================================================
## @private
class AmoebaAngleTorsionGenerator(object):
"""An AmoebaAngleTorsionGenerator constructs a AmoebaAngleTorsionForce."""
def __init__(self):
self.torsions = []
self.builder = amoebaforces.AmoebaAngleTorsionForceBuilder()
@staticmethod
def parseElement(element, forceField):
# <AmoebaAngleTorsionForce>
# <Torsion class1="44" class2="46" class3="68" class4="65" v11="3.3555680000000003" v12="0.0" v13="-13.903432" v21="0.0" v22="0.0" v23="-2.63592"/>
generator = AmoebaAngleTorsionGenerator()
forceField._forces.append(generator)
params = ('v11', 'v12', 'v13', 'v21', 'v22', 'v23')
for torsion in element.findall('Torsion'):
types = forceField._findAtomTypes(torsion.attrib, 4)
if None not in types:
v = [float(torsion.attrib[param]) for param in params]
generator.torsions.append((types, v))
try:
params = tuple(float(torsion.attrib[p]) for p in ('v11', 'v12', 'v13', 'v21', 'v22', 'v23'))
generator.builder.registerParams((torsion.attrib['class1'], torsion.attrib['class2'], torsion.attrib['class3'], torsion.attrib['class4']), params)
except:
outputString = "AmoebaAngleTorsionGenerator: error getting classes: %s %s %s %s" % (
torsion.attrib['class1'],
torsion.attrib['class2'],
torsion.attrib['class3'],
torsion.attrib['class4'])
raise ValueError(outputString)
def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
pass
......@@ -4314,103 +4012,35 @@ class AmoebaAngleTorsionGenerator(object):
def postprocessSystem(self, sys, data, args):
# We need to wait until after all angles and torsions have been added before adding the angle-torsions,
# since it needs parameters from them.
builder = amoebaforces.AmoebaAngleTorsionForceBuilder()
force = builder.getForce(sys)
# Record parameters for angles and torsions so we can look them up quickly.
angleForce = [f for f in sys.getForces() if type(f) == mm.CustomAngleForce and f.getName() == 'AmoebaAngle'][0]
inPlaneAngleForce = [f for f in sys.getForces() if type(f) == mm.CustomCompoundBondForce and f.getName() == 'AmoebaInPlaneAngle'][0]
torsionForce = [f for f in sys.getForces() if type(f) == mm.PeriodicTorsionForce][0]
equilAngle = {}
torsionPhase = defaultdict(lambda: [0.0, math.pi, 0.0])
angleScale = math.pi/180
for i in range(angleForce.getNumAngles()):
p1, p2, p3, params = angleForce.getAngleParameters(i)
equilAngle[(p1, p2, p3)] = params[0]*angleScale
equilAngle[(p3, p2, p1)] = params[0]*angleScale
for i in range(inPlaneAngleForce.getNumBonds()):
particles, params = inPlaneAngleForce.getBondParameters(i)
equilAngle[tuple(particles[:3])] = params[0]*angleScale
equilAngle[tuple(reversed(particles[:3]))] = params[0]*angleScale
for i in range(torsionForce.getNumTorsions()):
p1, p2, p3, p4, periodicity, phase, k = torsionForce.getTorsionParameters(i)
if periodicity < 4:
phase = phase.value_in_unit(unit.radian)
torsionPhase[(p1, p2, p3, p4)][periodicity-1] = phase
torsionPhase[(p4, p3, p2, p1)][periodicity-1] = phase
# Add stretch-torsions.
for torsion in data.propers:
type1 = data.atomType[data.atoms[torsion[0]]]
type2 = data.atomType[data.atoms[torsion[1]]]
type3 = data.atomType[data.atoms[torsion[2]]]
type4 = data.atomType[data.atoms[torsion[3]]]
for types, v in self.torsions:
if (type1 in types[3] and type2 in types[2] and type3 in types[1] and type4 in types[0]):
type1, type2, type3, type4 = type4, type3, type2, type1
torsion = tuple(reversed(torsion))
if (type1 in types[0] and type2 in types[1] and type3 in types[2] and type4 in types[3]):
params = list(v)
params.append(equilAngle[(torsion[0], torsion[1], torsion[2])])
params.append(equilAngle[(torsion[1], torsion[2], torsion[3])])
params += torsionPhase[torsion]
force.addBond(torsion, params)
break
force = self.builder.getForce(sys)
self.builder.addAngleTorsions(sys, force, data.atomClasses, data.propers)
parsers["AmoebaAngleTorsionForce"] = AmoebaAngleTorsionGenerator.parseElement
#=============================================================================================
## @private
class AmoebaTorsionTorsionGenerator(object):
#=============================================================================================
"""An AmoebaTorsionTorsionGenerator constructs a AmoebaTorsionTorsionForce."""
#=============================================================================================
def __init__(self):
self.types1 = []
self.types2 = []
self.types3 = []
self.types4 = []
self.types5 = []
self.gridIndex = []
self.grids = []
#=============================================================================================
self.builder = amoebaforces.AmoebaTorsionTorsionForceBuilder()
@staticmethod
def parseElement(element, forceField):
generator = AmoebaTorsionTorsionGenerator()
forceField._forces.append(generator)
maxGridIndex = -1
# <AmoebaTorsionTorsionForce >
# <TorsionTorsion class1="3" class2="1" class3="2" class4="3" class5="1" grid="0" nx="25" ny="25" />
generator = AmoebaTorsionTorsionGenerator()
forceField._forces.append(generator)
# Register torsion-torsion parameters
for torsionTorsion in element.findall('TorsionTorsion'):
types = forceField._findAtomTypes(torsionTorsion.attrib, 5)
if None not in types:
generator.types1.append(types[0])
generator.types2.append(types[1])
generator.types3.append(types[2])
generator.types4.append(types[3])
generator.types5.append(types[4])
try:
gridIndex = int(torsionTorsion.attrib['grid'])
if (gridIndex > maxGridIndex):
maxGridIndex = gridIndex
generator.gridIndex.append(gridIndex)
else:
outputString = "AmoebaTorsionTorsionGenerator: error getting types: %s %s %s %s %s" % (
torsionTorsionType = (torsionTorsion.attrib['class1'], torsionTorsion.attrib['class2'],
torsionTorsion.attrib['class3'], torsionTorsion.attrib['class4'],
torsionTorsion.attrib['class5'])
generator.builder.registerParams(torsionTorsionType, gridIndex)
except:
outputString = "AmoebaTorsionTorsionGenerator: error getting classes: %s %s %s %s %s" % (
torsionTorsion.attrib['class1'],
torsionTorsion.attrib['class2'],
torsionTorsion.attrib['class3'],
......@@ -4418,38 +4048,28 @@ class AmoebaTorsionTorsionGenerator(object):
torsionTorsion.attrib['class5'] )
raise ValueError(outputString)
# load grid
# xml source
# Register grid data
# <TorsionTorsionGrid grid="0" nx="25" ny="25" >
# <Grid angle1="-180.00" angle2="-180.00" f="0.0" fx="2.31064374824e-05" fy="0.0" fxy="-0.0052801799672" />
# <Grid angle1="-165.00" angle2="-180.00" f="-0.66600912" fx="-0.06983370052" fy="-0.075058725744" fxy="-0.0044462732032" />
# output grid:
# grid[x][y][0] = x value
# grid[x][y][1] = y value
# grid[x][y][2] = function value
# grid[x][y][3] = dfdx value
# grid[x][y][4] = dfdy value
# grid[x][y][5] = dfd(xy) value
# Grid format:
# grid[x][y][0] = x value (angle1)
# grid[x][y][1] = y value (angle2)
# grid[x][y][2] = function value (f)
# grid[x][y][3] = dfdx value (fx)
# grid[x][y][4] = dfdy value (fy)
# grid[x][y][5] = dfd(xy) value (fxy)
maxGridIndex += 1
generator.grids = maxGridIndex*[]
for torsionTorsionGrid in element.findall('TorsionTorsionGrid'):
gridIndex = int(torsionTorsionGrid.attrib[ "grid"])
nx = int(torsionTorsionGrid.attrib[ "nx"])
ny = int(torsionTorsionGrid.attrib[ "ny"])
gridIndex = int(torsionTorsionGrid.attrib["grid"])
nx = int(torsionTorsionGrid.attrib["nx"])
ny = int(torsionTorsionGrid.attrib["ny"])
grid = []
gridCol = []
gridColIndex = 0
for gridEntry in torsionTorsionGrid.findall('Grid'):
gridRow = []
gridRow.append(float(gridEntry.attrib['angle1']))
gridRow.append(float(gridEntry.attrib['angle2']))
......@@ -4460,78 +4080,47 @@ class AmoebaTorsionTorsionGenerator(object):
gridRow.append(float(gridEntry.attrib['fxy']))
gridCol.append(gridRow)
gridColIndex += 1
if (gridColIndex == nx):
gridColIndex += 1
if gridColIndex == nx:
grid.append(gridCol)
gridCol = []
gridColIndex = 0
if (gridIndex == len(generator.grids)):
generator.grids.append(grid)
else:
while(len(generator.grids) < gridIndex):
generator.grids.append([])
generator.grids[gridIndex] = grid
#=============================================================================================
generator.builder.registerGridData(gridIndex, grid)
def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
builder = amoebaforces.AmoebaTorsionTorsionForceBuilder()
force = builder.getForce(sys)
# Add torsion-torsion interactions
builder.addTorsionTorsionInteractions(force, data, self.types1, self.types2, self.types3,
self.types4, self.types5, self.gridIndex, sys)
# Set grids
for (index, grid) in enumerate(self.grids):
builder.setTorsionTorsionGrid(force, index, grid)
force = self.builder.getForce(sys)
masses = [sys.getParticleMass(i).value_in_unit(unit.amu) for i in range(sys.getNumParticles())]
self.builder.addTorsionTorsionInteractions(force, data.propers, data.atomClasses, masses)
parsers["AmoebaTorsionTorsionForce"] = AmoebaTorsionTorsionGenerator.parseElement
#=============================================================================================
## @private
class AmoebaStretchBendGenerator(object):
#=============================================================================================
"""An AmoebaStretchBendGenerator constructs a AmoebaStretchBendForce."""
#=============================================================================================
def __init__(self, forcefield):
self.forcefield = forcefield
self.types1 = []
self.types2 = []
self.types3 = []
self.k1 = []
self.k2 = []
#=============================================================================================
self.builder = amoebaforces.AmoebaStretchBendForceBuilder()
self.stretchBendParams = {}
@staticmethod
def parseElement(element, forceField):
generator = AmoebaStretchBendGenerator(forceField)
forceField._forces.append(generator)
# <AmoebaStretchBendForce stretchBendUnit="1.0">
# <StretchBend class1="2" class2="1" class3="3" k1="5.25776946506" k2="5.25776946506" />
# <StretchBend class1="2" class2="1" class3="4" k1="3.14005676385" k2="3.14005676385" />
generator = AmoebaStretchBendGenerator(forceField)
forceField._forces.append(generator)
stretchBendUnit = float(element.attrib['stretchBendUnit']) if 'stretchBendUnit' in element.attrib else 1.0
for stretchBend in element.findall('StretchBend'):
types = forceField._findAtomTypes(stretchBend.attrib, 3)
if None not in types:
generator.types1.append(types[0])
generator.types2.append(types[1])
generator.types3.append(types[2])
generator.k1.append(float(stretchBend.attrib['k1']))
generator.k2.append(float(stretchBend.attrib['k2']))
else:
outputString = "AmoebaStretchBendGenerator : error getting types: %s %s %s" % (
try:
class1 = stretchBend.attrib['class1']
class2 = stretchBend.attrib['class2']
class3 = stretchBend.attrib['class3']
k1 = float(stretchBend.attrib['k1'])*stretchBendUnit
k2 = float(stretchBend.attrib['k2'])*stretchBendUnit
generator.stretchBendParams[(class1, class2, class3)] = {'k1': k1, 'k2': k2}
except:
outputString = "AmoebaStretchBendGenerator : error getting classes: %s %s %s" % (
stretchBend.attrib['class1'],
stretchBend.attrib['class2'],
stretchBend.attrib['class3'])
......@@ -4558,108 +4147,42 @@ class AmoebaStretchBendGenerator(object):
#=============================================================================================
def createForcePostAmoebaBondForce(self, sys, data, nonbondedMethod, nonbondedCutoff, angleList, args):
builder = amoebaforces.AmoebaStretchBendForceBuilder()
force = builder.getForce(sys)
for angleDict in angleList:
angle = angleDict['angle']
type1 = data.atomType[data.atoms[angle[0]]]
type2 = data.atomType[data.atoms[angle[1]]]
type3 = data.atomType[data.atoms[angle[2]]]
radian = 57.2957795130
for i in range(len(self.types1)):
types1 = self.types1[i]
types2 = self.types2[i]
types3 = self.types3[i]
# match types
# get ideal bond lengths, bondAB, bondCB
# get ideal angle
if (type2 in types2 and ((type1 in types1 and type3 in types3) or (type3 in types1 and type1 in types3))):
bondAB = -1.0
bondCB = -1.0
swap = 0
for bond in data.atomBonds[angle[1]]:
atom1 = data.bonds[bond].atom1
atom2 = data.bonds[bond].atom2
length = data.bonds[bond].length
if (atom1 == angle[0]):
bondAB = length
if (atom1 == angle[2]):
bondCB = length
if (atom2 == angle[2]):
bondCB = length
if (atom2 == angle[0]):
bondAB = length
# check that ideal angle and bonds are set
if ('idealAngle' not in angleDict):
outputString = "AmoebaStretchBendGenerator: ideal angle is not set for following entry:\n"
outputString += " types: %5s %5s %5s atoms: " % (type1, type2, type3)
outputString += getAtomPrint( data, angle[0] ) + ' '
outputString += getAtomPrint( data, angle[1] ) + ' '
outputString += getAtomPrint( data, angle[2] )
raise ValueError(outputString)
elif (bondAB < 0 or bondCB < 0):
outputString = "AmoebaStretchBendGenerator: bonds not set: %15.7e %15.7e. for following entry:" % (bondAB, bondCB)
outputString += " types: [%5s %5s %5s] atoms: " % (type1, type2, type3)
outputString += getAtomPrint( data, angle[0] ) + ' '
outputString += getAtomPrint( data, angle[1] ) + ' '
outputString += getAtomPrint( data, angle[2] )
raise ValueError(outputString)
else:
if type1 in types1 and type3 in types3:
k1, k2 = self.k1[i], self.k2[i]
else:
k1, k2 = self.k2[i], self.k1[i]
force.addBond((angle[0], angle[1], angle[2]), (bondAB, bondCB, angleDict['idealAngle']/radian, k1, k2))
break
def createForcePostAmoebaBondForce(self, sys, data, nonbondedMethod, nonbondedCutoff, angles, idealAngles, args):
# Get ideal bond lengths from AmoebaBondForce
bondParams = {(data.atomClasses[bond.atom1], data.atomClasses[bond.atom2]): {"r0": bond.length} for bond in data.bonds}
processedAngles = self.builder.registerAllStretchBendParams(data.atomClasses, angles, self.stretchBendParams, bondParams, idealAngles)
force = self.builder.getForce(sys)
self.builder.addStretchBends(force, processedAngles)
parsers["AmoebaStretchBendForce"] = AmoebaStretchBendGenerator.parseElement
#=============================================================================================
## @private
class AmoebaVdwGenerator(object):
"""A AmoebaVdwGenerator constructs a AmoebaVdwForce."""
#=============================================================================================
"""An AmoebaVdwGenerator constructs a AmoebaVdwForce."""
def __init__(self, type, radiusrule, radiustype, radiussize, epsilonrule, vdw13Scale, vdw14Scale, vdw15Scale):
self.type = type
self.radiusrule = radiusrule
self.radiustype = radiustype
self.radiussize = radiussize
self.epsilonrule = epsilonrule
self.vdw13Scale = vdw13Scale
self.vdw14Scale = vdw14Scale
self.vdw15Scale = vdw15Scale
self.params = {}
self.pairs = []
#=============================================================================================
self.builder = amoebaforces.AmoebaVdwForceBuilder(str(type),
str(radiusrule),
str(radiustype),
str(radiussize),
str(epsilonrule),
float(vdw13Scale),
float(vdw14Scale),
float(vdw15Scale))
@staticmethod
def parseElement(element, forceField):
# <AmoebaVdwForce type="BUFFERED-14-7" radiusrule="CUBIC-MEAN" radiustype="R-MIN" radiussize="DIAMETER" epsilonrule="HHG" vdw-13-scale="0.0" vdw-14-scale="1.0" vdw-15-scale="1.0" >
# <Vdw class="1" sigma="0.371" epsilon="0.46024000000000004" reduction="1.0"/>
# <Vdw class="2" sigma="0.382" epsilon="0.42258400000000007" reduction="1.0"/>
existing = [f for f in forceField._forces if isinstance(f, AmoebaVdwGenerator)]
if len(existing) == 0:
generator = AmoebaVdwGenerator(element.attrib['type'], element.attrib['radiusrule'], element.attrib['radiustype'], element.attrib['radiussize'], element.attrib['epsilonrule'],
float(element.attrib['vdw-13-scale']), float(element.attrib['vdw-14-scale']), float(element.attrib['vdw-15-scale']))
element.attrib['vdw-13-scale'], element.attrib['vdw-14-scale'], element.attrib['vdw-15-scale'])
forceField.registerGenerator(generator)
else:
# Multiple <AmoebaVdwForce> tags were found, probably in different files. Simply add more types to the existing one.
# Multiple <AmoebaVdwForce> tags were found, probably in different files. Simply add more types to the existing one.
generator = existing[0]
if abs(generator.vdw13Scale - float(element.attrib['vdw-13-scale'])) > NonbondedGenerator.SCALETOL or \
abs(generator.vdw14Scale - float(element.attrib['vdw-14-scale'])) > NonbondedGenerator.SCALETOL or \
......@@ -4669,38 +4192,22 @@ class AmoebaVdwGenerator(object):
generator.radiustype != element.attrib['radiustype'] or generator.radiussize != element.attrib['radiussize']:
raise ValueError('Found multiple AmoebaVdwForce tags with different combining rules')
for vdw in element.findall('Vdw'):
generator.params[vdw.attrib['class']] = tuple(float(vdw.attrib[name]) for name in ('sigma', 'epsilon', 'reduction'))
generator.builder.registerClassParams(vdw.attrib['class'], float(vdw.attrib['sigma']), float(vdw.attrib['epsilon']), float(vdw.attrib['reduction']))
for pair in element.findall('Pair'):
generator.pairs.append((pair.attrib['class1'], pair.attrib['class2'], float(pair.attrib['sigma']), float(pair.attrib['epsilon'])))
generator.classNameForType = dict((t.name, int(t.atomClass)) for t in forceField._atomTypes.values())
#=============================================================================================
generator.builder.registerPairParams(pair.attrib['class1'], pair.attrib['class2'], float(pair.attrib['sigma']), float(pair.attrib['epsilon']))
def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
builder = amoebaforces.AmoebaVdwForceBuilder(self.type, self.radiusrule, self.radiustype, self.radiussize, self.epsilonrule, self.vdw13Scale, self.vdw14Scale, self.vdw15Scale)
for atomClass, params in self.params.items():
builder.registerClassParams(int(atomClass), *params)
for params in self.pairs:
builder.registerPairParams(int(params[0]), int(params[1]), params[2], params[3])
force = builder.getForce(sys, nonbondedMethod, args.get('vdwCutoff', nonbondedCutoff), args.get('useDispersionCorrection', True))
atomClasses = [self.classNameForType[data.atomType[atom]] for atom in data.atoms]
builder.addParticles(force, atomClasses, data.atoms, _findBondsForExclusions(data, sys))
force = self.builder.getForce(sys, nonbondedMethod, args.get('vdwCutoff', nonbondedCutoff), args.get('useDispersionCorrection', True))
self.builder.addParticles(force, data.atomClasses, data.atoms, _findBondsForExclusions(data, sys))
parsers["AmoebaVdwForce"] = AmoebaVdwGenerator.parseElement
#=============================================================================================
## @private
class AmoebaMultipoleGenerator(object):
#=============================================================================================
"""A AmoebaMultipoleGenerator constructs an AmoebaMultipoleForce."""
#=============================================================================================
def __init__(self, forceField):
self.multipoleType = defaultdict(list)
self.polarizationType = {}
self.builder = amoebaforces.AmoebaMultipoleForceBuilder()
@staticmethod
def parseElement(element, forceField):
......@@ -4709,35 +4216,32 @@ class AmoebaMultipoleGenerator(object):
generator = AmoebaMultipoleGenerator(forceField)
forceField.registerGenerator(generator)
else:
# Multiple <AmoebaMultipoleForce> tags were found, probably in different files. Simply add more types to the existing one.
# Multiple <AmoebaMultipoleForce> tags were found, probably in different files. Simply add more types to the existing one.
generator = existing[0]
# set type map: [ kIndices, multipoles, AMOEBA/OpenMM axis type]
# Register multipole parameters using the builder
for atom in element.findall('Multipole'):
types = forceField._findAtomTypes(atom.attrib, 1)
if None not in types:
# k-indices not provided default to 0
kIndices = [int(atom.attrib['type'])]
kStrings = [ 'kz', 'kx', 'ky' ]
for kString in kStrings:
kIndices.append(int(atom.attrib.get(kString,0)))
# set multipole
# Set multipole
charge = float(atom.attrib['c0'])
conversion = 1.0
dipole = [ conversion*float(atom.attrib['d1']), conversion*float(atom.attrib['d2']), conversion*float(atom.attrib['d3'])]
quadrupole = [conversion*float(atom.attrib[a]) for a in ['q11', 'q21', 'q31', 'q21', 'q22', 'q32', 'q31', 'q32', 'q33']]
for t in types[0]:
generator.multipoleType[t].append({'kIndices':kIndices, 'charge':charge, 'dipole':dipole, 'quadrupole':quadrupole})
generator.builder.registerMultipoleParams(int(t),
amoebaforces.MultipoleParams(kIndices, charge, dipole, quadrupole))
else:
outputString = "AmoebaMultipoleGenerator: error getting type for multipole: %s" % (atom.attrib['class'])
outputString = "AmoebaMultipoleGenerator: error getting type for multipole: %s" % (atom.attrib['type'])
raise ValueError(outputString)
# polarization parameters
# Register polarization parameters using the builder
for atom in element.findall('Polarize'):
types = forceField._findAtomTypes(atom.attrib, 1)
if None not in types:
......@@ -4748,57 +4252,45 @@ class AmoebaMultipoleGenerator(object):
pgrp = f'pgrp{index}'
if pgrp in atom.attrib:
groupTypes.append(int(atom.attrib[pgrp]))
for t in types[0]:
if t in generator.polarizationType:
raise ValueError(f'Multipole Polarize tags found for type {t}')
generator.polarizationType[t] = {'polarizability': polarizability, 'thole':thole, 'groupTypes':groupTypes}
if int(t) in generator.builder.polarizationParams:
raise ValueError(f'Multiple Polarize tags found for type {t}')
generator.builder.registerPolarizationParams(int(t),
amoebaforces.PolarizationParams(polarizability, thole, groupTypes))
else:
outputString = "AmoebaMultipoleGenerator: error getting type for polarize: %s" % (atom.attrib['class'])
outputString = "AmoebaMultipoleGenerator: error getting class for polarize: %s" % (atom.attrib.get('class', 'unknown'))
raise ValueError(outputString)
def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
builder = amoebaforces.AmoebaMultipoleForceBuilder()
for type, params in self.multipoleType.items():
for p in params:
builder.registerMultipoleParams(int(type), amoebaforces.MultipoleParams(p['kIndices'], p['charge'], p['dipole'], p['quadrupole']))
for type, params in self.polarizationType.items():
builder.registerPolarizationParams(int(type), amoebaforces.PolarizationParams(params['polarizability'], params['thole'], params['groupTypes']))
force = builder.getForce(sys, nonbondedMethod, nonbondedCutoff, args.get('ewaldErrorTolerance', 1e-4), args.get('polarization', 'mutual'),
args.get('mutualInducedTargetEpsilon', 1e-5), args.get('mutualInducedMaxIterations', 60))
force = self.builder.getForce(sys, nonbondedMethod, nonbondedCutoff,
args.get('ewaldErrorTolerance', 1e-4),
args.get('polarization', 'mutual'),
args.get('mutualInducedTargetEpsilon', 1e-5),
args.get('mutualInducedMaxIterations', 60))
atomTypes = [int(data.atomType[atom]) for atom in data.atoms]
builder.addMultipoles(force, atomTypes, data.atoms, _findBondsForExclusions(data, sys))
self.builder.addMultipoles(force, atomTypes, data.atoms, _findBondsForExclusions(data, sys))
parsers["AmoebaMultipoleForce"] = AmoebaMultipoleGenerator.parseElement
#=============================================================================================
## @private
class AmoebaWcaDispersionGenerator(object):
"""A AmoebaWcaDispersionGenerator constructs a AmoebaWcaDispersionForce."""
#=========================================================================================
def __init__(self, epso, epsh, rmino, rminh, awater, slevy, dispoff, shctd):
self.epso = epso
self.epsh = epsh
self.rmino = rmino
self.rminh = rminh
self.awater = awater
self.slevy = slevy
self.dispoff = dispoff
self.shctd = shctd
self.params = {}
#=========================================================================================
self.builder = amoebaforces.AmoebaWcaDispersionForceBuilder(float(epso),
float(epsh),
float(rmino),
float(rminh),
float(awater),
float(slevy),
float(dispoff),
float(shctd))
@staticmethod
def parseElement(element, forceField):
# <AmoebaWcaDispersionForce epso="0.46024" epsh="0.056484" rmino="0.17025" rminh="0.13275" awater="33.428" slevy="1.0" dispoff="0.026" shctd="0.81" >
# <WcaDispersion class="1" radius="0.1855" epsilon="0.46024" />
# <WcaDispersion class="2" radius="0.191" epsilon="0.422584" />
existing = [f for f in forceField._forces if isinstance(f, AmoebaWcaDispersionGenerator)]
if len(existing) == 0:
generator = AmoebaWcaDispersionGenerator(element.attrib['epso'],
......@@ -4811,60 +4303,36 @@ class AmoebaWcaDispersionGenerator(object):
element.attrib['shctd'])
forceField.registerGenerator(generator)
else:
# Multiple <AmoebaWcaDispersionForce> tags were found, probably in different files. Simply add more types to the existing one.
# Multiple <AmoebaWcaDispersionForce> tags were found, probably in different files. Simply add more types to the existing one.
generator = existing[0]
for wca in element.findall('WcaDispersion'):
generator.params[int(wca.attrib['class'])] = tuple(float(wca.attrib[name]) for name in ('radius', 'epsilon'))
generator.classNameForType = dict((t.name, int(t.atomClass)) for t in forceField._atomTypes.values())
#=========================================================================================
generator.builder.registerParams(wca.attrib['class'], float(wca.attrib['radius']), float(wca.attrib['epsilon']))
def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
builder = amoebaforces.AmoebaWcaDispersionForceBuilder(float(self.epso),
float(self.epsh),
float(self.rmino),
float(self.rminh),
float(self.awater),
float(self.slevy),
float(self.dispoff),
float(self.shctd))
force = builder.getForce(sys)
atomClasses = [self.classNameForType[data.atomType[atom]] for atom in data.atoms]
for atomClass, params in self.params.items():
builder.registerParams(atomClass, *params)
builder.addParticles(force, atomClasses)
force = self.builder.getForce(sys)
self.builder.addParticles(force, data.atomClasses)
parsers["AmoebaWcaDispersionForce"] = AmoebaWcaDispersionGenerator.parseElement
#=============================================================================================
## @private
class AmoebaGeneralizedKirkwoodGenerator(object):
#=========================================================================================
"""A AmoebaGeneralizedKirkwoodGenerator constructs a AmoebaGeneralizedKirkwoodForce."""
#=========================================================================================
def __init__(self, forceField, solventDielectric, soluteDielectric, includeCavityTerm, probeRadius, surfaceAreaFactor):
self.builder = amoebaforces.AmoebaGeneralizedKirkwoodForceBuilder(float(solventDielectric),
float(soluteDielectric),
bool(includeCavityTerm),
float(probeRadius),
float(surfaceAreaFactor))
self.forceField = forceField
self.solventDielectric = solventDielectric
self.soluteDielectric = soluteDielectric
self.includeCavityTerm = includeCavityTerm
self.probeRadius = probeRadius
self.surfaceAreaFactor = surfaceAreaFactor
#=========================================================================================
@staticmethod
def parseElement(element, forceField):
# <AmoebaGeneralizedKirkwoodForce solventDielectric="78.3" soluteDielectric="1.0" includeCavityTerm="1" probeRadius="0.14" surfaceAreaFactor="-170.351730663">
# <GeneralizedKirkwood type="1" charge="-0.22620" shct="0.79" />
# <GeneralizedKirkwood type="2" charge="-0.15245" shct="0.72" />
existing = [f for f in forceField._forces if isinstance(f, AmoebaGeneralizedKirkwoodGenerator)]
if len(existing) == 0:
generator = AmoebaGeneralizedKirkwoodGenerator(forceField, element.attrib['solventDielectric'],
......@@ -4877,110 +4345,81 @@ class AmoebaGeneralizedKirkwoodGenerator(object):
# Multiple <AmoebaGeneralizedKirkwoodFprce> tags were found, probably in different files. Simply add more types to the existing one.
generator = existing[0]
#=========================================================================================
def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
if( nonbondedMethod != NoCutoff ):
raise ValueError( "Only the nonbondedMethod=NoCutoff option is available for implicit solvent simulations." )
# check if AmoebaMultipoleForce exists since charges needed
# if it has not been created, raise an error
amoebaMultipoleForceList = [f for f in sys.getForces() if type(f) == mm.AmoebaMultipoleForce]
if (len(amoebaMultipoleForceList) > 0):
amoebaMultipoleForce = amoebaMultipoleForceList[0]
else:
# call AmoebaMultipoleForceGenerator.createForce() to ensure charges have been set
for force in self.forceField._forces:
if (force.__class__.__name__ == 'AmoebaMultipoleGenerator'):
force.createForce(sys, data, nonbondedMethod, nonbondedCutoff, args)
# Use the amoebaforces builder
solventDielectric = args.get('solventDielectric', self.solventDielectric)
soluteDielectric = args.get('soluteDielectric', self.soluteDielectric)
includeCavityTerm = args.get('includeCavityTerm', self.includeCavityTerm)
builder = amoebaforces.AmoebaGeneralizedKirkwoodForceBuilder(float(solventDielectric),
float(soluteDielectric),
int(includeCavityTerm),
float(self.probeRadius),
float(self.surfaceAreaFactor))
break
# Get the charges from the AmoebaMultipoleForce and register them with the builder
for atomIndex in range(0, amoebaMultipoleForce.getNumMultipoles()):
for atomIndex in range(amoebaMultipoleForce.getNumMultipoles()):
multipoleParameters = amoebaMultipoleForce.getMultipoleParameters(atomIndex)
builder.registerParams(multipoleParameters[0])
self.builder.registerParams(multipoleParameters[0])
force = builder.getForce(sys, implicitSolvent=True)
force = self.builder.getForce(sys, implicitSolvent=True)
if force is not None:
builder.addParticles(force, data.atoms, _findBondsForExclusions(data, sys))
self.builder.addParticles(force, data.atoms, _findBondsForExclusions(data, sys))
parsers["AmoebaGeneralizedKirkwoodForce"] = AmoebaGeneralizedKirkwoodGenerator.parseElement
#=============================================================================================
## @private
class AmoebaUreyBradleyGenerator(object):
#=============================================================================================
"""An AmoebaUreyBradleyGenerator constructs a AmoebaUreyBradleyForce."""
#=============================================================================================
def __init__(self):
self.anglesForAtom2Type = defaultdict(list)
self.types1 = []
self.types2 = []
self.types3 = []
self.length = []
self.k = []
#=============================================================================================
self.builder = amoebaforces.AmoebaUreyBradleyForceBuilder()
@staticmethod
def parseElement(element, forceField):
# <AmoebaUreyBradleyForce>
# <UreyBradley class1="74" class2="73" class3="74" k="16003.8" d="0.15537" />
generator = AmoebaUreyBradleyGenerator()
forceField._forces.append(generator)
for bond in element.findall('UreyBradley'):
types = forceField._findAtomTypes(bond.attrib, 3)
if None not in types:
index = len(generator.types1)
generator.types1.append(types[0])
generator.types2.append(types[1])
generator.types3.append(types[2])
for t in types[1]:
generator.anglesForAtom2Type[t].append(index)
generator.length.append(float(bond.attrib['d']))
generator.k.append(float(bond.attrib['k']))
#=============================================================================================
try:
isClass = any(k.startswith('class') for k in bond.attrib)
isType = any(k.startswith('type') for k in bond.attrib)
if isClass and isType:
raise ValueError(
"AmoebaUreyBradleyGenerator: Cannot mix classes and types in the same Urey-Bradley term definition. "
"Please contact the developers if you require this functionality."
)
prefix = 'class' if isClass else 'type'
keys = (prefix + "1", prefix + "2", prefix + "3")
k = float(bond.attrib['k'])
d = float(bond.attrib['d'])
generator.builder.registerParams(
(bond.attrib[keys[0]], bond.attrib[keys[1]], bond.attrib[keys[2]]),
(k, d),
isClass
)
except Exception as e:
outputString = f"AmoebaUreyBradleyGenerator : error getting {'classes' if isClass else 'types'}: " \
f"{bond.attrib[keys[0]]} {bond.attrib[keys[1]]} {bond.attrib[keys[2]]} " \
f"({e})"
raise ValueError(outputString)
def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
builder = amoebaforces.AmoebaUreyBradleyForceBuilder()
force = builder.getForce(sys)
for (angle, isConstrained) in zip(data.angles, data.isAngleConstrained):
if (isConstrained and not args.get('flexibleConstraints', False)):
continue
type1 = data.atomType[data.atoms[angle[0]]]
type2 = data.atomType[data.atoms[angle[1]]]
type3 = data.atomType[data.atoms[angle[2]]]
for i in self.anglesForAtom2Type[type2]:
types1 = self.types1[i]
types2 = self.types2[i]
types3 = self.types3[i]
if ((type1 in types1 and type2 in types2 and type3 in types3) or (type3 in types1 and type2 in types2 and type1 in types3)):
force.addBond(angle[0], angle[2], self.length[i], 2*self.k[i])
break
force = self.builder.getForce(sys)
# Match against parameters registered by class
self.builder.addUreyBradleys(force, data.atomClasses, data.angles, True, data.isAngleConstrained, args.get('flexibleConstraints', False))
# Match against parameters registered by type
atomTypes = [data.atomType[atom] for atom in data.atoms]
self.builder.addUreyBradleys(force, atomTypes, data.angles, False, data.isAngleConstrained, args.get('flexibleConstraints', False))
parsers["AmoebaUreyBradleyForce"] = AmoebaUreyBradleyGenerator.parseElement
#=============================================================================================
## @private
class HippoNonbondedGenerator(object):
"""A HippoNonbondedGenerator constructs a HippoNonbondedForce."""
......@@ -5178,6 +4617,4 @@ class DrudeGenerator(object):
sys.setParticleMass(particle, drudeMass)
sys.setParticleMass(parent, sys.getParticleMass(parent)-transferMass)
parsers["DrudeForce"] = DrudeGenerator.parseElement
#=============================================================================================
parsers["DrudeForce"] = DrudeGenerator.parseElement
\ No newline at end of file
wrappers/python/openmm/app/internal/amoebaforces.py
View file @ 489e2c46
......@@ -65,13 +65,35 @@ class BaseAmoebaForceBuilder:
ValueError
If the atom does not have an associated element or atomic number.
"""
if hasattr(atom, 'element'):
return atom.element.atomic_number
elif hasattr(atom, 'atomicNumber'):
if hasattr(atom, 'atomicNumber'):
return atom.atomicNumber
elif hasattr(atom, 'element'):
return atom.element.atomic_number
else:
raise ValueError(f"Atom {atom} does not have an associated element or atomic number.")
@staticmethod
def _getNeighbors(atom: Any) -> List[int]:
"""
Get the indices of neighboring atoms bonded to the given atom.
Parameters
----------
atom : Any
Atom object that has either a bonds attribute or a bondedAtoms attribute.
Returns
-------
List[int]
List of indices of neighboring atoms.
"""
if hasattr(atom, 'bondedAtoms'):
return [a.index for a in atom.bondedAtoms]
elif hasattr(atom, 'bonds'):
return [a for a in atom.bonds]
else:
raise ValueError(f"Atom {atom} does not have bonded atoms information.")
def _findExistingForce(self, sys: mm.System, forceType: type, energyFunction: Optional[str] = None, name: Optional[str] = None) -> Optional[Any]:
"""
Find existing force in system by type and optionally by energy function or name.
......@@ -178,37 +200,99 @@ class BaseAmoebaForceBuilder:
class AmoebaBondForceBuilder(BaseAmoebaForceBuilder):
"""Builder for Bond force for AMOEBA force field"""
"""
Builder for AmoebaBondForce for AMOEBA force field
Attributes
----------
cubic : float
Cubic term coefficient
quartic : float
Quartic term coefficient
bondParams : list
List of bond parameters as tuples of (bondType, (r0, k0)),
where bondType is a tuple of atom classes.
Parameters
----------
cubic : float
Cubic term coefficient
quartic : float
Quartic term coefficient
"""
def __init__(self, cubic, quartic):
def __init__(self, cubic: float, quartic: float) -> None:
super().__init__()
self.cubic = cubic
self.quartic = quartic
self.bondParams = []
def registerParams(self, bondType, r0, k0):
"""Register bond parameters"""
def registerParams(self, bondType: tuple, r0: float, k0: float) -> None:
"""
Register bond parameters.
Parameters
----------
bondType : tuple
A tuple of atom classes defining the bond type.
r0 : float
The equilibrium bond length.
k0 : float
The force constant.
"""
self.bondParams.append((bondType, (r0, k0)))
def getForce(self, sys):
def getForce(self, sys: mm.System) -> mm.CustomBondForce:
"""
Get or create the AmoebaBondForce in the system.
Parameters
----------
sys : mm.System
The system to get the force from or add it to.
Returns
-------
mm.CustomBondForce
The AmoebaBondForce instance.
"""
energy = f"k*(d^2 + {self.cubic}*d^3 + {self.quartic}*d^4); d=r-r0"
def createForce():
force = mm.CustomBondForce(energy)
force.addPerBondParameter("r0")
force.addPerBondParameter("k")
force.setName("AmoebaBond")
force.setName("AmoebaBondForce")
return force
return self._createOrGetForce(sys, mm.CustomBondForce, createForce, energyFunction=energy)
def addBonds(self, force, atomClasses, bonds):
"""Add bonds to the force"""
for atom1, atom2 in bonds:
atomTypes = (atomClasses[atom1], atomClasses[atom2])
params = self._findMatchingParams(self.bondParams, atomTypes)
if params is not None and params[1] != 0:
force.addBond(atom1, atom2, params)
def addBonds(self, force: mm.CustomBondForce, atomClasses: list, bonds: list, bondsConstraints: Optional[list] = None, flexibleConstraints: bool = False) -> None:
"""
Add bonds to the AmoebaBondForce.
Parameters
----------
force : mm.CustomBondForce
The AmoebaBondForce instance to add bonds to.
atomClasses : list
List of atom classes indexed by atom index.
bonds : list
List of bonds indices as tuples of (atom1, atom2).
isConstrained : list
List of flags indicating if a given bond is constrainted.
flexibleConstraints : bool
If True, constrained bonds will still be added to the system.
"""
for i, (atom1, atom2) in enumerate(bonds):
isConstrained = False if bondsConstraints is None else bondsConstraints[i]
if not isConstrained or flexibleConstraints:
atomTypes = (atomClasses[atom1], atomClasses[atom2])
params = self._findMatchingParams(self.bondParams, atomTypes)
if params is not None:
if params[1] != 0.0:
force.addBond(atom1, atom2, params)
class AmoebaAngleForceBuilder(BaseAmoebaForceBuilder):
"""
......@@ -289,12 +373,49 @@ class AmoebaAngleForceBuilder(BaseAmoebaForceBuilder):
force = mm.CustomAngleForce(energy)
force.addPerAngleParameter("theta0")
force.addPerAngleParameter("k")
force.setName("AmoebaAngle")
force.setName("AmoebaAngleForce")
return force
return self._createOrGetForce(sys, mm.CustomAngleForce, createForce, energyFunction=energy)
def addAngles(self, force: mm.CustomAngleForce, angles: List[Tuple[int, int, int]]) -> None:
def getIdealAngle(self, angle: Tuple[int, int, int], thetaList: List[float], atoms: List[Any], bondedToAtom: List[Any]) -> Optional[float]:
"""
Get the ideal angle for a given angle.
Notes
-----
Some angle definitions have multiple equilibrium values. For example:
angle 7 7 8 48.20 109.50 110.20 111.00
The appropriate equilibrium angle is determined by the number of hydrogens attached to the central atom that do not participate in the angle.
Parameters
----------
angle : Tuple[int, int, int]
A tuple of atom indices defining the angle.
thetaList : List[float]
List of ideal angles for different k-indices.
atoms : List[Any]
List of atoms indexed by atom index.
bondedToAtom : List[Any]
List of bonded atoms indexed by atom index.
Returns
-------
Optional[float]
The ideal angle if found, None otherwise.
"""
nTheta = len(thetaList)
if nTheta > 1:
partners = [at for at in bondedToAtom[angle[1]] if at not in angle]
nHyd = sum(1 for i in partners if self._getAtomicNumber(atoms[i]) == 1)
if nHyd < nTheta:
theta = thetaList[nHyd]
else:
raise ValueError(f"Angle parameters out of range for angle with indices {angle}")
else:
theta = thetaList[0]
return theta
def addAngles(self, force: mm.CustomAngleForce, angles: List[Tuple[int, int, int]], anglesConstraints: Optional[list] = None, flexibleConstraints: bool = False) -> None:
"""
Add angles to the force.
......@@ -304,20 +425,25 @@ class AmoebaAngleForceBuilder(BaseAmoebaForceBuilder):
The AmoebaAngleForce instance to add angles to.
angles : List[Tuple[int, int, int]]
List of angle indices as tuples of (atom1, atom2, atom3).
anglesConstraints : Optional[list]
List of flags indicating if a given angle is constrainted.
flexibleConstraints : bool
If True, constrained angles will still be added to the system.
"""
for atom1, atom2, atom3 in angles:
for angleType, params in self.angleParams:
if params[1] != 0:
if self._matchParams((atom1, atom2, atom3), angleType):
for i, (atom1, atom2, atom3) in enumerate(angles):
isConstrained = False if anglesConstraints is None else anglesConstraints[i]
if not isConstrained or flexibleConstraints:
angle = (atom1, atom2, atom3)
params = self._findMatchingParams(self.angleParams, angle)
if params is not None:
if params[1] != 0.0:
force.addAngle(atom1, atom2, atom3, params)
break
class AmoebaInPlaneAngleForceBuilder(BaseAmoebaForceBuilder):
"""
Builder for InPlaneAngleForce force for AMOEBA force field.
Attributes
----------
cubic : float
......@@ -401,12 +527,12 @@ class AmoebaInPlaneAngleForceBuilder(BaseAmoebaForceBuilder):
force = mm.CustomCompoundBondForce(4, energy)
force.addPerBondParameter("theta0")
force.addPerBondParameter("k")
force.setName("AmoebaInPlaneAngle")
force.setName("AmoebaInPlaneAngleForce")
return force
return self._createOrGetForce(sys, mm.CustomCompoundBondForce, createForce, energyFunction=energy)
def addInPlaneAngles(self, force: mm.CustomCompoundBondForce, atomClasses: List[Any], inPlaneAngles: List[Tuple[int, int, int, int]]) -> None:
def addInPlaneAngles(self, force: mm.CustomCompoundBondForce, atomClasses: List[Any], inPlaneAngles: List[Tuple[int, int, int, int]], anglesConstraints: Optional[list] = None, flexibleConstraints: bool = False) -> None:
"""
Add in-plane angles to the force.
......@@ -418,14 +544,19 @@ class AmoebaInPlaneAngleForceBuilder(BaseAmoebaForceBuilder):
List of atom classes indexed by atom index.
inPlaneAngles : List[Tuple[int, int, int, int]]
List of in-plane angle indices as tuples of (atom1, atom2, atom3, atom4).
anglesConstraints : Optional[list]
List of flags indicating if a given angle is constrainted.
flexibleConstraints : bool
If True, constrained angles will still be added to the system.
"""
for atom1, atom2, atom3, atom4 in inPlaneAngles:
for inPlaneAngleType, params in self.inPlaneAngleParams:
if params[1] != 0:
angleClasses = (atomClasses[atom1], atomClasses[atom2], atomClasses[atom3])
if self._matchParams(angleClasses, inPlaneAngleType[:3], reverseMatch=True):
for i, (atom1, atom2, atom3, atom4) in enumerate(inPlaneAngles):
isConstrained = False if anglesConstraints is None else anglesConstraints[i]
if not isConstrained or flexibleConstraints:
angleClasses = (atomClasses[atom1], atomClasses[atom2], atomClasses[atom3])
params = self._findMatchingParams(self.inPlaneAngleParams, angleClasses, reverseMatch=True)
if params is not None:
if params[1] != 0.0:
force.addBond([atom1, atom2, atom3, atom4], params)
break
class AmoebaOutOfPlaneBendForceBuilder(BaseAmoebaForceBuilder):
......@@ -486,7 +617,91 @@ class AmoebaOutOfPlaneBendForceBuilder(BaseAmoebaForceBuilder):
"""
at1, at2, at3, at4 = atomClasses
p1, p2, p3, p4 = paramClasses
return at2 == p2 and at4 == p1 and {at1, at3} == {at1 if p3 == "0" else p3, at3 if p4 == "0" else p4}
return at2 == p2 and at4 == p1 and {at1, at3} == {at1 if p3 == '' else p3, at3 if p4 == '' else p4}
@staticmethod
def classifyAngles(angles: List[Tuple[int, int, int]], atomClasses: List[Any],
bondedToAtom: List[List[int]], opbendTypes: List[Tuple[Any, ...]]) -> Tuple[
List[Tuple[int, ...]],
List[Tuple[int, int, int, int]],
List[Tuple[int, int, int]]]:
"""
Classify angles into in-plane, out-of-plane, and generic categories.
Notes
-----
For atoms with covalency 3, if the central atom and all partners match
out-of-plane bend parameters, the angles are classified as in-plane
with corresponding out-of-plane bends.
Parameters
----------
angles : List[Tuple[int, int, int]]
List of angles as tuples of (atom1, central_atom, atom2).
atomClasses : List[Any]
List of atom classes indexed by atom index.
bondedToAtom : List[List[int]]
List of bonded atoms for each atom index.
opbendTypes : List[Tuple[Any, ...]]
List of out-of-plane bend parameter types as (type1, type2, type3, type4).
Returns
-------
Tuple[List[Tuple[int, ...]], List[Tuple[int, int, int, int]], List[Tuple[int, int, int]]]
A tuple containing:
- inPlaneAngles: List of in-plane angles (atom1, central, atom2) or (atom1, central, atom2, atom3)
- outOfPlaneAngles: List of out-of-plane bends (atom1, central, atom3, atom4)
- genericAngles: List of generic angles (atom1, central, atom2)
"""
inPlaneAngles = []
outOfPlaneAngles = []
genericAngles = []
skipAtoms = {}
for angle in angles:
middleAtom = angle[1]
middleClass = atomClasses[middleAtom]
middleCovalency = len(bondedToAtom[middleAtom])
if middleCovalency == 3 and middleAtom not in skipAtoms:
partners = []
for partner in bondedToAtom[middleAtom]:
partnerClass = atomClasses[partner]
for opBendType in opbendTypes:
if middleClass in opBendType[1:2] and partnerClass in opBendType[0:1]:
partners.append(partner)
break
if len(partners) == 3:
outOfPlaneAngles.append((partners[0], middleAtom, partners[1], partners[2]))
outOfPlaneAngles.append((partners[2], middleAtom, partners[0], partners[1]))
outOfPlaneAngles.append((partners[1], middleAtom, partners[2], partners[0]))
skipAtoms[middleAtom] = set(partners[:3])
angleList = list(angle[:3])
for atomIndex in partners:
if atomIndex not in angleList:
angleList.append(atomIndex)
inPlaneAngles.append(tuple(angleList))
else:
angleList = list(angle[:3])
for atomIndex in partners:
if atomIndex not in angleList:
angleList.append(atomIndex)
genericAngles.append(tuple(angleList))
elif middleCovalency == 3 and middleAtom in skipAtoms:
angleList = list(angle[:3])
for atomIndex in skipAtoms[middleAtom]:
if atomIndex not in angleList:
angleList.append(atomIndex)
inPlaneAngles.append(tuple(angleList))
else:
genericAngles.append(angle)
return inPlaneAngles, outOfPlaneAngles, genericAngles
def registerParams(self, outOfPlaneBendType: Tuple[Any, ...], params: Tuple[float, ...]) -> None:
"""
......@@ -534,7 +749,7 @@ class AmoebaOutOfPlaneBendForceBuilder(BaseAmoebaForceBuilder):
def createForce():
force = mm.CustomCompoundBondForce(4, energy)
force.addPerBondParameter("k")
force.setName("AmoebaOutOfPlaneBend")
force.setName("AmoebaOutOfPlaneBendForce")
return force
return self._createOrGetForce(sys, mm.CustomCompoundBondForce, createForce, energyFunction=energy)
......@@ -554,11 +769,11 @@ class AmoebaOutOfPlaneBendForceBuilder(BaseAmoebaForceBuilder):
"""
for atom1, atom2, atom3, atom4 in outOfPlaneBends:
for outOfPlaneBendType, params in self.outOfPlaneBendParams:
if params[0] != 0:
angleClasses = (atomClasses[atom1], atomClasses[atom2], atomClasses[atom3], atomClasses[atom4])
if self._matchParams(angleClasses, outOfPlaneBendType):
angleClasses = (atomClasses[atom1], atomClasses[atom2], atomClasses[atom3], atomClasses[atom4])
if self._matchParams(angleClasses, outOfPlaneBendType):
if params[0] != 0.0:
force.addBond([atom1, atom2, atom3, atom4], params)
break
break
class AmoebaStretchBendForceBuilder(BaseAmoebaForceBuilder):
......@@ -568,11 +783,112 @@ class AmoebaStretchBendForceBuilder(BaseAmoebaForceBuilder):
super().__init__()
self.stretchBendParams = []
def registerParams(self, stretchBendType, params):
"""Register stretch-bend parameters"""
def registerParams(self, stretchBendType: tuple, params: tuple) -> None:
"""
Register stretch-bend parameters.
Parameters
----------
stretchBendType : tuple
A tuple identifying the type of stretch-bend interaction.
params : tuple
A tuple containing the parameters for the stretch-bend interaction.
"""
self.stretchBendParams.append((stretchBendType, params))
def getForce(self, sys):
def registerAllStretchBendParams(self,
atomClasses: list,
angles: list,
stretchBendParams: dict,
bondParams: dict,
idealAngles: dict
):
"""
Register all stretch-bend parameters for the given angles.
Notes
-----
For each angle defined by three atoms (1-2-3), this method checks if there are
corresponding stretch-bend parameters for the atom classes of these atoms. If such
parameters exist, it retrieves the ideal angle and bond parameters for the bonds
(1-2) and (2-3). It then registers the parameters for the stretch-bend interaction.
The method returns a list of angles for which stretch-bend parameters were successfully
registered.
Parameters
----------
atomClasses : list
List of atom classes indexed by atom index.
angles : list
List of angles defined by tuples of (atom1, atom2, atom3).
stretchBendParams : dict
Dictionary mapping atom class triplets to stretch-bend parameters.
bondParams : dict
Dictionary mapping atom class pairs to bond parameters.
idealAngles : dict
Dictionary mapping angle tuples to ideal angles.
Returns
-------
list
List of angles for which stretch-bend parameters were registered.
"""
processedAngles = []
for angle in angles:
class1 = atomClasses[angle[0]]
class2 = atomClasses[angle[1]]
class3 = atomClasses[angle[2]]
params = stretchBendParams.get((class1, class2, class3))
if params is not None:
swap = False
else:
params = stretchBendParams.get((class3, class2, class1))
swap = params is not None
if params:
idealAngle = idealAngles.get(tuple(angle), None)
if idealAngle is None:
continue
if not swap:
# 1-2-3
bondAB = bondParams.get((class1, class2)) or bondParams.get((class2, class1))
bondCB = bondParams.get((class2, class3)) or bondParams.get((class3, class2))
if bondAB and bondCB:
bondAB, bondCB = bondAB['r0'], bondCB['r0']
k1, k2 = params["k1"], params["k2"]
else:
# 3-2-1
bondAB = bondParams.get((class3, class2)) or bondParams.get((class2, class3))
bondCB = bondParams.get((class2, class1)) or bondParams.get((class1, class2))
if bondAB and bondCB:
bondAB, bondCB = bondCB['r0'], bondAB['r0']
k1, k2 = params["k2"], params["k1"]
if bondAB and bondCB:
# Because for the same class triplet, depending on the number of hydrogens on the central atom,
# the ideal angle can be different, we need to register parameters for each angle.
paramKey = (angle[0], angle[1], angle[2])
self.registerParams(paramKey, (bondAB, bondCB, idealAngle, k1, k2))
processedAngles.append((angle[0], angle[1], angle[2]))
return processedAngles
def getForce(self, sys: mm.System) -> mm.CustomCompoundBondForce:
"""
Get the force for the system.
Parameters
----------
sys : mm.System
The system to get the force for.
Returns
-------
mm.CustomCompoundBondForce
The force for the system.
"""
energy = (
"(k1*(distance(p1,p2)-r12) + k2*(distance(p2,p3)-r23))*(%.15g*(angle(p1,p2,p3)-theta0))"
% (180 / math.pi)
......@@ -585,19 +901,28 @@ class AmoebaStretchBendForceBuilder(BaseAmoebaForceBuilder):
force.addPerBondParameter("theta0")
force.addPerBondParameter("k1")
force.addPerBondParameter("k2")
force.setName("AmoebaStretchBend")
force.setName("AmoebaStretchBendForce")
return force
return self._createOrGetForce(sys, mm.CustomCompoundBondForce, createForce, energyFunction=energy)
def addStretchBends(self, force, angles):
"""Add stretch-bend terms to the force"""
def addStretchBends(self, force: mm.CustomCompoundBondForce, angles: list) -> None:
"""
Add stretch-bend terms to the force
Parameters
----------
force : mm.CustomCompoundBondForce
The force to add the stretch-bend terms to.
angles : list of tuples
The angles to add the stretch-bend terms for.
"""
for atom1, atom2, atom3 in angles:
for stretchBendType, params in self.stretchBendParams:
angle = (atom1, atom2, atom3)
if self._matchParams(angle, stretchBendType, reverseMatch=True):
force.addBond(angle, params)
break
angle = (atom1, atom2, atom3)
params = self._findMatchingParams(self.stretchBendParams, (atom1, atom2, atom3), reverseMatch=True)
if params[3] != 0.0 or params[4] != 0.0:
force.addBond(angle, params)
class AmoebaStretchTorsionForceBuilder(BaseAmoebaForceBuilder):
"""
......@@ -661,7 +986,7 @@ class AmoebaStretchTorsionForceBuilder(BaseAmoebaForceBuilder):
force.addPerBondParameter(f'length{i+1}')
for i in range(3):
force.addPerBondParameter(f'phase{i+1}')
force.setName('AmoebaStretchTorsion')
force.setName('AmoebaStretchTorsionForce')
return force
return self._createOrGetForce(sys, mm.CustomCompoundBondForce, createForce, energyFunction=energy)
......@@ -681,10 +1006,8 @@ class AmoebaStretchTorsionForceBuilder(BaseAmoebaForceBuilder):
torsions : List[Tuple[int, int, int, int]]
List of torsion indices as tuples of (atom1, atom2, atom3, atom4).
"""
# Record parameters for bonds and torsions so we can look them up quickly.
bondForce = [f for f in sys.getForces() if type(f) == mm.CustomBondForce and f.getName() == 'AmoebaBond'][0]
bondForce = [f for f in sys.getForces() if type(f) == mm.CustomBondForce and f.getName() == 'AmoebaBondForce'][0]
torsionForce = [f for f in sys.getForces() if type(f) == mm.PeriodicTorsionForce][0]
bondLength = {}
torsionPhase = defaultdict(lambda: [0.0, math.pi, 0.0])
......@@ -700,7 +1023,6 @@ class AmoebaStretchTorsionForceBuilder(BaseAmoebaForceBuilder):
torsionPhase[(p4, p3, p2, p1)][periodicity-1] = phase
# Add stretch-torsions.
for torsion in torsions:
atom1, atom2, atom3, atom4 = torsion
atomTypes = (atomClasses[atom1], atomClasses[atom2], atomClasses[atom3], atomClasses[atom4])
......@@ -713,7 +1035,8 @@ class AmoebaStretchTorsionForceBuilder(BaseAmoebaForceBuilder):
params.append(bondLength[(atom2, atom3)])
params.append(bondLength[(atom3, atom4)])
params += torsionPhase[torsion]
force.addBond((atom1, atom2, atom3, atom4), params)
if any(p != 0.0 for p in params[:9]):
force.addBond((atom1, atom2, atom3, atom4), params)
break
......@@ -776,7 +1099,7 @@ class AmoebaAngleTorsionForceBuilder(BaseAmoebaForceBuilder):
force.addPerBondParameter(f'angle{i+1}')
for i in range(3):
force.addPerBondParameter(f'phase{i+1}')
force.setName('AmoebaAngleTorsion')
force.setName('AmoebaAngleTorsionForce')
return force
return self._createOrGetForce(sys, mm.CustomCompoundBondForce, createForce, energyFunction=energy)
......@@ -796,11 +1119,9 @@ class AmoebaAngleTorsionForceBuilder(BaseAmoebaForceBuilder):
torsions : List[Tuple[int, int, int, int]]
List of torsion indices as tuples of (atom1, atom2, atom3, atom4).
"""
# Record parameters for angles and torsions so we can look them up quickly.
angleForce = [f for f in sys.getForces() if type(f) == mm.CustomAngleForce and f.getName() == 'AmoebaAngle'][0]
inPlaneAngleForce = [f for f in sys.getForces() if type(f) == mm.CustomCompoundBondForce and f.getName() == 'AmoebaInPlaneAngle'][0]
angleForce = [f for f in sys.getForces() if type(f) == mm.CustomAngleForce and f.getName() == 'AmoebaAngleForce'][0]
inPlaneAngleForce = [f for f in sys.getForces() if type(f) == mm.CustomCompoundBondForce and f.getName() == 'AmoebaInPlaneAngleForce'][0]
torsionForce = [f for f in sys.getForces() if type(f) == mm.PeriodicTorsionForce][0]
equilAngle = {}
torsionPhase = defaultdict(lambda: [0.0, math.pi, 0.0])
......@@ -821,7 +1142,6 @@ class AmoebaAngleTorsionForceBuilder(BaseAmoebaForceBuilder):
torsionPhase[(p4, p3, p2, p1)][periodicity-1] = phase
# Add angle-torsions.
for torsion in torsions:
atom1, atom2, atom3, atom4 = torsion
atomTypes = (atomClasses[atom1], atomClasses[atom2], atomClasses[atom3], atomClasses[atom4])
......@@ -833,7 +1153,8 @@ class AmoebaAngleTorsionForceBuilder(BaseAmoebaForceBuilder):
params.append(equilAngle[(atom1, atom2, atom3)])
params.append(equilAngle[(atom2, atom3, atom4)])
params += torsionPhase[torsion]
force.addBond((atom1, atom2, atom3, atom4), params)
if any(p != 0.0 for p in params[:6]):
force.addBond((atom1, atom2, atom3, atom4), params)
break
......@@ -895,18 +1216,17 @@ class AmoebaTorsionForceBuilder(BaseAmoebaForceBuilder):
List of torsion indices as tuples of (atom1, atom2, atom3, atom4).
"""
for atom1, atom2, atom3, atom4 in torsions:
for torsionType, params in self.torsionParams:
atomTypes = (atomClasses[atom1], atomClasses[atom2],
atomClasses[atom3], atomClasses[atom4])
if self._matchParams(atomTypes, torsionType):
t1, t2, t3 = params
if t1[0] != 0:
force.addTorsion(atom1, atom2, atom3, atom4, 1, t1[1], t1[0])
if t2[0] != 0:
force.addTorsion(atom1, atom2, atom3, atom4, 2, t2[1], t2[0])
if t3[0] != 0:
force.addTorsion(atom1, atom2, atom3, atom4, 3, t3[1], t3[0])
break
torsionType = (atomClasses[atom1], atomClasses[atom2], atomClasses[atom3], atomClasses[atom4])
params = self._findMatchingParams(self.torsionParams, torsionType)
if params is not None:
t1, t2, t3 = params
if t1[0] != 0:
force.addTorsion(atom1, atom2, atom3, atom4, 1, t1[1], t1[0])
if t2[0] != 0:
force.addTorsion(atom1, atom2, atom3, atom4, 2, t2[1], t2[0])
if t3[0] != 0:
force.addTorsion(atom1, atom2, atom3, atom4, 3, t3[1], t3[0])
continue
class AmoebaPiTorsionForceBuilder(BaseAmoebaForceBuilder):
......@@ -963,11 +1283,47 @@ class AmoebaPiTorsionForceBuilder(BaseAmoebaForceBuilder):
def createForce():
force = mm.CustomCompoundBondForce(6, energy)
force.addPerBondParameter("k")
force.setName("AmoebaPiTorsion")
force.setName("AmoebaPiTorsionForce")
return force
return self._createOrGetForce(sys, mm.CustomCompoundBondForce, createForce, energyFunction=energy)
def getAllPiTorsions(self, atomClasses: List[str], bondedToAtom: List[Tuple[int, ...]], bonds: List[Tuple[int, int]]) -> List[Tuple[int, int, int, int, int, int]]:
"""
Get all pi-torsion indices from the bonds.
Parameters
----------
atomClasses : List[str]
List of atom classes indexed by atom index.
bondedToAtom : List[Any]
List of atoms bonded to each atom.
bonds : List[Tuple[int, int]]
List of bonds defined by tuples of (atom1, atom2).
Returns
-------
List[Tuple[int, int, int, int, int, int]]
List of pi-torsion indices as tuples of six atoms.
"""
processedPiTorsions = []
for (at1, at2) in bonds:
valence1 = len(bondedToAtom[at1])
valence2 = len(bondedToAtom[at2])
if valence1 == 3 and valence2 == 3:
class1 = atomClasses[at1]
class2 = atomClasses[at2]
params = self._findMatchingParams(self.piTorsionParams, (class1, class2))
if params:
piTorsionAtom3 = at1
piTorsionAtom4 = at2
# piTorsionAtom1, piTorsionAtom2 are the atoms bonded to atom1, excluding atom2
# piTorsionAtom5, piTorsionAtom6 are the atoms bonded to atom2, excluding atom1
piTorsionAtom1, piTorsionAtom2 = [at for at in bondedToAtom[at1] if at != piTorsionAtom4]
piTorsionAtom5, piTorsionAtom6 = [at for at in bondedToAtom[at2] if at != piTorsionAtom3]
processedPiTorsions.append((piTorsionAtom1, piTorsionAtom2, piTorsionAtom3, piTorsionAtom4, piTorsionAtom5, piTorsionAtom6))
return processedPiTorsions
def addPiTorsions(self, force: mm.CustomCompoundBondForce, atomClasses: List[Any], piTorsions: List[Tuple[int, int, int, int, int, int]]) -> None:
"""
Add pi-torsions to the force.
......@@ -982,12 +1338,11 @@ class AmoebaPiTorsionForceBuilder(BaseAmoebaForceBuilder):
List of pi-torsion indices as tuples of six atoms.
"""
for atom1, atom2, atom3, atom4, atom5, atom6 in piTorsions:
for piTorsionType, params in self.piTorsionParams:
if params[0] != 0:
types = (atomClasses[atom3], atomClasses[atom4])
if self._matchParams(types, piTorsionType, reverseMatch=True):
force.addBond([atom1, atom2, atom3, atom4, atom5, atom6], params)
break
piTorsionType = (atomClasses[atom3], atomClasses[atom4])
params = self._findMatchingParams(self.piTorsionParams, piTorsionType)
if params is not None:
if params[0] != 0.0:
force.addBond([atom1, atom2, atom3, atom4, atom5, atom6], params)
class AmoebaUreyBradleyForceBuilder(BaseAmoebaForceBuilder):
......@@ -1003,9 +1358,9 @@ class AmoebaUreyBradleyForceBuilder(BaseAmoebaForceBuilder):
def __init__(self) -> None:
"""Initialize the AmoebaUreyBradleyForceBuilder."""
super().__init__()
self.ureyBradleyParams = []
self.ureyBradleyParams = {True: [], False: []} # True = class-based, False = type-based
def registerParams(self, ureyBradleyType: Tuple[Any, ...], params: Tuple[float, float]) -> None:
def registerParams(self, ureyBradleyType: Tuple[Any, ...], params: Tuple[float, float], isClass: bool) -> None:
"""
Register Urey-Bradley parameters.
......@@ -1015,8 +1370,10 @@ class AmoebaUreyBradleyForceBuilder(BaseAmoebaForceBuilder):
A tuple of atom classes defining the Urey-Bradley type.
params : Tuple[float, float]
The Urey-Bradley parameters (force constant, equilibrium distance).
isClass : bool
Indicates if these parameters are defined for atom classes (True) or atom types (False).
"""
self.ureyBradleyParams.append((ureyBradleyType, params))
self.ureyBradleyParams[isClass].append((ureyBradleyType, params))
def getForce(self, sys: mm.System) -> mm.HarmonicBondForce:
"""
......@@ -1034,7 +1391,7 @@ class AmoebaUreyBradleyForceBuilder(BaseAmoebaForceBuilder):
"""
return self._createOrGetForce(sys, mm.HarmonicBondForce, mm.HarmonicBondForce)
def addUreyBradleys(self, force: mm.HarmonicBondForce, atomClasses: List[Any], angles: List[Tuple[int, int, int]], anglesConstraints: Optional[List[bool]] = None, flexibleConstraints: bool = False) -> None:
def addUreyBradleys(self, force: mm.HarmonicBondForce, atomClassesOrTypes: List[Any], angles: List[Tuple[int, int, int]], isClass: bool, anglesConstraints: Optional[List[bool]] = None, flexibleConstraints: bool = False) -> None:
"""
Add Urey-Bradley terms to the force.
......@@ -1042,42 +1399,88 @@ class AmoebaUreyBradleyForceBuilder(BaseAmoebaForceBuilder):
----------
force : mm.HarmonicBondForce
The HarmonicBondForce instance to add Urey-Bradley terms to.
atomClasses : List[Any]
List of atom classes indexed by atom index.
atomClassesOrTypes : List[Any]
List of atom classes or types indexed by atom index.
angles : List[Tuple[int, int, int]]
List of angle indices as tuples of (atom1, atom2, atom3).
isClass : bool
Indicates if the parameters are defined for atom classes (True) or atom types (False).
anglesConstraints : Optional[List[bool]], default=None
List of flags indicating if a given angle is constrained.
flexibleConstraints : bool, default=False
If True, constrained angles will still be added to the system.
"""
paramSet = self.ureyBradleyParams[isClass]
for i, (atom1, atom2, atom3) in enumerate(angles):
isConstrained = False if anglesConstraints is None else anglesConstraints[i]
if not isConstrained or flexibleConstraints:
ubType = (atomClasses[atom1], atomClasses[atom2], atomClasses[atom3])
params = self._findMatchingParams(self.ureyBradleyParams, ubType, reverseMatch=True)
ubIdentifier = (atomClassesOrTypes[atom1], atomClassesOrTypes[atom2], atomClassesOrTypes[atom3])
params = self._findMatchingParams(paramSet, ubIdentifier, reverseMatch=True)
if params is not None:
k, d = params
force.addBond(atom1, atom3, d, 2 * k)
if k != 0.0:
force.addBond(atom1, atom3, d, 2 * k)
class AmoebaTorsionTorsionForceBuilder(BaseAmoebaForceBuilder):
"""Builder for TorsionTorsion force for AMOEBA force field"""
"""
Builder for TorsionTorsion force for AMOEBA force field.
Attributes
----------
torsionTorsionParams : List[Tuple[Tuple[str, ...], int]]
List of registered torsion-torsion parameters as (atom_class_pattern, grid_index) tuples.
gridData : Dict[int, Any]
Dictionary mapping grid indices to grid data for torsion-torsion interactions.
"""
def __init__(self):
def __init__(self) -> None:
"""Initialize the AmoebaTorsionTorsionForceBuilder."""
super().__init__()
self.torsionTorsionParams = []
self.gridData = {}
self.torsionTorsionParams: List[Tuple[Tuple[str, ...], int]] = []
self.gridData: Dict[int, Any] = {}
def registerParams(self, torsionTorsionType, params):
"""Register torsion-torsion parameters"""
def registerParams(self, torsionTorsionType: Tuple[str, ...], params: int) -> None:
"""
Register torsion-torsion parameters during XML parsing.
Parameters
----------
torsionTorsionType : Tuple[str, ...]
Tuple of atom class strings defining the torsion-torsion type pattern.
params : int
Grid index for the torsion-torsion interaction.
"""
self.torsionTorsionParams.append((torsionTorsionType, params))
def registerGridData(self, gridIndex, grid):
"""Register grid data for torsion-torsion interactions"""
def registerGridData(self, gridIndex: int, grid: Any) -> None:
"""
Register grid data for torsion-torsion interactions.
Parameters
----------
gridIndex : int
Index of the grid in the force.
grid : Any
Grid data object for the torsion-torsion interaction.
"""
self.gridData[gridIndex] = grid
def getForce(self, sys):
def getForce(self, sys: mm.System) -> mm.AmoebaTorsionTorsionForce:
"""
Get or create the AmoebaTorsionTorsionForce in the system.
Parameters
----------
sys : mm.System
The OpenMM System to get the force from or add it to.
Returns
-------
mm.AmoebaTorsionTorsionForce
The AmoebaTorsionTorsionForce instance with grid data registered.
"""
def createForce():
force = mm.AmoebaTorsionTorsionForce()
return force
......@@ -1089,245 +1492,106 @@ class AmoebaTorsionTorsionForceBuilder(BaseAmoebaForceBuilder):
return force
@staticmethod
def setTorsionTorsionGrid(force, gridIndex, grid):
"""Set the torsion-torsion grid for the force."""
force.setTorsionTorsionGrid(gridIndex, grid)
@staticmethod
def addTorsionTorsion(force, atom1, atom2, atom3, atom4, atom5, chiralIndex, gridIndex):
"""Add a torsion-torsion interaction to the force."""
force.addTorsionTorsion(atom1, atom2, atom3, atom4, atom5, chiralIndex, gridIndex)
@staticmethod
def addTorsionTorsionInteractions(force, data, types1, types2, types3, types4, types5, gridIndex, sys):
"""Add torsion-torsion interactions based on TINKER subroutine bitors()"""
for angle in data.angles:
# search for bitorsions; based on TINKER subroutine bitors()
ib = angle[0]
ic = angle[1]
id = angle[2]
for bondIndex in data.atomBonds[ib]:
bondedAtom1 = data.bonds[bondIndex].atom1
bondedAtom2 = data.bonds[bondIndex].atom2
if (bondedAtom1 != ib):
ia = bondedAtom1
else:
ia = bondedAtom2
if (ia != ic and ia != id):
for bondIndex2 in data.atomBonds[id]:
bondedAtom1 = data.bonds[bondIndex2].atom1
bondedAtom2 = data.bonds[bondIndex2].atom2
if (bondedAtom1 != id):
ie = bondedAtom1
else:
ie = bondedAtom2
if (ie != ic and ie != ib and ie != ia):
# found candidate set of atoms
# check if types match in order or reverse order
type1 = data.atomType[data.atoms[ia]]
type2 = data.atomType[data.atoms[ib]]
type3 = data.atomType[data.atoms[ic]]
type4 = data.atomType[data.atoms[id]]
type5 = data.atomType[data.atoms[ie]]
for i in range(len(types1)):
types1_i = types1[i]
types2_i = types2[i]
types3_i = types3[i]
types4_i = types4[i]
types5_i = types5[i]
# match in order
if (type1 in types1_i and type2 in types2_i and type3 in types3_i and
type4 in types4_i and type5 in types5_i):
chiralAtomIndex = AmoebaTorsionTorsionForceBuilder._getChiralAtomIndex(data, ib, ic, id)
force.addTorsionTorsion(ia, ib, ic, id, ie, chiralAtomIndex, gridIndex[i])
# match in reverse order
elif (type5 in types1_i and type4 in types2_i and type3 in types3_i and
type2 in types4_i and type1 in types5_i):
chiralAtomIndex = AmoebaTorsionTorsionForceBuilder._getChiralAtomIndex(data, ib, ic, id)
force.addTorsionTorsion(ie, id, ic, ib, ia, chiralAtomIndex, gridIndex[i])
@staticmethod
def createTorsionTorsionInteractions(force, angles, atoms, tortorParams):
def addTorsionTorsionInteractions(self, force: mm.AmoebaTorsionTorsionForce, torsions: List[Tuple[int, ...]], atomClasses: List[str], masses: List[float]) -> None:
"""
Create torsion-torsion interactions based on the angles and tortor parameters.
Create torsion-torsion interactions for both ForceField and TinkerFiles usage.
Parameters
----------
force : mm.AmoebaTorsionTorsionForce
The OpenMM AmoebaTorsionTorsionForce object
angles : list
List of angles in the system
atoms : list
List of TinkerAtom objects
tortorParams : list
List of torsion-torsion parameter sets
"""
The OpenMM AmoebaTorsionTorsionForce to add interactions to.
torsions : List[Tuple[int, ...]]
List of torsion tuples containing atom indices.
atomClasses : List[str]
List of atom class strings indexed by atom index.
masses : List[float]
List of atom masses indexed by atom index.
"""
# Determine bonded atoms from the torsions
atomBonds = [[] for _ in range(len(atomClasses))]
for torsion in torsions:
for i in range(len(torsion) - 1):
atom1, atom2 = torsion[i], torsion[i+1]
if atom2 not in atomBonds[atom1]:
atomBonds[atom1].append(atom2)
if atom1 not in atomBonds[atom2]:
atomBonds[atom2].append(atom1)
bonds = [(atom1, atom2) for atom1, bonded in enumerate(atomBonds) for atom2 in bonded if atom1 < atom2]
# Get angles from torsions
angles_set = set()
for torsion in torsions:
for i in range(len(torsion) - 2):
angles_set.add((torsion[i], torsion[i+1], torsion[i+2]))
angles = list(angles_set)
# Create torsion-torsion interactions
addedInteractions = set()
for angle in angles:
# angle = (atom1, atom2, atom3) where atom2 is the central atom
ib = angle[0] # first atom of angle
ic = angle[1] # central atom of angle
id = angle[2] # last atom of angle
# Find atoms bonded to ib (excluding ic and id)
for ia in atoms[ib].bonds:
if ia != ic and ia != id:
# Find atoms bonded to id (excluding ic and ib and ia)
for ie in atoms[id].bonds:
ib, ic, id = angle
for ia in atomBonds[ib]:
if ia != ic:
for ie in atomBonds[id]:
if ie != ic and ie != ib and ie != ia:
# We have a potential torsion-torsion pattern: ia-ib-ic-id-ie
# Check if atom types match any tortor parameters
for gridIndex, (tortorInfo, gridData) in enumerate(tortorParams):
class1_param = tortorInfo[0]
class2_param = tortorInfo[1]
class3_param = tortorInfo[2]
class4_param = tortorInfo[3]
class5_param = tortorInfo[4]
type1 = atoms[ia].atomClass
type2 = atoms[ib].atomClass
type3 = atoms[ic].atomClass
type4 = atoms[id].atomClass
type5 = atoms[ie].atomClass
# Check forward direction
if (type1 == class1_param and type2 == class2_param and
type3 == class3_param and type4 == class4_param and type5 == class5_param):
# Find chiral atom index
chiralIndex = AmoebaTorsionTorsionForceBuilder._getChiralAtomIndex(atoms, ib, ic, id)
AmoebaTorsionTorsionForceBuilder.addTorsionTorsion(force, ia, ib, ic, id, ie, chiralIndex, gridIndex)
# Check reverse direction
elif (type5 == class1_param and type4 == class2_param and
type3 == class3_param and type2 == class4_param and type1 == class5_param):
# Find chiral atom index
chiralIndex = AmoebaTorsionTorsionForceBuilder._getChiralAtomIndex(atoms, ib, ic, id)
AmoebaTorsionTorsionForceBuilder.addTorsionTorsion(force, ie, id, ic, ib, ia, chiralIndex, gridIndex)
interaction = tuple(sorted([(ia, ib, ic, id, ie), (ie, id, ic, ib, ia)]))
if interaction in addedInteractions:
continue
atomTypes = (atomClasses[ia], atomClasses[ib], atomClasses[ic], atomClasses[id], atomClasses[ie])
for torsionTorsionType, gridIdx in self.torsionTorsionParams:
if atomTypes == torsionTorsionType:
chiralIndex = AmoebaTorsionTorsionForceBuilder._getChiralIndex(ib, ic, id, atomClasses, bonds, masses)
force.addTorsionTorsion(ia, ib, ic, id, ie, chiralIndex, gridIdx)
addedInteractions.add(interaction)
break
elif atomTypes == tuple(reversed(torsionTorsionType)):
chiralIndex = AmoebaTorsionTorsionForceBuilder._getChiralIndex(ib, ic, id, atomClasses, bonds, masses)
force.addTorsionTorsion(ie, id, ic, ib, ia, chiralIndex, gridIdx)
addedInteractions.add(interaction)
break
@staticmethod
def _getChiralAtomIndex(data, atomB, atomC, atomD):
def _getChiralIndex(atomB: int, atomC: int, atomD: int, atomClasses: List[str], bonds: List[Tuple[int, int]], masses: List[float]) -> int:
"""
Get the chiral atom index based on the TINKER algorithm.
Get chiral atom index.
Parameters
----------
data : ForceFieldData or list
ForceFieldData object (for ForceField usage) or list of TinkerAtom objects (for TinkerFiles usage)
atomB : int
Index of atom B
Index of first atom in the central angle.
atomC : int
Index of atom C (central atom)
Index of central atom that should have 4 bonds for chiral determination.
atomD : int
Index of atom D
Index of third atom in the central angle.
atomClasses : List[str]
List of atom class strings indexed by atom index.
bonds : List[Tuple[int, int]]
List of bonds defined by tuples of (atom1, atom2).
masses : List[float]
List of particle masses indexed by atom index.
Returns
-------
int
Index of the chiral atom, or -1 if not found
"""
chiralAtomIndex = -1
# Check if we're dealing with ForceField data (has atomBonds attribute) or TinkerAtom list
if hasattr(data, 'atomBonds'):
# ForceField usage - data is ForceFieldData object
atoms = data.atoms
atomBonds = data.atomBonds
bonds = data.bonds
atomType = data.atomType
# Get atoms bonded to atomC
bondedAtoms = []
for bondIndex in atomBonds[atomC]:
bond = bonds[bondIndex]
if bond.atom1 == atomC:
bondedAtoms.append(bond.atom2)
else:
bondedAtoms.append(bond.atom1)
# If atomC has four bonds, find the two bonds that do not include atomB and atomD
if len(bondedAtoms) == 4:
atomE = -1
atomF = -1
for bondedAtom in bondedAtoms:
if bondedAtom != atomB and bondedAtom != atomD:
if atomE == -1:
atomE = bondedAtom
else:
atomF = bondedAtom
# Raise error if atoms E or F not found
if atomE == -1 or atomF == -1:
raise ValueError(f"getChiralAtomIndex: error getting bonded partners of atomC={atomC}")
# Check for different type between atoms E & F
typeE = atomType[atoms[atomE]]
typeF = atomType[atoms[atomF]]
if typeE and typeF:
# Compare atom types as strings/identifiers
if str(typeE) > str(typeF):
chiralAtomIndex = atomE
elif str(typeF) > str(typeE):
chiralAtomIndex = atomF
# If types are same, check masses
if chiralAtomIndex == -1:
massE = atoms[atomE].element.mass.value_in_unit(amu) if atoms[atomE].element else 0.0
massF = atoms[atomF].element.mass.value_in_unit(amu) if atoms[atomF].element else 0.0
if massE > massF:
chiralAtomIndex = atomE
elif massF > massE:
chiralAtomIndex = atomF
else:
# TinkerFiles usage - data is list of TinkerAtom objects
atoms = data
# If atomC has four bonds, find the two bonds that do not include atomB and atomD
# Set chiralAtomIndex to one of these, if they are not the same atom (type/mass)
atomC_bonds = atoms[atomC].bonds
if len(atomC_bonds) == 4:
atomE = -1
atomF = -1
for bondedAtom in atomC_bonds:
if bondedAtom != atomB and bondedAtom != atomD:
if atomE == -1:
atomE = bondedAtom
else:
atomF = bondedAtom
# Raise error if atoms E or F not found
if atomE == -1 or atomF == -1:
raise ValueError(f"getChiralAtomIndex: error getting bonded partners of atomC={atomC}")
# Check for different type/mass between atoms E & F
typeE = atoms[atomE].atomClass
typeF = atoms[atomF].atomClass
if typeE and typeF:
if int(typeE) > int(typeF):
chiralAtomIndex = atomE
elif int(typeF) > int(typeE):
chiralAtomIndex = atomF
# If types are same, check masses
if chiralAtomIndex == -1:
massE = atoms[atomE].mass if atoms[atomE].mass else 0.0
massF = atoms[atomF].mass if atoms[atomF].mass else 0.0
if massE > massF:
chiralAtomIndex = atomE
elif massF > massE:
chiralAtomIndex = atomF
return chiralAtomIndex
Index of the chiral atom (atomE or atomF) based on TINKER algorithm,
or -1 if atomC doesn't have exactly 4 bonds or other conditions aren't met.
"""
if len(bonds[atomC]) == 4:
otherAtoms = [atom for atom in bonds(atomC) if atom != atomB and atom != atomD]
if len(otherAtoms) == 2:
atomE, atomF = otherAtoms
typeE, typeF = atomClasses[atomE], atomClasses[atomF]
if typeE > typeF:
return atomE
elif typeF > typeE:
return atomF
massE, massF = masses[atomE], masses[atomF]
if massE > massF:
return atomE
elif massF > massE:
return atomF
return -1
class AmoebaWcaDispersionForceBuilder(BaseAmoebaForceBuilder):
......@@ -1933,17 +2197,15 @@ class AmoebaMultipoleForceBuilder(BaseAmoebaForceBuilder):
savedMultipoleParams = multipoleParams
hit = True
# assign multipole parameters via 1-2 and 1-3 connected atoms
# Assign multipole parameters via 1-2 and 1-3 connected atoms
for multipoleParams in multipoleList:
if hit:
break
kz, kx, ky = multipoleParams.kIndices[1:4]
# assign multipole parameters
# Assign multipole parameters
# (1) get bonded partners
# (2) match parameter types
bondedAtom12Indices = self.bonded12ParticleSets[atomIndex]
bondedAtom13Indices = self.bonded13ParticleSets[atomIndex]
zaxis = -1
......@@ -1963,8 +2225,7 @@ class AmoebaMultipoleForceBuilder(BaseAmoebaForceBuilder):
zaxis = bondedAtomZIndex
xaxis = bondedAtomXIndex
# select xaxis w/ smallest index
# Select xaxis w/ smallest index
for bondedAtomXIndex2 in bondedAtom13Indices:
bondedAtomX1Type = atomTypes[bondedAtomXIndex2]
if bondedAtomX1Type == kx and bondedAtomXIndex2 != bondedAtomZIndex and bondedAtomZIndex in self.bonded12ParticleSets[bondedAtomXIndex2] and bondedAtomXIndex2 < xaxis:
......@@ -1984,8 +2245,7 @@ class AmoebaMultipoleForceBuilder(BaseAmoebaForceBuilder):
savedMultipoleParams = multipoleParams
hit = True
# assign multipole parameters via only a z-defining atom
# Assign multipole parameters via only a z-defining atom
for multipoleParams in multipoleList:
if hit:
break
......@@ -2002,8 +2262,7 @@ class AmoebaMultipoleForceBuilder(BaseAmoebaForceBuilder):
savedMultipoleParams = multipoleParams
hit = True
# assign multipole parameters via no connected atoms
# Assign multipole parameters via no connected atoms
for multipoleParams in multipoleList:
if hit:
break
......@@ -2036,21 +2295,18 @@ class AmoebaMultipoleForceBuilder(BaseAmoebaForceBuilder):
raise ValueError('No multipole type for atom %d (%s %s %d)' % (atomIndex, atom.name, atom.residue.name, atom.residue.index))
# Set up the polarization groups.
self.setPolarGroups(force, atomTypes)
def buildBondedParticleSets(self, numAtoms, bonds):
"""Identify sets of particles that are separated by various numbers of bonds."""
# 1-2
bonded12ParticleSets = [set() for _ in range(numAtoms)]
for atom1, atom2 in bonds:
bonded12ParticleSets[atom1].add(atom2)
bonded12ParticleSets[atom2].add(atom1)
# 1-3
bonded13ParticleSets = []
for i in range(numAtoms):
bonded13Set = set()
......@@ -2059,7 +2315,6 @@ class AmoebaMultipoleForceBuilder(BaseAmoebaForceBuilder):
bonded13Set = bonded13Set.union(bonded12ParticleSets[j])
# remove 1-2 and self from set
bonded13Set = bonded13Set - bonded12ParticleSet
selfSet = set()
selfSet.add(i)
......@@ -2068,7 +2323,6 @@ class AmoebaMultipoleForceBuilder(BaseAmoebaForceBuilder):
bonded13ParticleSets.append(bonded13Set)
# 1-4
bonded14ParticleSets = []
for i in range(numAtoms):
bonded14Set = set()
......@@ -2077,7 +2331,6 @@ class AmoebaMultipoleForceBuilder(BaseAmoebaForceBuilder):
bonded14Set = bonded14Set.union(bonded12ParticleSets[j])
# remove 1-3, 1-2 and self from set
bonded14Set = bonded14Set - bonded12ParticleSets[i]
bonded14Set = bonded14Set - bonded13ParticleSet
selfSet = set()
......@@ -2087,7 +2340,6 @@ class AmoebaMultipoleForceBuilder(BaseAmoebaForceBuilder):
bonded14ParticleSets.append(bonded14Set)
# 1-5
bonded15ParticleSets = []
for i in range(numAtoms):
bonded15Set = set()
......@@ -2096,7 +2348,6 @@ class AmoebaMultipoleForceBuilder(BaseAmoebaForceBuilder):
bonded15Set = bonded15Set.union(bonded12ParticleSets[j])
# remove 1-4, 1-3, 1-2 and self from set
bonded15Set = bonded15Set - bonded12ParticleSets[i]
bonded15Set = bonded15Set - bonded13ParticleSets[i]
bonded15Set = bonded15Set - bonded14ParticleSet
......@@ -2118,7 +2369,6 @@ class AmoebaMultipoleForceBuilder(BaseAmoebaForceBuilder):
for atomIndex in range(len(atomTypes)):
# assign multipole parameters via only 1-2 connected atoms
groupTypes = self.polarizationParams[atomTypes[atomIndex]].groupAtomTypes
bondedAtomIndices = self.bonded12ParticleSets[atomIndex]
polarizationGroupForAtom[atomIndex].add(atomIndex)
......@@ -2129,7 +2379,6 @@ class AmoebaMultipoleForceBuilder(BaseAmoebaForceBuilder):
polarizationGroupForAtom[bondedAtomIndex].add(atomIndex)
# pgrp11
for atomIndex in range(len(atomTypes)):
if len(groupSetsForAtom[atomIndex]) > 0:
continue
......@@ -2157,7 +2406,6 @@ class AmoebaMultipoleForceBuilder(BaseAmoebaForceBuilder):
force.setCovalentMap(atomIndex, mm.AmoebaMultipoleForce.PolarizationCovalent11, groupSetsForAtom[atomIndex][0])
# pgrp12
for atomIndex in range(len(atomTypes)):
if len(groupSetsForAtom[atomIndex]) > 1:
continue
......@@ -2175,7 +2423,6 @@ class AmoebaMultipoleForceBuilder(BaseAmoebaForceBuilder):
force.setCovalentMap(atomIndex, mm.AmoebaMultipoleForce.PolarizationCovalent12, groupSetsForAtom[atomIndex][1])
# pgrp13
for atomIndex in range(len(atomTypes)):
if len(groupSetsForAtom[atomIndex]) > 2:
continue
......@@ -2195,7 +2442,6 @@ class AmoebaMultipoleForceBuilder(BaseAmoebaForceBuilder):
force.setCovalentMap(atomIndex, mm.AmoebaMultipoleForce.PolarizationCovalent13, groupSetsForAtom[atomIndex][2])
# pgrp14
for atomIndex in range(len(atomTypes)):
if len(groupSetsForAtom[atomIndex]) > 3:
continue
......@@ -2376,9 +2622,7 @@ class AmoebaVdwForceBuilder(BaseAmoebaForceBuilder):
bonds : List[Tuple[int, int]]
List of bond indices as tuples of (atom1, atom2).
"""
# Define types
sigmaScale = 1
if self.radiustype == 'SIGMA':
sigmaScale = 1.122462048309372
......@@ -2390,19 +2634,17 @@ class AmoebaVdwForceBuilder(BaseAmoebaForceBuilder):
classToTypeMap[className] = force.addParticleType(sigma*sigmaScale, epsilon)
# Record what other atoms each atom is bonded to.
bondedParticleSets = [set() for _ in range(len(atoms))]
for atom1, atom2 in bonds:
bondedParticleSets[atom1].add(atom2)
bondedParticleSets[atom2].add(atom1)
# add particles to force
# Add particles to force
for i, atom in enumerate(atoms):
className = atomClasses[i]
if className not in self.classParams:
raise ValueError(f"AmoebaVdwForce: No VdW parameters found for atom class {className}. "
f"Available classes: {list(self.classParams.keys())}")
f"Available classes: {list(self.classParams.keys())}")
_, _, reduction = self.classParams[className]
# ivIndex = index of bonded partner for hydrogens; otherwise ivIndex = particle index
ivIndex = i
......@@ -2413,28 +2655,22 @@ class AmoebaVdwForceBuilder(BaseAmoebaForceBuilder):
force.addParticle(ivIndex, classToTypeMap[className], reduction)
# Add pairs
for c1, c2, sigma, epsilon in self.pairParams:
force.addTypePair(classToTypeMap[c1], classToTypeMap[c2], sigma, epsilon)
# set combining rules
# set particle exclusions: self, 1-2 and 1-3 bonds
# (1) collect in bondedParticleSets[i], 1-2 indices for all bonded partners of particle i
# (2) add 1-2,1-3 and self to exclusion set
for i in range(len(atoms)):
# 1-2 partners
exclusionSet = bondedParticleSets[i].copy()
# 1-3 partners
if self.vdw13Scale == 0.0:
for bondedParticle in bondedParticleSets[i]:
exclusionSet = exclusionSet.union(bondedParticleSets[bondedParticle])
# self
exclusionSet.add(i)
force.setParticleExclusions(i, tuple(exclusionSet))
wrappers/python/openmm/app/tinkerfiles.py
View file @ 489e2c46
......@@ -523,22 +523,23 @@ class TinkerFiles:
# List of bonds as (atom1_index, atom2_index) used by various force builders
bonds = [(at1.index, at2.index) for at1, at2 in self.topology.bonds()]
# List of atoms bonded to each atom used by various force builders
bondedToAtom = [[] for _ in range(self.topology.getNumAtoms())]
for (at1, at2) in bonds:
bondedToAtom[at1].append(at2)
bondedToAtom[at2].append(at1)
# Add AmoebaBondForce
if "bond" in self._forces:
bondParams = {(at1, at2): {"k": float(k), "r0": float(r0)} for at1, at2, k, r0 in self._forces["bond"]}
bondParams = {(at1, at2): {"k": float(k)*100.0*4.184, "r0": float(r0)*0.1} for at1, at2, k, r0 in self._forces["bond"]}
bondForceBuilder = AmoebaBondForceBuilder(float(self._scalars["bond-cubic"])*10, float(self._scalars["bond-quartic"])*100)
for (class1, class2), params in bondParams.items():
bondForceBuilder.registerParams((class1, class2), params["r0"]*0.1, params["k"]*100.0*4.184)
bondForceBuilder.registerParams((class1, class2), params["r0"], params["k"])
bondForce = bondForceBuilder.getForce(sys)
bondForceBuilder.addBonds(bondForce, atomClasses, bonds)
if "angle" in self._forces or "opbend" in self._forces or "anglep" in self._forces:
# Find all unique angles
bondedToAtom = [[] for _ in range(self.topology.getNumAtoms())]
for (at1, at2) in bonds:
bondedToAtom[at1].append(at2)
bondedToAtom[at2].append(at1)
uniqueAngles = set()
for atom in range(len(bondedToAtom)):
neighbors = bondedToAtom[atom]
......@@ -548,63 +549,18 @@ class TinkerFiles:
uniqueAngles.add(angle)
angles = sorted(list(uniqueAngles))
# Need to initialize here as used to create out-of-plane angles
if "opbend" in self._forces:
opbendParams = {(at1, at2, at3, at4): {"k": float(k)} for at1, at2, at3, at4, k in self._forces["opbend"]}
opbendParams = {(at1, at2, at3 if at3 != '0' else '', at4 if at4 != '0' else ''): {"k": float(k)} for at1, at2, at3, at4, k in self._forces["opbend"]}
else:
opbendParams = {}
# Classify angles into in-plane, out-of-plane, and generic
opbendTypes = list(opbendParams.keys())
inPlaneAngles, outOfPlaneAngles, genericAngles = AmoebaOutOfPlaneBendForceBuilder.classifyAngles(
angles, atomClasses, bondedToAtom, opbendTypes
)
# Classify angles into in-plane, out-of-plane and generic.
# If middle atom has covalency of 3 and the types of the middle atom and the partner atom (atom bonded to
# middle atom, but not in angle) match types1 and types2, then three out-of-plane bend angles are generated.
# Three in-plane angle are also generated.
# If the conditions are not satisfied, the angle is marked as 'generic' angle (not a in-plane angle).
genericAngles = []
inPlaneAngles = []
outOfPlaneAngles = []
skipAtoms = {}
for angle in angles:
middleAtom = angle[1]
middleClass = atomClasses[middleAtom]
middleValence = len(self.atoms[middleAtom].bonds)
if middleValence == 3 and middleAtom not in skipAtoms:
partners = []
for partner in self.atoms[middleAtom].bonds:
partnerClass = atomClasses[partner]
for opBendClass in opbendParams:
if middleClass == opBendClass[1] and partnerClass == opBendClass[0]:
partners.append(partner)
break
if len(partners) == 3:
outOfPlaneAngles.append([partners[0], middleAtom, partners[1], partners[2]])
outOfPlaneAngles.append([partners[2], middleAtom, partners[0], partners[1]])
outOfPlaneAngles.append([partners[1], middleAtom, partners[2], partners[0]])
skipAtoms[middleAtom] = set(partners[:3])
angleList = list(angle[:3])
for atomIndex in partners:
if atomIndex not in angleList:
angleList.append(atomIndex)
inPlaneAngles.append(angleList)
else:
angleList = list(angle[:3])
for atomIndex in partners:
if atomIndex not in angleList:
angleList.append(atomIndex)
genericAngles.append(angleList)
elif middleValence == 3 and middleAtom in skipAtoms:
angleList = list(angle[:3])
for atomIndex in skipAtoms[middleAtom]:
if atomIndex not in angleList:
angleList.append(atomIndex)
inPlaneAngles.append(angleList)
else:
genericAngles.append(list(angle))
# Sanity checks for angle classification
assert len(inPlaneAngles) + len(genericAngles)== len(angles), (
assert len(inPlaneAngles) + len(genericAngles) == len(angles), (
f"Angle classification mismatch:\n"
f" in-plane: {len(inPlaneAngles)}\n"
f" generic: {len(genericAngles)}\n"
......@@ -630,138 +586,52 @@ class TinkerFiles:
angleParams = {(at1, at2, at3): {"k": float(k), "theta0": [float(theta) for theta in theta]} for at1, at2, at3, k, *theta in self._forces["angle"]}
if genericAngles:
angleForceBuilder = AmoebaAngleForceBuilder(self._scalars["angle-cubic"], self._scalars["angle-quartic"], self._scalars["angle-pentic"], self._scalars["angle-sextic"])
processedParams = set()
processedAngles = []
for angle in genericAngles:
class1 = atomClasses[angle[0]]
class2 = atomClasses[angle[1]]
class3 = atomClasses[angle[2]]
params = angleParams.get((class1, class2, class3)) or angleParams.get((class3, class2, class1))
if params:
if len(params["theta0"]) > 1:
# Get k-index by counting number of non-angle hydrogens on the central atom
partners = [at for at in bondedToAtom[angle[1]] if at not in angle]
nHyd = sum(1 for i in partners if self.atoms[i].atomicNumber == 1)
if nHyd < len(params["theta0"]):
theta0 = params["theta0"][nHyd]
else:
raise ValueError(f"Angle parameters out of range for atom classes {class1}-{class2}-{class3}")
else:
theta0 = params["theta0"][0]
angleTuple = (angle[0], angle[1], angle[2])
idealAngles[angleTuple] = theta0 # Store ideal angle for stretch-bend
angleForceBuilder.registerParams(angleTuple, (theta0, params["k"]*4.184*(math.pi/180)**2))
processedAngles.append(angleTuple)
if processedAngles:
angleForce = angleForceBuilder.getForce(sys)
angleForceBuilder.addAngles(angleForce, processedAngles)
angleKey = (atomClasses[angle[0]], atomClasses[angle[1]], atomClasses[angle[2]])
angleTuple = (angle[0], angle[1], angle[2])
params = angleParams.get(angleKey) or angleParams.get(angleKey[::-1])
theta0 = angleForceBuilder.getIdealAngle(angleTuple, params["theta0"], self.atoms, bondedToAtom)
idealAngles[angleTuple] = theta0*math.pi/180 # Store ideal angle for stretch-bend
angleForceBuilder.registerParams(angleTuple, (theta0, params["k"]*4.184*(math.pi/180)**2))
angleForce = angleForceBuilder.getForce(sys)
angleForceBuilder.addAngles(angleForce, genericAngles)
# Add AmoebaInPlaneAngleForce
if "anglep" in self._forces:
inPlaneAngleParams = {(at1, at2, at3): {"k": float(k), "theta0": [float(theta) for theta in theta]} for at1, at2, at3, k, *theta in self._forces["anglep"]}
else:
inPlaneAngleParams = {}
inPlaneAngleParams.update(angleParams) # Poltype e.g. does not use the "anglep" keyword, but uses "angle" keyword for in-plane angles
if inPlaneAngles:
inPlaneAngleForceBuilder = AmoebaInPlaneAngleForceBuilder(self._scalars["angle-cubic"], self._scalars["angle-quartic"], self._scalars["angle-pentic"], self._scalars["angle-sextic"])
processedParams = set()
for angle in inPlaneAngles:
class1 = atomClasses[angle[0]]
class2 = atomClasses[angle[1]]
class3 = atomClasses[angle[2]]
params = inPlaneAngleParams.get((class1, class2, class3)) or inPlaneAngleParams.get((class3, class2, class1))
if params:
class4 = atomClasses[angle[3]]
paramKey = (class1, class2, class3, class4)
if len(params["theta0"]) > 1:
# Get k-index by counting number of non-angle hydrogens on the central atom
partners = [at for at in bondedToAtom[angle[1]] if at not in angle[:]]
nHyd = sum(1 for i in partners if self.atoms[i].atomicNumber == 1)
if nHyd < len(params["theta0"]):
theta0 = params["theta0"][nHyd]
else:
raise ValueError(f"Angle parameters out of range for atom classes {class1}-{class2}-{class3}")
else:
theta0 = params["theta0"][0]
idealAngles[tuple(angle)] = theta0 # Store ideal angle for stretch-bend
if paramKey not in processedParams:
inPlaneAngleForceBuilder.registerParams(paramKey, (theta0, params["k"]*4.184*(math.pi/180)**2))
processedParams.add(paramKey)
if processedParams:
inPlaneAngleForce = inPlaneAngleForceBuilder.getForce(sys)
inPlaneAngleForceBuilder.addInPlaneAngles(inPlaneAngleForce, atomClasses, inPlaneAngles)
angleKey = (atomClasses[angle[0]], atomClasses[angle[1]], atomClasses[angle[2]], atomClasses[angle[3]])
params = inPlaneAngleParams.get(angleKey[:3]) or inPlaneAngleParams.get(angleKey[:3][::-1])
theta0 = angleForceBuilder.getIdealAngle(tuple(angle), params["theta0"], self.atoms, bondedToAtom)
idealAngles[tuple(angle)] = theta0*math.pi/180 # Store ideal angle for stretch-bend
inPlaneAngleForceBuilder.registerParams(angleKey[:3], (theta0, params["k"]*4.184*(math.pi/180)**2))
inPlaneAngleForce = inPlaneAngleForceBuilder.getForce(sys)
inPlaneAngleForceBuilder.addInPlaneAngles(inPlaneAngleForce, atomClasses, inPlaneAngles)
# Add AmoebaStretchBendForce
if "strbnd" in self._forces:
stretchBendParams = {(at1, at2, at3): {"k1": float(k1), "k2": float(k2)} for at1, at2, at3, k1, k2 in self._forces["strbnd"]}
stretchBendParams = {(at1, at2, at3): {"k1": float(k1)*41.84*math.pi/180, "k2": float(k2)*41.84*math.pi/180} for at1, at2, at3, k1, k2 in self._forces["strbnd"]}
stretchBendForceBuilder = AmoebaStretchBendForceBuilder()
processedAngles = []
for angle in genericAngles + inPlaneAngles:
class1 = atomClasses[angle[0]]
class2 = atomClasses[angle[1]]
class3 = atomClasses[angle[2]]
params = stretchBendParams.get((class1, class2, class3))
if params is not None:
swap = False
else:
params = stretchBendParams.get((class3, class2, class1))
swap = params is not None
if params:
idealAngle = idealAngles.get(tuple(angle), None)
if idealAngle is None:
continue
if not swap:
# 1-2-3
bondAB = bondParams.get((class1, class2)) or bondParams.get((class2, class1))
bondCB = bondParams.get((class2, class3)) or bondParams.get((class3, class2))
if bondAB and bondCB:
bondAB, bondCB = bondAB['r0'], bondCB['r0']
k1, k2 = params["k1"], params["k2"]
else:
# 3-2-1
bondAB = bondParams.get((class3, class2)) or bondParams.get((class2, class3))
bondCB = bondParams.get((class2, class1)) or bondParams.get((class1, class2))
if bondAB and bondCB:
bondAB, bondCB = bondCB['r0'], bondAB['r0']
k1, k2 = params["k2"], params["k1"]
if bondAB and bondCB:
# Because for the same class triplet, depending on the number of hydrogens on the central atom,
# the ideal angle can be different, we need to register parameters for each angle.
paramKey = (angle[0], angle[1], angle[2])
stretchBendForceBuilder.registerParams(paramKey, (bondAB*0.1, bondCB*0.1, idealAngle*math.pi/180, k1*41.84*math.pi/180, k2*41.84*math.pi/180))
processedAngles.append((angle[0], angle[1], angle[2]))
if processedAngles:
stretchBendForce = stretchBendForceBuilder.getForce(sys)
stretchBendForceBuilder.addStretchBends(stretchBendForce, processedAngles)
processedAngles = stretchBendForceBuilder.registerAllStretchBendParams(atomClasses, genericAngles + inPlaneAngles, stretchBendParams, bondParams, idealAngles)
stretchBendForce = stretchBendForceBuilder.getForce(sys)
stretchBendForceBuilder.addStretchBends(stretchBendForce, processedAngles)
# Add AMOEBA Urey-Bradley force
# Add AmoebaUreyBradleyForce
if "ureybrad" in self._forces:
ureyBradleyParams = {(at1, at2, at3): {"k": float(k), "d": float(d)} for at1, at2, at3, k, d in self._forces["ureybrad"]}
ureyBradleyForceBuilder = AmoebaUreyBradleyForceBuilder()
for (class1, class2, class3), params in ureyBradleyParams.items():
ureyBradleyForceBuilder.registerParams((class1, class2, class3), (params["k"]*4.184*100.0, params["d"]*0.1))
processedAngles = []
for at1, at2, at3 in angles:
class1 = atomClasses[at1]
class2 = atomClasses[at2]
class3 = atomClasses[at3]
params = ureyBradleyParams.get((class1, class2, class3)) or ureyBradleyParams.get((class3, class2, class1))
if params:
processedAngles.append((at1, at2, at3))
if processedAngles:
ureyBradleyForce = ureyBradleyForceBuilder.getForce(sys)
ureyBradleyForceBuilder.addUreyBradleys(ureyBradleyForce, atomClasses, processedAngles)
# Find all unique proper torsions in the system
ureyBradleyForceBuilder.registerParams((class1, class2, class3), (params["k"]*4.184*100.0, params["d"]*0.1), isClass=True)
ureyBradleyForce = ureyBradleyForceBuilder.getForce(sys)
ureyBradleyForceBuilder.addUreyBradleys(ureyBradleyForce, atomClasses, angles, isClass=True)
# Find all unique proper torsions
uniquePropers = set()
for angle in angles:
for atom in self.atoms[angle[0]].bonds:
......@@ -784,58 +654,23 @@ class TinkerFiles:
torsionParams = {(at1, at2, at3, at4): {"t1": [float(t11)*torsionScale, float(t12)*math.pi/180.0, int(t13)],
"t2": [float(t21)*torsionScale, float(t22)*math.pi/180.0, int(t23)],
"t3": [float(t31)*torsionScale, float(t32)*math.pi/180.0, int(t33)]}
for at1, at2, at3, at4, t11, t12, t13, t21, t22, t23, t31, t32, t33 in self._forces["torsion"]}
for at1, at2, at3, at4, t11, t12, t13, t21, t22, t23, t31, t32, t33 in self._forces["torsion"]}
torsionForceBuilder = AmoebaTorsionForceBuilder()
for (class1, class2, class3, class4), params in torsionParams.items():
torsionForceBuilder.registerParams((class1, class2, class3, class4), (params["t1"], params["t2"], params["t3"]))
processedTorsions = []
for at1, at2, at3, at4 in propers:
class1 = atomClasses[at1]
class2 = atomClasses[at2]
class3 = atomClasses[at3]
class4 = atomClasses[at4]
params = torsionParams.get((class1, class2, class3, class4)) or torsionParams.get((class4, class3, class2, class1))
if params:
processedTorsions.append((at1, at2, at3, at4))
if processedTorsions:
torsionForce = torsionForceBuilder.getForce(sys)
atomClassesDict = {}
for atom_idx in range(len(self.atoms)):
atomClassesDict[atom_idx] = atomClasses[atom_idx]
torsionForceBuilder.addTorsions(torsionForce, atomClasses, processedTorsions)
torsionForce = torsionForceBuilder.getForce(sys)
torsionForceBuilder.addTorsions(torsionForce, atomClasses, propers)
# Add AmoebaPiTorsionForce
if "pitors" in self._forces:
piTorsionParams = {(at1, at2): {"k": float(k)} for at1, at2, k in self._forces["pitors"]}
piTorsionForceBuilder = AmoebaPiTorsionForceBuilder()
piTorsionScale = float(self._scalars["pitorsunit"])*4.184
piTorsionParams = {(at1, at2): {"k": float(k)*piTorsionScale} for at1, at2, k in self._forces["pitors"]}
piTorsionForceBuilder = AmoebaPiTorsionForceBuilder()
for (class1, class2), params in piTorsionParams.items():
piTorsionForceBuilder.registerParams((class1, class2), (params["k"]*piTorsionScale,))
processedPiTorsions = []
for (at1, at2) in bonds:
valence1 = len(self.atoms[at1].bonds)
valence2 = len(self.atoms[at2].bonds)
if valence1 == 3 and valence2 == 3:
class1 = atomClasses[at1]
class2 = atomClasses[at2]
params = piTorsionParams.get((class1, class2)) or piTorsionParams.get((class2, class1))
if params:
piTorsionAtom3 = at1
piTorsionAtom4 = at2
# piTorsionAtom1, piTorsionAtom2 are the atoms bonded to atom1, excluding atom2
# piTorsionAtom5, piTorsionAtom6 are the atoms bonded to atom2, excluding atom1
piTorsionAtom1, piTorsionAtom2 = [bond for bond in self.atoms[at1].bonds if bond != piTorsionAtom4]
piTorsionAtom5, piTorsionAtom6 = [bond for bond in self.atoms[at2].bonds if bond != piTorsionAtom3]
processedPiTorsions.append((piTorsionAtom1, piTorsionAtom2, piTorsionAtom3, piTorsionAtom4, piTorsionAtom5, piTorsionAtom6))
if processedPiTorsions:
piTorsionForce = piTorsionForceBuilder.getForce(sys)
piTorsionForceBuilder.addPiTorsions(piTorsionForce, atomClassesDict, processedPiTorsions)
piTorsionForceBuilder.registerParams((class1, class2), (params["k"],))
processedPiTorsions = piTorsionForceBuilder.getAllPiTorsions(atomClasses, bondedToAtom, bonds)
piTorsionForce = piTorsionForceBuilder.getForce(sys)
piTorsionForceBuilder.addPiTorsions(piTorsionForce, atomClasses, processedPiTorsions)
# Add AmoebaStretchTorsionForce
if "strtors" in self._forces:
......@@ -844,18 +679,8 @@ class TinkerFiles:
stretchTorsionForceBuilder = AmoebaStretchTorsionForceBuilder()
for classes, params in stretchTorsionParams.items():
stretchTorsionForceBuilder.registerParams(classes, params)
processedStretchTorsions = []
for at1, at2, at3, at4 in propers:
class1 = atomClasses[at1]
class2 = atomClasses[at2]
class3 = atomClasses[at3]
class4 = atomClasses[at4]
params = stretchTorsionParams.get((class1, class2, class3, class4)) or stretchTorsionParams.get((class4, class3, class2, class1))
if params:
processedStretchTorsions.append((at1, at2, at3, at4))
if processedStretchTorsions:
stretchTorsionForce = stretchTorsionForceBuilder.getForce(sys)
stretchTorsionForceBuilder.addStretchTorsions(sys, stretchTorsionForce, atomClasses, processedStretchTorsions)
stretchTorsionForce = stretchTorsionForceBuilder.getForce(sys)
stretchTorsionForceBuilder.addStretchTorsions(sys, stretchTorsionForce, atomClasses, propers)
# Add AmoebaAngleTorsionForce
if "angtors" in self._forces:
......@@ -864,36 +689,28 @@ class TinkerFiles:
angleTorsionForceBuilder = AmoebaAngleTorsionForceBuilder()
for classes, params in angleTorsionParams.items():
angleTorsionForceBuilder.registerParams(classes, params)
processedAngleTorsions = []
for at1, at2, at3, at4 in propers:
class1 = atomClasses[at1]
class2 = atomClasses[at2]
class3 = atomClasses[at3]
class4 = atomClasses[at4]
params = angleTorsionParams.get((class1, class2, class3, class4)) or angleTorsionParams.get((class4, class3, class2, class1))
if params:
processedAngleTorsions.append((at1, at2, at3, at4))
if processedAngleTorsions:
angleTorsionForce = angleTorsionForceBuilder.getForce(sys)
angleTorsionForceBuilder.addAngleTorsions(sys, angleTorsionForce, atomClasses, processedAngleTorsions)
angleTorsionForce = angleTorsionForceBuilder.getForce(sys)
angleTorsionForceBuilder.addAngleTorsions(sys, angleTorsionForce, atomClasses, propers)
# Add AmoebaTorsionTorsionForce
if "tortors" in self._forces:
torsionTorsionForce = AmoebaTorsionTorsionForceBuilder()
torsionTorsionForceBuilder = AmoebaTorsionTorsionForceBuilder()
for gridIndex, (tortorInfo, gridData) in enumerate(self._forces["tortors"]):
tortorType = (tortorInfo[0], tortorInfo[1], tortorInfo[2], tortorInfo[3], tortorInfo[4])
torsionTorsionForceBuilder.registerParams(tortorType, gridIndex)
nx = int(tortorInfo[5])
ny = int(tortorInfo[6])
grid = np.array(gridData, dtype=np.float64).reshape((nx, ny, -1))
grid[:, :, 2] *= 4.184
torsionTorsionForce.registerGridData(gridIndex, grid)
force = torsionTorsionForce.getForce(sys)
AmoebaTorsionTorsionForceBuilder.createTorsionTorsionInteractions(force, angles, self.atoms, self._forces["tortors"])
torsionTorsionForceBuilder.registerGridData(gridIndex, grid)
torsionTorsionForce = torsionTorsionForceBuilder.getForce(sys)
torsionTorsionForceBuilder.addTorsionTorsionInteractions(torsionTorsionForce, propers, atomClasses, [atom.mass for atom in self.atoms])
# Add AmoebaVdwForce
if "vdw" in self._forces:
vdwForceBuilder = AmoebaVdwForceBuilder(self._scalars['vdwtype'], self._scalars['radiusrule'], self._scalars['radiustype'],
self._scalars['radiussize'], self._scalars['epsilonrule'], float(self._scalars['vdw-13-scale']),
float(self._scalars['vdw-14-scale']), float(self._scalars['vdw-15-scale']))
self._scalars['radiussize'], self._scalars['epsilonrule'], float(self._scalars['vdw-13-scale']),
float(self._scalars['vdw-14-scale']), float(self._scalars['vdw-15-scale']))
for params in self._forces['vdw']:
reduction = 0.0 if len(params) < 4 else float(params[3])
vdwForceBuilder.registerClassParams(int(params[0]), float(params[1])*0.1, float(params[2])*4.184, reduction)
......@@ -923,7 +740,6 @@ class TinkerFiles:
# Add AmoebaGeneralizedKirkwoodForce
if implicitSolvent and "multipole" in self._forces:
gkForceBuilder = AmoebaGeneralizedKirkwoodForceBuilder(**self.GK_PARAMS)
force = gkForceBuilder.getForce(sys, implicitSolvent=implicitSolvent)
multipoleForce = [f for f in sys.getForces() if isinstance(f, mm.AmoebaMultipoleForce)][0]
for atomIndex in range(0, multipoleForce.getNumMultipoles()):
multipoleParameters = multipoleForce.getMultipoleParameters(atomIndex)
......@@ -939,15 +755,12 @@ class TinkerFiles:
convert *= 1.122462048309372
if self._scalars["radiussize"] == "DIAMETER":
convert *= 0.5
# Register VdW parameters for each atom class
for vdw in self._forces["vdw"]:
atomClass = int(vdw[0])
sigma = float(vdw[1])*convert
epsilon = float(vdw[2])*4.184
wcaDispersionForceBuilder.registerParams(atomClass, sigma, epsilon)
wcaDispersionForce = wcaDispersionForceBuilder.getForce(sys)
atomClasses = [int(atom.atomClass) for atom in self.atoms]
wcaDispersionForceBuilder.addParticles(wcaDispersionForce, atomClasses)
# Set periodic boundary conditions
......@@ -973,7 +786,6 @@ class TinkerFiles:
return sys
# ------------------------------------------------------------------------------------------ #
# TOPOLOGY FUNCTIONS #
# ------------------------------------------------------------------------------------------ #
......
wrappers/python/tests/TestForceField.py
View file @ 489e2c46
......@@ -1958,12 +1958,12 @@ class AmoebaTestForceField(unittest.TestCase):
# Compare to values computed with Tinker.
self.assertAlmostEqual(290.2445, energies['AmoebaBond'], 4)
self.assertAlmostEqual(496.4300, energies['AmoebaAngle']+energies['AmoebaInPlaneAngle'], 4)
self.assertAlmostEqual(51.2913, energies['AmoebaOutOfPlaneBend'], 4)
self.assertAlmostEqual(5.7695, energies['AmoebaStretchBend'], 4)
self.assertAlmostEqual(290.2445, energies['AmoebaBondForce'], 4)
self.assertAlmostEqual(496.4300, energies['AmoebaAngleForce']+energies['AmoebaInPlaneAngleForce'], 4)
self.assertAlmostEqual(51.2913, energies['AmoebaOutOfPlaneBendForce'], 4)
self.assertAlmostEqual(5.7695, energies['AmoebaStretchBendForce'], 4)
self.assertAlmostEqual(75.6890, energies['PeriodicTorsionForce'], 4)
self.assertAlmostEqual(19.3364, energies['AmoebaPiTorsion'], 4)
self.assertAlmostEqual(19.3364, energies['AmoebaPiTorsionForce'], 4)
self.assertAlmostEqual(-32.6689, energies['AmoebaTorsionTorsionForce'], 4)
self.assertAlmostEqual(383.8705, energies['AmoebaVdwForce'], 4)
self.assertAlmostEqual(-1323.5640-225.3660, energies['AmoebaMultipoleForce'], 2)
......@@ -1975,14 +1975,14 @@ class AmoebaTestForceField(unittest.TestCase):
# Compare to values computed with Tinker.
self.assertAlmostEqual(749.6953, energies['AmoebaBond'], 4)
self.assertAlmostEqual(579.9971, energies['AmoebaAngle']+energies['AmoebaInPlaneAngle'], 4)
self.assertAlmostEqual(10.6630, energies['AmoebaOutOfPlaneBend'], 4)
self.assertAlmostEqual(5.2225, energies['AmoebaStretchBend'], 4)
self.assertAlmostEqual(749.6953, energies['AmoebaBondForce'], 4)
self.assertAlmostEqual(579.9971, energies['AmoebaAngleForce']+energies['AmoebaInPlaneAngleForce'], 4)
self.assertAlmostEqual(10.6630, energies['AmoebaOutOfPlaneBendForce'], 4)
self.assertAlmostEqual(5.2225, energies['AmoebaStretchBendForce'], 4)
self.assertAlmostEqual(166.7233, energies['PeriodicTorsionForce'], 4)
self.assertAlmostEqual(57.2066, energies['AmoebaPiTorsion'], 4)
self.assertAlmostEqual(-4.2538, energies['AmoebaStretchTorsion'], 4)
self.assertAlmostEqual(-5.0402, energies['AmoebaAngleTorsion'], 4)
self.assertAlmostEqual(57.2066, energies['AmoebaPiTorsionForce'], 4)
self.assertAlmostEqual(-4.2538, energies['AmoebaStretchTorsionForce'], 4)
self.assertAlmostEqual(-5.0402, energies['AmoebaAngleTorsionForce'], 4)
self.assertAlmostEqual(187.1103, energies['AmoebaVdwForce'], 4)
self.assertAlmostEqual(1635.1289-236.1484, energies['AmoebaMultipoleForce'], 3)
self.assertAlmostEqual(3146.3046, sum(list(energies.values())), 3)
......
wrappers/python/tests/TestTinkerFiles.py
View file @ 489e2c46
......@@ -14,8 +14,6 @@ class TestTinkerFiles(unittest.TestCase):
polarization="mutual",
mutualInducedTargetEpsilon=1e-5,
nonbondedMethod=NoCutoff,
constraints=None,
useDispersionCorrection=False,
)
# Compute the energy with OpenMM.
......@@ -91,10 +89,10 @@ class TestTinkerFiles(unittest.TestCase):
energies, _, _ = self.computeAmoebaEnergies(xyzFile, keyFiles)
# Compare to values computed with Tinker.
self.assertEnergyEqual(1104.0455, energies["AmoebaBond"])
self.assertEnergyEqual(602.7082, energies["AmoebaAngle"] + energies["AmoebaInPlaneAngle"])
self.assertEnergyEqual(2.0572, energies["AmoebaOutOfPlaneBend"], 1e-4)
self.assertEnergyEqual(-0.1361, energies["AmoebaStretchBend"], 1e-3)
self.assertEnergyEqual(1104.0455, energies["AmoebaBondForce"])
self.assertEnergyEqual(602.7082, energies["AmoebaAngleForce"] + energies["AmoebaInPlaneAngleForce"])
self.assertEnergyEqual(2.0572, energies["AmoebaOutOfPlaneBendForce"], 1e-4)
self.assertEnergyEqual(-0.1361, energies["AmoebaStretchBendForce"], 1e-3)
self.assertEnergyEqual(-0.8625, energies["PeriodicTorsionForce"], 1e-4)
self.assertEnergyEqual(-33.8595, energies["HarmonicBondForce"])
self.assertEnergyEqual(5908.1343, energies["AmoebaVdwForce"])
......@@ -156,12 +154,12 @@ class TestTinkerFiles(unittest.TestCase):
'LYS', 'ARG'], f'Unexpected residues: {residues}'
# Compare to values computed with Tinker.
self.assertEnergyEqual(19.6519, energies["AmoebaBond"])
self.assertEnergyEqual(58.2509, energies["AmoebaAngle"] + energies["AmoebaInPlaneAngle"])
self.assertEnergyEqual(1.9697, energies["AmoebaOutOfPlaneBend"], 1e-4)
self.assertEnergyEqual(-0.4384, energies["AmoebaStretchBend"], 1e-3)
self.assertEnergyEqual(19.6519, energies["AmoebaBondForce"])
self.assertEnergyEqual(58.2509, energies["AmoebaAngleForce"] + energies["AmoebaInPlaneAngleForce"])
self.assertEnergyEqual(1.9697, energies["AmoebaOutOfPlaneBendForce"], 1e-4)
self.assertEnergyEqual(-0.4384, energies["AmoebaStretchBendForce"], 1e-3)
self.assertEnergyEqual(-2.3514, energies["PeriodicTorsionForce"], 1e-4)
self.assertEnergyEqual(1.2115, energies["AmoebaPiTorsion"], 1e-4)
self.assertEnergyEqual(1.2115, energies["AmoebaPiTorsionForce"], 1e-4)
self.assertEnergyEqual(-3.2958, energies["AmoebaTorsionTorsionForce"])
self.assertEnergyEqual(1509.1915, energies["AmoebaVdwForce"])
self.assertEnergyEqual(-488.0403 - 110.9042, energies["AmoebaMultipoleForce"], 1e-3)
......@@ -224,14 +222,14 @@ class TestTinkerFiles(unittest.TestCase):
'G', 'C', 'G', 'U', 'U', 'A', 'A', 'G', 'U', 'C', 'G', 'C', 'A'], f'Unexpected residues: {residues}'
# Compare to values computed with Tinker.
self.assertEnergyEqual(749.6953, energies["AmoebaBond"])
self.assertEnergyEqual(579.9971, energies["AmoebaAngle"] + energies["AmoebaInPlaneAngle"])
self.assertEnergyEqual(10.6630, energies["AmoebaOutOfPlaneBend"])
self.assertEnergyEqual(5.2225, energies["AmoebaStretchBend"])
self.assertEnergyEqual(749.6953, energies["AmoebaBondForce"])
self.assertEnergyEqual(579.9971, energies["AmoebaAngleForce"] + energies["AmoebaInPlaneAngleForce"])
self.assertEnergyEqual(10.6630, energies["AmoebaOutOfPlaneBendForce"])
self.assertEnergyEqual(5.2225, energies["AmoebaStretchBendForce"])
self.assertEnergyEqual(166.7233, energies["PeriodicTorsionForce"])
self.assertEnergyEqual(57.2066, energies["AmoebaPiTorsion"])
self.assertEnergyEqual(-4.2538, energies["AmoebaStretchTorsion"])
self.assertEnergyEqual(-0.0880, energies["AmoebaAngleTorsion"], 1e-4)
self.assertEnergyEqual(57.2066, energies["AmoebaPiTorsionForce"])
self.assertEnergyEqual(-4.2538, energies["AmoebaStretchTorsionForce"])
self.assertEnergyEqual(-0.0880, energies["AmoebaAngleTorsionForce"], 1e-4)
self.assertEnergyEqual(187.1103, energies["AmoebaVdwForce"])
self.assertEnergyEqual(1635.1289 - 236.1484, energies["AmoebaMultipoleForce"])
self.assertEnergyEqual(3151.2568, sum(list(energies.values())))
......
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