Skip to content

GitLab

Commit cd4031f9 authored by Robert T. McGibbon's avatar Robert T. McGibbon
Browse files

Merge branch 'master' of github.com:pandegroup/openmm into sphinx

parents 9af3e56b 4d32047c
Hide whitespace changes
Inline Side-by-side
Showing with 103 additions and 138 deletions
wrappers/python/simtk/openmm/app/data/tip4pew.xml
View file @ cd4031f9
......@@ -10,9 +10,9 @@
<Atom name="H1" type="tip4pew-H"/>
<Atom name="H2" type="tip4pew-H"/>
<Atom name="M" type="tip4pew-M"/>
<VirtualSite type="average3" index="3" atom1="0" atom2="1" atom3="2" weight1="0.786646558" weight2="0.106676721" weight3="0.106676721"/>
<Bond from="0" to="1"/>
<Bond from="0" to="2"/>
<VirtualSite type="average3" siteName="M" atomName1="O" atomName2="H1" atomName3="H2" weight1="0.786646558" weight2="0.106676721" weight3="0.106676721"/>
<Bond atomName1="O" atomName2="H1"/>
<Bond atomName1="O" atomName2="H2"/>
</Residue>
</Residues>
<HarmonicBondForce>
......
wrappers/python/simtk/openmm/app/data/tip5p.xml
View file @ cd4031f9
......@@ -11,10 +11,10 @@
<Atom name="H2" type="tip5p-H"/>
<Atom name="M1" type="tip5p-M"/>
<Atom name="M2" type="tip5p-M"/>
<VirtualSite type="outOfPlane" index="3" atom1="0" atom2="1" atom3="2" weight12="-0.34490826" weight13="-0.34490826" weightCross="-6.4437903"/>
<VirtualSite type="outOfPlane" index="4" atom1="0" atom2="1" atom3="2" weight12="-0.34490826" weight13="-0.34490826" weightCross="6.4437903"/>
<Bond from="0" to="1"/>
<Bond from="0" to="2"/>
<VirtualSite type="outOfPlane" siteName="M1" atomName1="O" atomName2="H1" atomName3="H2" weight12="-0.34490826" weight13="-0.34490826" weightCross="-6.4437903"/>
<VirtualSite type="outOfPlane" siteName="M2" atomName1="O" atomName2="H1" atomName3="H2" weight12="-0.34490826" weight13="-0.34490826" weightCross="6.4437903"/>
<Bond atomName1="O" atomName2="H1"/>
<Bond atomName1="O" atomName2="H2"/>
</Residue>
</Residues>
<HarmonicBondForce>
......
wrappers/python/simtk/openmm/app/forcefield.py
View file @ cd4031f9
......@@ -153,18 +153,29 @@ class ForceField(object):
for residue in root.find('Residues').findall('Residue'):
resName = residue.attrib['name']
template = ForceField._TemplateData(resName)
atomIndices = {}
for atom in residue.findall('Atom'):
params = {}
for key in atom.attrib:
if key not in ('name', 'type'):
params[key] = _convertParameterToNumber(atom.attrib[key])
template.atoms.append(ForceField._TemplateAtomData(atom.attrib['name'], atom.attrib['type'], self._atomTypes[atom.attrib['type']][2], params))
atomName = atom.attrib['name']
if atomName in atomIndices:
raise ValueError('Residue '+resName+' contains multiple atoms named '+atomName)
atomIndices[atomName] = len(template.atoms)
template.atoms.append(ForceField._TemplateAtomData(atomName, atom.attrib['type'], self._atomTypes[atom.attrib['type']][2], params))
for site in residue.findall('VirtualSite'):
template.virtualSites.append(ForceField._VirtualSiteData(site))
template.virtualSites.append(ForceField._VirtualSiteData(site, atomIndices))
for bond in residue.findall('Bond'):
template.addBond(int(bond.attrib['from']), int(bond.attrib['to']))
if 'atomName1' in bond.attrib:
template.addBond(atomIndices[bond.attrib['atomName1']], atomIndices[bond.attrib['atomName2']])
else:
template.addBond(int(bond.attrib['from']), int(bond.attrib['to']))
for bond in residue.findall('ExternalBond'):
b = int(bond.attrib['from'])
if 'atomName' in bond.attrib:
b = atomIndices[bond.attrib['atomName']]
else:
b = int(bond.attrib['from'])
template.externalBonds.append(b)
template.atoms[b].externalBonds += 1
self.registerResidueTemplate(template)
......@@ -320,27 +331,32 @@ class ForceField(object):
class _VirtualSiteData:
"""Inner class used to encapsulate data about a virtual site."""
def __init__(self, node):
def __init__(self, node, atomIndices):
attrib = node.attrib
self.index = int(attrib['index'])
self.type = attrib['type']
if self.type == 'average2':
self.atoms = [int(attrib['atom1']), int(attrib['atom2'])]
numAtoms = 2
self.weights = [float(attrib['weight1']), float(attrib['weight2'])]
elif self.type == 'average3':
self.atoms = [int(attrib['atom1']), int(attrib['atom2']), int(attrib['atom3'])]
numAtoms = 3
self.weights = [float(attrib['weight1']), float(attrib['weight2']), float(attrib['weight3'])]
elif self.type == 'outOfPlane':
self.atoms = [int(attrib['atom1']), int(attrib['atom2']), int(attrib['atom3'])]
numAtoms = 3
self.weights = [float(attrib['weight12']), float(attrib['weight13']), float(attrib['weightCross'])]
elif self.type == 'localCoords':
self.atoms = [int(attrib['atom1']), int(attrib['atom2']), int(attrib['atom3'])]
numAtoms = 3
self.originWeights = [float(attrib['wo1']), float(attrib['wo2']), float(attrib['wo3'])]
self.xWeights = [float(attrib['wx1']), float(attrib['wx2']), float(attrib['wx3'])]
self.yWeights = [float(attrib['wy1']), float(attrib['wy2']), float(attrib['wy3'])]
self.localPos = [float(attrib['p1']), float(attrib['p2']), float(attrib['p3'])]
else:
raise ValueError('Unknown virtual site type: %s' % self.type)
if 'siteName' in attrib:
self.index = atomIndices[attrib['siteName']]
self.atoms = [atomIndices[attrib['atomName%d'%(i+1)]] for i in range(numAtoms)]
else:
self.index = int(attrib['index'])
self.atoms = [int(attrib['atom%d'%(i+1)]) for i in range(numAtoms)]
if 'excludeWith' in attrib:
self.excludeWith = int(attrib['excludeWith'])
else:
......@@ -1377,92 +1393,6 @@ class GBSAOBCGenerator:
parsers["GBSAOBCForce"] = GBSAOBCGenerator.parseElement
## @private
class GBVIGenerator:
"""A GBVIGenerator constructs a GBVIForce."""
def __init__(self,ff):
self.ff = ff
self.fixedParameters = {}
self.fixedParameters['soluteDielectric'] = 1.0
self.fixedParameters['solventDielectric'] = 78.3
self.fixedParameters['scalingMethod'] = 1
self.fixedParameters['quinticUpperBornRadiusLimit'] = 5.0
self.fixedParameters['quinticLowerLimitFactor'] = 0.8
self.typeMap = {}
@staticmethod
def parseElement(element, ff):
generator = GBVIGenerator(ff)
for key in generator.fixedParameters.iterkeys():
if (key in element.attrib):
generator.fixedParameters[key] = float(element.attrib[key])
ff.registerGenerator(generator)
for atom in element.findall('Atom'):
types = ff._findAtomTypes(atom.attrib, 1)
if None not in types:
values = (float(atom.attrib['charge']), float(atom.attrib['radius']), float(atom.attrib['gamma']))
for t in types[0]:
generator.typeMap[t] = values
def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
methodMap = {NoCutoff:mm.NonbondedForce.NoCutoff,
CutoffNonPeriodic:mm.NonbondedForce.CutoffNonPeriodic,
CutoffPeriodic:mm.NonbondedForce.CutoffPeriodic}
if nonbondedMethod not in methodMap:
raise ValueError('Illegal nonbonded method for GB/VI Force')
# add particles
force = mm.GBVIForce()
for atom in data.atoms:
t = data.atomType[atom]
if t in self.typeMap:
values = self.typeMap[t]
force.addParticle(values[0], values[1], values[2])
else:
raise ValueError('No GB/VI parameters defined for atom type '+t)
# get HarmonicBond generator -- exit if not found
hbGenerator = 0
for generator in self.ff._forces:
if (generator.__class__.__name__ == 'HarmonicBondGenerator'):
hbGenerator = generator
break
if (hbGenerator == 0):
raise ValueError('HarmonicBondGenerator not found.')
# add bonds
for bond in data.bonds:
type1 = data.atomType[data.atoms[bond.atom1]]
type2 = data.atomType[data.atoms[bond.atom2]]
for i in range(len(hbGenerator.types1)):
types1 = hbGenerator.types1[i]
types2 = hbGenerator.types2[i]
if (type1 in types1 and type2 in types2) or (type1 in types2 and type2 in types1):
#bond.length = hbGenerator.length[i]
force.addBond(bond.atom1, bond.atom2, hbGenerator.length[i])
force.setNonbondedMethod(methodMap[nonbondedMethod])
force.setCutoffDistance(nonbondedCutoff)
force.setSolventDielectric(self.fixedParameters['solventDielectric'])
force.setSoluteDielectric(self.fixedParameters['soluteDielectric'])
force.setBornRadiusScalingMethod(self.fixedParameters['scalingMethod'])
force.setQuinticLowerLimitFactor(self.fixedParameters['quinticLowerLimitFactor'])
force.setQuinticUpperBornRadiusLimit(self.fixedParameters['quinticUpperBornRadiusLimit'])
sys.addForce(force)
parsers["GBVIForce"] = GBVIGenerator.parseElement
## @private
class CustomBondGenerator:
"""A CustomBondGenerator constructs a CustomBondForce."""
......@@ -2695,7 +2625,7 @@ class AmoebaPiTorsionGenerator:
atom1 = bond.atom1
atom2 = bond.atom2
if (len(data.atomBonds[atom1]) == 3 and len(data.atomBonds[atom1]) == 3):
if (len(data.atomBonds[atom1]) == 3 and len(data.atomBonds[atom2]) == 3):
type1 = data.atomType[data.atoms[atom1]]
type2 = data.atomType[data.atoms[atom2]]
......@@ -2980,7 +2910,7 @@ class AmoebaTorsionTorsionGenerator:
# match in reverse order
if (type5 in types1 and type4 in types2 and type3 in types3 and type2 in types4 and type1 in types5):
elif (type5 in types1 and type4 in types2 and type3 in types3 and type2 in types4 and type1 in types5):
chiralAtomIndex = self.getChiralAtomIndex(data, sys, ib, ic, id)
force.addTorsionTorsion(ie, id, ic, ib, ia, chiralAtomIndex, self.gridIndex[i])
......
wrappers/python/simtk/openmm/app/internal/amber_file_parser.py
View file @ cd4031f9
......@@ -674,7 +674,7 @@ def readAmberSystem(topology, prmtop_filename=None, prmtop_loader=None, shake=No
OPTIONAL ARGUMENTS
shake (String) - if 'h-bonds', will SHAKE all bonds to hydrogen and water; if 'all-bonds', will SHAKE all bonds and water (default: None)
gbmodel (String) - if 'OBC', OBC GBSA will be used; if 'GBVI', GB/VI will be used (default: None)
gbmodel (String) - if 'OBC', OBC GBSA will be used (default: None)
soluteDielectric (float) - The solute dielectric constant to use in the implicit solvent model (default: 1.0)
solventDielectric (float) - The solvent dielectric constant to use in the implicit solvent model (default: 78.5)
implicitSolventKappa (float) - The Debye screening parameter corresponding to implicit solvent ionic strength
......@@ -1097,8 +1097,6 @@ def readAmberSystem(topology, prmtop_filename=None, prmtop_loader=None, shake=No
# created above. Do not bind force to another name before this!
force.setReactionFieldDielectric(1.0)
# TODO: Add GBVI terms?
return system
#=============================================================================================
......
wrappers/python/simtk/openmm/app/modeller.py
View file @ cd4031f9
......@@ -1109,7 +1109,7 @@ class Modeller(object):
# Now run a few iterations of SHAKE to try to select reasonable positions.
for iteration in range(10):
for iteration in range(15):
for atom1, atom2, distance in bonds:
if atom1 in missingPositions:
if atom2 in missingPositions:
......
wrappers/python/simtk/openmm/app/pdbfile.py
View file @ cd4031f9
......@@ -154,7 +154,7 @@ class PDBFile(object):
# Add bonds based on CONECT records.
connectBonds = []
for connect in pdb.models[0].connects:
for connect in pdb.models[-1].connects:
i = connect[0]
for j in connect[1:]:
if i in atomByNumber and j in atomByNumber:
......
wrappers/python/src/swig_doxygen/swigInputBuilder.py
View file @ cd4031f9
......@@ -378,7 +378,6 @@ class SwigInputBuilder:
self.fOut.write("\n%s};\n" % INDENT)
if len(enumNodes)>0: self.fOut.write("\n")
def writeMethods(self, classNode):
methodList=getClassMethodList(classNode, self.skipMethods)
......@@ -465,35 +464,44 @@ class SwigInputBuilder:
(shortClassName, memberNode,
shortMethDefinition, methName,
isConstructors, isDestructor, templateType, templateName) = items
paramList=findNodes(memberNode, 'param')
paramList = findNodes(memberNode, 'param')
#write pythonprepend blocks
# write pythonprepend blocks
mArgsstring = getText("argsstring", memberNode)
if self.fOutPythonprepend and \
len(paramList) and \
mArgsstring.find('=0')<0:
key=(shortClassName, methName)
if key in self.configModule.STEAL_OWNERSHIP:
for argNum in self.configModule.STEAL_OWNERSHIP[key]:
if self.SWIG_COMPACT_ARGUMENTS:
argName = 'args[%s]' % argNum
else:
argName = getText('declname', paramList[argNum])
mArgsstring.find('=0') < 0:
text = '''
%pythonprepend OpenMM::{shortClassName}::{methName}{mArgsstring} %{{{{{{0}}
%}}}}'''.format(shortClassName=shortClassName, methName=methName, mArgsstring=mArgsstring)
textInside = ''
key = (shortClassName, methName)
for argNum in self.configModule.STEAL_OWNERSHIP.get(key, []):
if self.SWIG_COMPACT_ARGUMENTS:
argName = 'args[%s]' % argNum
else:
argName = getText('declname', paramList[argNum])
text = '''
%pythonprepend OpenMM::{shortClassName}::{methName}{mArgsstring} %{{
textInside += '''
if not {argName}.thisown:
s = ("the %s object does not own its corresponding OpenMM object"
% self.__class__.__name__)
raise Exception(s)
%}}'''.format(argName=argName, shortClassName=shortClassName, methName=methName, mArgsstring=mArgsstring)
self.fOutPythonprepend.write(text)
raise Exception(s)'''.format(argName=argName)
for argNum in self.configModule.REQUIRE_ORDERED_SET.get(key, []):
if self.SWIG_COMPACT_ARGUMENTS:
argName = 'args[%s]' % argNum
else:
argName = getText('declname', paramList[argNum])
textInside += '''
{argName} = list({argName})'''.format(argName=argName)
if textInside:
self.fOutPythonprepend.write(text.format(textInside))
#write pythonappend blocks
# write pythonappend blocks
if self.fOutPythonappend \
and mArgsstring.find('=0')<0:
key=(shortClassName, methName)
and mArgsstring.find('=0') < 0:
key = (shortClassName, methName)
#print "key %s %s \n" % (shortClassName, methName)
addText=''
returnType = getText("type", memberNode)
......
wrappers/python/src/swig_doxygen/swigInputConfig.py
View file @ cd4031f9
......@@ -43,7 +43,6 @@ SKIP_METHODS = [('State',),
('CalcCustomTorsionForceKernel',),
('CalcForcesAndEnergyKernel',),
('CalcGBSAOBCForceKernel',),
('CalcGBVIForceKernel',),
('CalcHarmonicAngleForceKernel',),
('CalcHarmonicBondForceKernel',),
('CalcKineticEnergyKernel',),
......@@ -144,6 +143,12 @@ STEAL_OWNERSHIP = {("Platform", "registerPlatform") : [0],
("CompoundIntegrator", "addIntegrator") : [0],
}
REQUIRE_ORDERED_SET = {("CustomNonbondedForce", "addInteractionGroup") : [0, 1],
("CustomNonbondedForce", "setInteractionGroupParameters") : [1, 2],
}
# This is a list of units to attach to return values and method args.
# Indexed by (ClassName, MethodsName)
UNITS = {
......@@ -387,14 +392,6 @@ UNITS = {
: (None, ('unit.elementary_charge',
'unit.nanometer', None)),
("GBSAOBCForce", "getSurfaceAreaEnergy") : ('unit.kilojoule_per_mole/unit.nanometer/unit.nanometer', ()),
("GBVIForce", "getBornRadiusScalingMethod") : (None, ()),
("GBVIForce", "getQuinticLowerLimitFactor") : (None, ()),
("GBVIForce", "getQuinticUpperBornRadiusLimit") : ('unit.nanometer', ()),
("GBVIForce", "getBondParameters")
: (None, (None, None, 'unit.nanometer')),
("GBVIForce", "getParticleParameters")
: (None, ('unit.elementary_charge',
'unit.nanometer', 'unit.kilojoule_per_mole')),
("HarmonicAngleForce", "getAngleParameters")
: (None, (None, None, None, 'unit.radian',
'unit.kilojoule_per_mole/(unit.radian*unit.radian)')),
......
wrappers/python/tests/TestAmberPrmtopFile.py
View file @ cd4031f9
......@@ -296,6 +296,38 @@ class TestAmberPrmtopFile(unittest.TestCase):
diff = norm(f1-f2)
self.assertTrue(diff < 0.1 or diff/norm(f1) < 1e-4)
def testSwitchFunction(self):
""" Tests the switching function option in AmberPrmtopFile """
system = prmtop1.createSystem(nonbondedMethod=PME,
nonbondedCutoff=1*nanometer,
switchDistance=0.8*nanometer)
for force in system.getForces():
if isinstance(force, NonbondedForce):
self.assertTrue(force.getUseSwitchingFunction())
self.assertEqual(force.getSwitchingDistance(), 0.8*nanometer)
break
else:
assert False, 'Did not find expected nonbonded force!'
# Check error handling
system = prmtop1.createSystem(nonbondedMethod=PME,
nonbondedCutoff=1*nanometer)
for force in system.getForces():
if isinstance(force, NonbondedForce):
self.assertFalse(force.getUseSwitchingFunction())
break
else:
assert False, 'Did not find expected nonbonded force!'
self.assertRaises(ValueError, lambda:
prmtop1.createSystem(nonbondedMethod=PME,
nonbondedCutoff=1*nanometer, switchDistance=-1)
)
self.assertRaises(ValueError, lambda:
prmtop1.createSystem(nonbondedMethod=PME,
nonbondedCutoff=1*nanometer, switchDistance=1.2)
)
def test_with_dcd_reporter(self):
"""Check that an amber simulation like the docs example works with a DCD reporter."""
......
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