Skip to content

GitLab

Commit 3e4f3866 authored by Peter Eastman's avatar Peter Eastman
Browse files

Added missing support for LocalCoordinatesSite

parent 58f5168e
Hide whitespace changes
Inline Side-by-side
Showing with 102 additions and 30 deletions
wrappers/python/simtk/openmm/app/modeller.py
View file @ 3e4f3866
...@@ -1007,6 +1007,15 @@ class Modeller(object): ...@@ -1007,6 +1007,15 @@ class Modeller(object):
v2 = templateAtomPositions[site.atoms[2]] - templateAtomPositions[site.atoms[0]] v2 = templateAtomPositions[site.atoms[2]] - templateAtomPositions[site.atoms[0]]
cross = Vec3(v1[1]*v2[2]-v1[2]*v2[1], v1[2]*v2[0]-v1[0]*v2[2], v1[0]*v2[1]-v1[1]*v2[0]) cross = Vec3(v1[1]*v2[2]-v1[2]*v2[1], v1[2]*v2[0]-v1[0]*v2[2], v1[0]*v2[1]-v1[1]*v2[0])
position = templateAtomPositions[site.atoms[0]] + site.weights[0]*v1 + site.weights[1]*v2 + site.weights[2]*cross position = templateAtomPositions[site.atoms[0]] + site.weights[0]*v1 + site.weights[1]*v2 + site.weights[2]*cross
elif site.type == 'localCoords':
origin = templateAtomPositions[site.atoms[0]]*site.originWeights[0] + templateAtomPositions[site.atoms[1]]*site.originWeights[1] + templateAtomPositions[site.atoms[2]]*site.originWeights[2];
xdir = templateAtomPositions[site.atoms[0]]*site.xWeights[0] + templateAtomPositions[site.atoms[1]]*site.xWeights[1] + templateAtomPositions[site.atoms[2]]*site.xWeights[2];
ydir = templateAtomPositions[site.atoms[0]]*site.yWeights[0] + templateAtomPositions[site.atoms[1]]*site.yWeights[1] + templateAtomPositions[site.atoms[2]]*site.yWeights[2];
zdir = Vec3(xdir[1]*ydir[2]-xdir[2]*ydir[1], xdir[2]*ydir[0]-xdir[0]*ydir[2], xdir[0]*ydir[1]-xdir[1]*ydir[0])
xdir /= norm(xdir);
zdir /= norm(zdir);
ydir = Vec3(zdir[1]*xdir[2]-zdir[2]*xdir[1], zdir[2]*xdir[0]-zdir[0]*xdir[2], zdir[0]*xdir[1]-zdir[1]*xdir[0])
position = origin + xdir*site.localPos[0] + ydir*site.localPos[1] + zdir*site.localPos[2];
if position is None and atom.type in drudeTypeMap: if position is None and atom.type in drudeTypeMap:
# This is a Drude particle. Put it on top of its parent atom. # This is a Drude particle. Put it on top of its parent atom.
......
wrappers/python/src/swig_doxygen/swigInputBuilder.py
View file @ 3e4f3866
...@@ -279,7 +279,7 @@ class SwigInputBuilder: ...@@ -279,7 +279,7 @@ class SwigInputBuilder:
self.fOut.write(",\n OpenMM::%s" % name) self.fOut.write(",\n OpenMM::%s" % name)
self.fOut.write(");\n\n") self.fOut.write(");\n\n")
self.fOut.write("%factory(OpenMM::VirtualSite& OpenMM::System::getVirtualSite, OpenMM::TwoParticleAverageSite, OpenMM::ThreeParticleAverageSite, OpenMM::OutOfPlaneSite);\n\n") self.fOut.write("%factory(OpenMM::VirtualSite& OpenMM::System::getVirtualSite, OpenMM::TwoParticleAverageSite, OpenMM::ThreeParticleAverageSite, OpenMM::OutOfPlaneSite, OpenMM::LocalCoordinatesSite);\n\n")
self.fOut.write("\n") self.fOut.write("\n")
def writeGlobalConstants(self): def writeGlobalConstants(self):
......
wrappers/python/tests/TestModeller.py
View file @ 3e4f3866
...@@ -923,38 +923,101 @@ class TestModeller(unittest.TestCase): ...@@ -923,38 +923,101 @@ class TestModeller(unittest.TestCase):
self.assertEqual(1, len(ep)) self.assertEqual(1, len(ep))
def test_addVirtualSites(self):
"""Test adding extra particles defined by virtual sites."""
xml = """
<ForceField>
<AtomTypes>
<Type name="C" class="C" element="C" mass="10"/>
<Type name="N" class="N" element="N" mass="10"/>
<Type name="O" class="O" element="O" mass="10"/>
<Type name="V" class="V" mass="0.0"/>
</AtomTypes>
<Residues>
<Residue name="Test">
<Atom name="C" type="C"/>
<Atom name="N" type="N"/>
<Atom name="O" type="O"/>
<Atom name="V1" type="V"/>
<Atom name="V2" type="V"/>
<Atom name="V3" type="V"/>
<Atom name="V4" type="V"/>
<VirtualSite type="average2" index="3" atom1="0" atom2="1" weight1="0.7" weight2="0.3"/>
<VirtualSite type="average3" index="4" atom1="0" atom2="1" atom3="2" weight1="0.2" weight2="0.3" weight3="0.5"/>
<VirtualSite type="outOfPlane" index="5" atom1="0" atom2="1" atom3="2" weight12="0.1" weight13="-0.2" weightCross="0.8"/>
<VirtualSite type="localCoords" index="6" atom1="0" atom2="1" atom3="2" wo1="0.1" wo2="0.5" wo3="0.4" wx1="1" wx2="-0.6" wx3="-0.4" wy1="0.1" wy2="0.9" wy3="-1" p1="-0.5" p2="0.4" p3="1.1"/>
</Residue>
</Residues>
</ForceField>"""
ff = ForceField(StringIO(xml))
# Create the three real atoms.
topology = Topology()
chain = topology.addChain()
residue = topology.addResidue('Test', chain)
topology.addAtom('C', element.carbon, residue)
topology.addAtom('N', element.nitrogen, residue)
topology.addAtom('V', element.oxygen, residue)
# Add the virtual sites.
modeller = Modeller(topology, [Vec3(0.1, 0.2, 0.3), Vec3(1.0, 0.9, 0.8), Vec3(1.5, 1.1, 0.7)]*nanometers)
modeller.addExtraParticles(ff)
top = modeller.topology
pos = modeller.positions
# Check that the correct particles were added.
self.assertEqual(len(pos), 7)
for atom, elem in zip(top.atoms(), [element.carbon, element.nitrogen, element.oxygen, None, None, None, None]):
self.assertEqual(elem, atom.element)
# Check that the positions were calculated correctly.
system = ff.createSystem(top)
integ = VerletIntegrator(1.0)
context = Context(system, integ)
context.setPositions(pos)
context.computeVirtualSites()
state = context.getState(getPositions=True)
for p1, p2 in zip (pos, state.getPositions()):
self.assertVecAlmostEqual(p1.value_in_unit(nanometers), p2.value_in_unit(nanometers), 1e-6)
def test_multiSiteIon(self): def test_multiSiteIon(self):
"""Test adding extra particles whose positions are determined based on bonds.""" """Test adding extra particles whose positions are determined based on bonds."""
xml = """ xml = """
<ForceField> <ForceField>
<AtomTypes> <AtomTypes>
<Type name="Zn" class="Zn" element="Zn" mass="53.380"/> <Type name="Zn" class="Zn" element="Zn" mass="53.380"/>
<Type name="DA" class="DA" mass="3.0"/> <Type name="DA" class="DA" mass="3.0"/>
</AtomTypes> </AtomTypes>
<Residues> <Residues>
<Residue name="ZN"> <Residue name="ZN">
<Atom name="ZN" type="Zn"/> <Atom name="ZN" type="Zn"/>
<Atom name="D1" type="DA"/> <Atom name="D1" type="DA"/>
<Atom name="D2" type="DA"/> <Atom name="D2" type="DA"/>
<Atom name="D3" type="DA"/> <Atom name="D3" type="DA"/>
<Atom name="D4" type="DA"/> <Atom name="D4" type="DA"/>
<Bond from="0" to="2"/> <Bond from="0" to="2"/>
<Bond from="0" to="1"/> <Bond from="0" to="1"/>
<Bond from="0" to="3"/> <Bond from="0" to="3"/>
<Bond from="0" to="4"/> <Bond from="0" to="4"/>
<Bond from="1" to="2"/> <Bond from="1" to="2"/>
<Bond from="1" to="3"/> <Bond from="1" to="3"/>
<Bond from="1" to="4"/> <Bond from="1" to="4"/>
<Bond from="2" to="4"/> <Bond from="2" to="4"/>
<Bond from="2" to="3"/> <Bond from="2" to="3"/>
<Bond from="3" to="4"/> <Bond from="3" to="4"/>
</Residue> </Residue>
</Residues> </Residues>
<HarmonicBondForce> <HarmonicBondForce>
<Bond class1="DA" class2="Zn" length="0.09" k="535552.0"/> <Bond class1="DA" class2="Zn" length="0.09" k="535552.0"/>
<Bond class1="DA" class2="DA" length="0.147" k="535552.0"/> <Bond class1="DA" class2="DA" length="0.147" k="535552.0"/>
</HarmonicBondForce> </HarmonicBondForce>
</ForceField>""" </ForceField>"""
ff = ForceField(StringIO(xml)) ff = ForceField(StringIO(xml))
# Create two zinc atoms. # Create two zinc atoms.
......
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