Skip to content

GitLab

Commit 4570189f authored by Jason Swails's avatar Jason Swails
Browse files

Merge branch 'master' into fix-charmm-box-vecs

parents d127d7c2 26790496
Hide whitespace changes
Inline Side-by-side
Showing with 364 additions and 23 deletions
wrappers/python/simtk/openmm/app/topology.py
View file @ 4570189f
......@@ -6,7 +6,7 @@ Simbios, the NIH National Center for Physics-Based Simulation of
Biological Structures at Stanford, funded under the NIH Roadmap for
Medical Research, grant U54 GM072970. See https://simtk.org.
Portions copyright (c) 2012-2013 Stanford University and the Authors.
Portions copyright (c) 2012-2015 Stanford University and the Authors.
Authors: Peter Eastman
Contributors:
......@@ -33,7 +33,9 @@ __version__ = "1.0"
import os
import xml.etree.ElementTree as etree
from simtk.unit import nanometers, sqrt
from simtk.openmm.vec3 import Vec3
from simtk.unit import nanometers, sqrt, is_quantity
from copy import deepcopy
class Topology(object):
"""Topology stores the topological information about a system.
......@@ -55,7 +57,7 @@ class Topology(object):
self._numResidues = 0
self._numAtoms = 0
self._bonds = []
self._unitCellDimensions = None
self._periodicBoxVectors = None
def addChain(self):
"""Create a new Chain and add it to the Topology.
......@@ -123,16 +125,48 @@ class Topology(object):
"""Iterate over all bonds (each represented as a tuple of two Atoms) in the Topology."""
return iter(self._bonds)
def getPeriodicBoxVectors(self):
"""Get the vectors defining the periodic box.
The return value may be None if this Topology does not represent a periodic structure."""
return self._periodicBoxVectors
def setPeriodicBoxVectors(self, vectors):
"""Set the vectors defining the periodic box."""
if vectors is not None:
if not is_quantity(vectors[0][0]):
vectors = vectors*nanometers
if vectors[0][1] != 0*nanometers or vectors[0][2] != 0*nanometers:
raise ValueError("First periodic box vector must be parallel to x.");
if vectors[1][2] != 0*nanometers:
raise ValueError("Second periodic box vector must be in the x-y plane.");
if vectors[0][0] <= 0*nanometers or vectors[1][1] <= 0*nanometers or vectors[2][2] <= 0*nanometers or vectors[0][0] < 2*abs(vectors[1][0]) or vectors[0][0] < 2*abs(vectors[2][0]) or vectors[1][1] < 2*abs(vectors[2][1]):
raise ValueError("Periodic box vectors must be in reduced form.");
self._periodicBoxVectors = deepcopy(vectors)
def getUnitCellDimensions(self):
"""Get the dimensions of the crystallographic unit cell.
The return value may be None if this Topology does not represent a periodic structure.
"""
return self._unitCellDimensions
if self._periodicBoxVectors is None:
return None
xsize = self._periodicBoxVectors[0][0].value_in_unit(nanometers)
ysize = self._periodicBoxVectors[1][1].value_in_unit(nanometers)
zsize = self._periodicBoxVectors[2][2].value_in_unit(nanometers)
return Vec3(xsize, ysize, zsize)*nanometers
def setUnitCellDimensions(self, dimensions):
"""Set the dimensions of the crystallographic unit cell."""
self._unitCellDimensions = dimensions
"""Set the dimensions of the crystallographic unit cell.
This method is an alternative to setPeriodicBoxVectors() for the case of a rectangular box. It sets
the box vectors to be orthogonal to each other and to have the specified lengths."""
if dimensions is None:
self._periodicBoxVectors = None
else:
if is_quantity(dimensions):
dimensions = dimensions.value_in_unit(nanometers)
self._periodicBoxVectors = (Vec3(dimensions[0], 0, 0), Vec3(0, dimensions[1], 0), Vec3(0, 0, dimensions[2]))*nanometers
@staticmethod
def loadBondDefinitions(file):
......
wrappers/python/tests/TestForceField.py
View file @ 4570189f
......@@ -20,7 +20,7 @@ class TestForceField(unittest.TestCase):
self.topology1 = self.pdb1.topology
self.topology1.setUnitCellDimensions(Vec3(2, 2, 2))
# alalnine dipeptide with implicit water
# alanine dipeptide with implicit water
self.pdb2 = PDBFile('systems/alanine-dipeptide-implicit.pdb')
self.forcefield2 = ForceField('amber99sb.xml', 'amber99_obc.xml')
......@@ -186,6 +186,18 @@ class TestForceField(unittest.TestCase):
system2 = ff2.createSystem(modeller.topology)
self.assertEqual(XmlSerializer.serialize(system1), XmlSerializer.serialize(system2))
def test_PeriodicBoxVectors(self):
"""Test setting the periodic box vectors."""
vectors = (Vec3(5, 0, 0), Vec3(-1.5, 4.5, 0), Vec3(0.4, 0.8, 7.5))*nanometers
self.pdb1.topology.setPeriodicBoxVectors(vectors)
self.assertEqual(Vec3(5, 4.5, 7.5)*nanometers, self.pdb1.topology.getUnitCellDimensions())
system = self.forcefield1.createSystem(self.pdb1.topology)
for i in range(3):
self.assertEqual(vectors[i], self.pdb1.topology.getPeriodicBoxVectors()[i])
self.assertEqual(vectors[i], system.getDefaultPeriodicBoxVectors()[i])
class AmoebaTestForceField(unittest.TestCase):
"""Test the ForceField.createSystem() method with the AMOEBA forcefield."""
......
wrappers/python/tests/TestGromacsGroFile.py 0 → 100644
View file @ 4570189f
import unittest
from simtk.openmm.app import *
from simtk.openmm import *
from simtk.unit import *
import simtk.openmm.app.element as elem
class TestGromacsGroFile(unittest.TestCase):
"""Test the Gromacs GRO file parser"""
def test_Triclinic(self):
"""Test parsing a file that describes a triclinic box."""
gro = GromacsGroFile('systems/triclinic.gro')
self.assertEqual(len(gro.positions), 8)
expectedPositions = [
Vec3(1.744, 2.788, 3.162),
Vec3(1.048, 0.762, 2.340),
Vec3(2.489, 1.570, 2.817),
Vec3(1.027, 1.893, 3.271),
Vec3(0.937, 0.825, 0.009),
Vec3(2.290, 1.887, 3.352),
Vec3(1.266, 1.111, 2.894),
Vec3(0.933, 1.862, 3.490)]*nanometers
for (p1, p2) in zip(expectedPositions, gro.positions):
self.assertEqual(p1, p2)
expectedVectors = [
Vec3(2.5, 0, 0),
Vec3(0.5, 3.0, 0),
Vec3(0.7, 0.9, 3.5)]*nanometers
for (v1, v2) in zip(expectedVectors, gro.getPeriodicBoxVectors()):
self.assertEqual(v1, v2)
self.assertEqual(Vec3(2.5, 3.0, 3.5)*nanometers, gro.getUnitCellDimensions())
for i in range(4):
self.assertEqual(elem.chlorine, gro.elements[i])
self.assertEqual('Cl', gro.atomNames[i])
self.assertEqual('Cl', gro.residueNames[i])
for i in range(4, 8):
self.assertEqual(elem.sodium, gro.elements[i])
self.assertEqual('Na', gro.atomNames[i])
self.assertEqual('Na', gro.residueNames[i])
if __name__ == '__main__':
unittest.main()
wrappers/python/tests/TestModeller.py
View file @ 4570189f
......@@ -308,44 +308,56 @@ class TestModeller(unittest.TestCase):
self.assertTrue(len(matoms)==0 and len(m1atoms)==1 and len(m2atoms)==1)
def test_addSolventPeriodicBox(self):
""" Test the addSolvent() method; test that the three ways of passing in the periodic box all work. """
""" Test the addSolvent() method; test that the four ways of passing in the periodic box all work. """
# First way of passing in periodic box dimensions: set it in the original topology.
# First way of passing in periodic box vectors: set it in the original topology.
topology_start = self.pdb.topology
topology_start.setUnitCellDimensions(Vec3(3.5, 4.5, 5.5)*nanometers)
modeller = Modeller(topology_start, self.positions)
modeller.deleteWater()
modeller.addSolvent(self.forcefield)
topology_after = modeller.getTopology()
dim3 = topology_after.getUnitCellDimensions()
dim3 = topology_after.getPeriodicBoxVectors()
self.assertAlmostEqual(dim3[0]/nanometers, 3.5)
self.assertAlmostEqual(dim3[1]/nanometers, 4.5)
self.assertAlmostEqual(dim3[2]/nanometers, 5.5)
self.assertVecAlmostEqual(dim3[0]/nanometers, Vec3(3.5, 0, 0))
self.assertVecAlmostEqual(dim3[1]/nanometers, Vec3(0, 4.5, 0))
self.assertVecAlmostEqual(dim3[2]/nanometers, Vec3(0, 0, 5.5))
# Second way of passing in the periodic box dimensions: as a parameter to addSolvent()
# Second way of passing in the periodic box vectors: with the boxSize parameter to addSolvent()
topology_start = self.pdb.topology
modeller = Modeller(topology_start, self.positions)
modeller.deleteWater()
modeller.addSolvent(self.forcefield, boxSize = Vec3(3.6, 4.6, 5.6)*nanometers)
topology_after = modeller.getTopology()
dim3 = topology_after.getUnitCellDimensions()
dim3 = topology_after.getPeriodicBoxVectors()
self.assertAlmostEqual(dim3[0]/nanometers, 3.6)
self.assertAlmostEqual(dim3[1]/nanometers, 4.6)
self.assertAlmostEqual(dim3[2]/nanometers, 5.6)
self.assertVecAlmostEqual(dim3[0]/nanometers, Vec3(3.6, 0, 0))
self.assertVecAlmostEqual(dim3[1]/nanometers, Vec3(0, 4.6, 0))
self.assertVecAlmostEqual(dim3[2]/nanometers, Vec3(0, 0, 5.6))
# Third way of passing in the periodic box dimensions: pass a 'padding' value to addSolvent()
# Third way of passing in the periodic box vectors: with the boxVectors parameter to addSolvent()
topology_start = self.pdb.topology
modeller = Modeller(topology_start, self.positions)
modeller.deleteWater()
modeller.addSolvent(self.forcefield, boxVectors = (Vec3(3.4, 0, 0), Vec3(0.5, 4.4, 0), Vec3(-1.0, -1.5, 5.4))*nanometers)
topology_after = modeller.getTopology()
dim3 = topology_after.getPeriodicBoxVectors()
self.assertVecAlmostEqual(dim3[0]/nanometers, Vec3(3.4, 0, 0))
self.assertVecAlmostEqual(dim3[1]/nanometers, Vec3(0.5, 4.4, 0))
self.assertVecAlmostEqual(dim3[2]/nanometers, Vec3(-1.0, -1.5, 5.4))
# Fourth way of passing in the periodic box vectors: pass a 'padding' value to addSolvent()
topology_start = self.pdb.topology
modeller = Modeller(topology_start, self.positions)
modeller.deleteWater()
modeller.addSolvent(self.forcefield, padding = 1.0*nanometers)
topology_after = modeller.getTopology()
dim3 = topology_after.getUnitCellDimensions()
dim3 = topology_after.getPeriodicBoxVectors()
self.assertAlmostEqual(dim3[0]/nanometers, 2.8802)
self.assertAlmostEqual(dim3[1]/nanometers, 2.8802)
self.assertAlmostEqual(dim3[2]/nanometers, 2.8802)
self.assertVecAlmostEqual(dim3[0]/nanometers, Vec3(2.8802, 0, 0))
self.assertVecAlmostEqual(dim3[1]/nanometers, Vec3(0, 2.8802, 0))
self.assertVecAlmostEqual(dim3[2]/nanometers, Vec3(0, 0, 2.8802))
def test_addSolventNeutralSolvent(self):
""" Test the addSolvent() method; test adding ions to neutral solvent. """
......@@ -874,6 +886,12 @@ class TestModeller(unittest.TestCase):
validate_equivalence(self, topology_LYN, topology_after)
def assertVecAlmostEqual(self, p1, p2, tol=1e-7):
scale = max(1.0, norm(p1),)
for i in range(3):
diff = abs(p1[i]-p2[i])/scale
self.assertTrue(diff < tol)
if __name__ == '__main__':
unittest.main()
......
wrappers/python/tests/TestPdbFile.py 0 → 100644
View file @ 4570189f
import unittest
from simtk.openmm.app import *
from simtk.openmm import *
from simtk.unit import *
import simtk.openmm.app.element as elem
import cStringIO
class TestPdbFile(unittest.TestCase):
"""Test the PDB file parser"""
def test_Triclinic(self):
"""Test parsing a file that describes a triclinic box."""
pdb = PDBFile('systems/triclinic.pdb')
self.assertEqual(len(pdb.positions), 8)
expectedPositions = [
Vec3(1.744, 2.788, 3.162),
Vec3(1.048, 0.762, 2.340),
Vec3(2.489, 1.570, 2.817),
Vec3(1.027, 1.893, 3.271),
Vec3(0.937, 0.825, 0.009),
Vec3(2.290, 1.887, 3.352),
Vec3(1.266, 1.111, 2.894),
Vec3(0.933, 1.862, 3.490)]*nanometers
for (p1, p2) in zip(expectedPositions, pdb.positions):
self.assertVecAlmostEqual(p1, p2)
expectedVectors = [
Vec3(2.5, 0, 0),
Vec3(0.5, 3.0, 0),
Vec3(0.7, 0.9, 3.5)]*nanometers
for (v1, v2) in zip(expectedVectors, pdb.topology.getPeriodicBoxVectors()):
self.assertVecAlmostEqual(v1, v2, 1e-4)
self.assertVecAlmostEqual(Vec3(2.5, 3.0, 3.5)*nanometers, pdb.topology.getUnitCellDimensions(), 1e-4)
for atom in pdb.topology.atoms():
if atom.index < 4:
self.assertEqual(elem.chlorine, atom.element)
self.assertEqual('Cl', atom.name)
self.assertEqual('Cl', atom.residue.name)
else:
self.assertEqual(elem.sodium, atom.element)
self.assertEqual('Na', atom.name)
self.assertEqual('Na', atom.residue.name)
def test_WriteFile(self):
"""Write a file, read it back, and make sure it matches the original."""
pdb1 = PDBFile('systems/triclinic.pdb')
output = cStringIO.StringIO()
PDBFile.writeFile(pdb1.topology, pdb1.positions, output)
input = cStringIO.StringIO(output.getvalue())
pdb2 = PDBFile(input)
output.close();
input.close();
self.assertEqual(len(pdb2.positions), 8)
for (p1, p2) in zip(pdb1.positions, pdb2.positions):
self.assertVecAlmostEqual(p1, p2)
for (v1, v2) in zip(pdb1.topology.getPeriodicBoxVectors(), pdb2.topology.getPeriodicBoxVectors()):
self.assertVecAlmostEqual(v1, v2, 1e-4)
self.assertVecAlmostEqual(pdb1.topology.getUnitCellDimensions(), pdb2.topology.getUnitCellDimensions(), 1e-4)
for atom1, atom2 in zip(pdb1.topology.atoms(), pdb2.topology.atoms()):
self.assertEqual(atom1.element, atom2.element)
self.assertEqual(atom1.name, atom2.name)
self.assertEqual(atom1.residue.name, atom2.residue.name)
def assertVecAlmostEqual(self, p1, p2, tol=1e-7):
unit = p1.unit
p1 = p1.value_in_unit(unit)
p2 = p2.value_in_unit(unit)
scale = max(1.0, norm(p1),)
for i in range(3):
diff = abs(p1[i]-p2[i])/scale
self.assertTrue(diff < tol)
if __name__ == '__main__':
unittest.main()
wrappers/python/tests/TestPdbxFile.py 0 → 100644
View file @ 4570189f
import unittest
from simtk.openmm.app import *
from simtk.openmm import *
from simtk.unit import *
import simtk.openmm.app.element as elem
class TestPdbxFile(unittest.TestCase):
"""Test the PDBx/mmCIF file parser"""
def test_Triclinic(self):
"""Test parsing a file that describes a triclinic box."""
pdb = PDBxFile('systems/triclinic.pdbx')
self.assertEqual(len(pdb.positions), 8)
expectedPositions = [
Vec3(1.744, 2.788, 3.162),
Vec3(1.048, 0.762, 2.340),
Vec3(2.489, 1.570, 2.817),
Vec3(1.027, 1.893, 3.271),
Vec3(0.937, 0.825, 0.009),
Vec3(2.290, 1.887, 3.352),
Vec3(1.266, 1.111, 2.894),
Vec3(0.933, 1.862, 3.490)]*nanometers
for (p1, p2) in zip(expectedPositions, pdb.positions):
self.assertVecAlmostEqual(p1, p2)
expectedVectors = [
Vec3(2.5, 0, 0),
Vec3(0.5, 3.0, 0),
Vec3(0.7, 0.9, 3.5)]*nanometers
for (v1, v2) in zip(expectedVectors, pdb.topology.getPeriodicBoxVectors()):
self.assertVecAlmostEqual(v1, v2, 1e-6)
self.assertVecAlmostEqual(Vec3(2.5, 3.0, 3.5)*nanometers, pdb.topology.getUnitCellDimensions(), 1e-6)
for atom in pdb.topology.atoms():
if atom.index < 4:
self.assertEqual(elem.chlorine, atom.element)
self.assertEqual('Cl', atom.name)
self.assertEqual('Cl', atom.residue.name)
else:
self.assertEqual(elem.sodium, atom.element)
self.assertEqual('Na', atom.name)
self.assertEqual('Na', atom.residue.name)
def assertVecAlmostEqual(self, p1, p2, tol=1e-7):
unit = p1.unit
p1 = p1.value_in_unit(unit)
p2 = p2.value_in_unit(unit)
scale = max(1.0, norm(p1),)
for i in range(3):
diff = abs(p1[i]-p2[i])/scale
self.assertTrue(diff < tol)
if __name__ == '__main__':
unittest.main()
wrappers/python/tests/systems/triclinic.gro 0 → 100644
View file @ 4570189f
Triclinic test system
8
1Cl Cl 1 1.744 2.788 3.162
2Cl Cl 2 1.048 0.762 2.340
3Cl Cl 3 2.489 1.570 2.817
4Cl Cl 4 1.027 1.893 3.271
5Na Na 5 0.937 0.825 0.009
6Na Na 6 2.290 1.887 3.352
7Na Na 7 1.266 1.111 2.894
8Na Na 8 0.933 1.862 3.490
2.50000 3.00000 3.50000 0.00000 0.00000 0.50000 0.00000 0.70000 0.90000
wrappers/python/tests/systems/triclinic.pdb 0 → 100644
View file @ 4570189f
CRYST1 25.000 30.414 36.810 74.19 79.04 80.54 P 1 1
ATOM 1 Cl Cl A 1 17.440 27.880 31.620 1.00 0.00 Cl
ATOM 2 Cl Cl A 2 10.480 7.620 23.400 1.00 0.00 Cl
ATOM 3 Cl Cl A 3 24.890 15.700 28.170 1.00 0.00 Cl
ATOM 4 Cl Cl A 4 10.270 18.930 32.710 1.00 0.00 Cl
ATOM 5 Na Na A 5 9.370 8.250 0.090 1.00 0.00 Na
ATOM 6 Na Na A 6 22.900 18.870 33.520 1.00 0.00 Na
ATOM 7 Na Na A 7 12.660 11.110 28.940 1.00 0.00 Na
ATOM 8 Na Na A 8 9.330 18.620 34.900 1.00 0.00 Na
TER 9 Na A 8
END
wrappers/python/tests/systems/triclinic.pdbx 0 → 100644
View file @ 4570189f
data_PHPWKVVWM
#
loop_
_pdbx_nonpoly_scheme.asym_id
_pdbx_nonpoly_scheme.entity_id
_pdbx_nonpoly_scheme.mon_id
_pdbx_nonpoly_scheme.ndb_seq_num
_pdbx_nonpoly_scheme.pdb_seq_num
_pdbx_nonpoly_scheme.auth_seq_num
_pdbx_nonpoly_scheme.pdb_mon_id
_pdbx_nonpoly_scheme.auth_mon_id
_pdbx_nonpoly_scheme.pdb_strand_id
_pdbx_nonpoly_scheme.pdb_ins_code
A 1 Cl 1 1 1 Cl Cl A .
B 1 Cl 1 2 2 Cl Cl A .
C 1 Cl 1 3 3 Cl Cl A .
D 1 Cl 1 4 4 Cl Cl A .
E 2 Na 1 5 5 Na Na A .
F 2 Na 1 6 6 Na Na A .
G 2 Na 1 7 7 Na Na A .
H 2 Na 1 8 8 Na Na A .
#
_cell.entry_id PHPWKVVWM
_cell.length_a 25.000
_cell.length_b 30.4138126515
_cell.length_c 36.8103246386
_cell.angle_alpha 74.1909185129
_cell.angle_beta 79.0376425369
_cell.angle_gamma 80.537677792
_cell.Z_PDB 1
_cell.pdbx_unique_axis ?
#
loop_
_atom_type.symbol
Cl
Na
#
loop_
_pdbx_entity_nonpoly.entity_id
_pdbx_entity_nonpoly.name
_pdbx_entity_nonpoly.comp_id
1 ? Cl
2 ? Na
#
loop_
_atom_site.group_PDB
_atom_site.id
_atom_site.type_symbol
_atom_site.label_atom_id
_atom_site.label_alt_id
_atom_site.label_comp_id
_atom_site.label_asym_id
_atom_site.label_entity_id
_atom_site.label_seq_id
_atom_site.pdbx_label_seq_num
_atom_site.pdbx_PDB_ins_code
_atom_site.Cartn_x
_atom_site.Cartn_y
_atom_site.Cartn_z
_atom_site.occupancy
_atom_site.B_iso_or_equiv
_atom_site.Cartn_x_esd
_atom_site.Cartn_y_esd
_atom_site.Cartn_z_esd
_atom_site.occupancy_esd
_atom_site.B_iso_or_equiv_esd
_atom_site.footnote_id
_atom_site.pdbx_formal_charge
_atom_site.auth_seq_id
_atom_site.auth_comp_id
_atom_site.auth_asym_id
_atom_site.auth_atom_id
_atom_site.pdbx_auth_seq_id
_atom_site.pdbx_auth_comp_id
_atom_site.pdbx_auth_asym_id
_atom_site.pdbx_auth_atom_name
_atom_site.pdbx_PDB_model_num
HETATM 1 Cl Cl . Cl A 1 . 1 ? 17.440 27.880 31.620 1.00 0.00 ? ? ? ? ? ? ? 1 Cl A Cl 1 Cl A Cl 1
HETATM 2 Cl Cl . Cl B 1 . 1 ? 10.480 7.620 23.400 1.00 0.00 ? ? ? ? ? ? ? 2 Cl A Cl 2 Cl A Cl 1
HETATM 3 Cl Cl . Cl C 1 . 1 ? 24.890 15.700 28.170 1.00 0.00 ? ? ? ? ? ? ? 3 Cl A Cl 3 Cl A Cl 1
HETATM 4 Cl Cl . Cl D 1 . 1 ? 10.270 18.930 32.710 1.00 0.00 ? ? ? ? ? ? ? 4 Cl A Cl 4 Cl A Cl 1
HETATM 5 Na Na . Na E 2 . 1 ? 9.370 8.250 0.090 1.00 0.00 ? ? ? ? ? ? ? 5 Na A Na 5 Na A Na 1
HETATM 6 Na Na . Na F 2 . 1 ? 22.900 18.870 33.520 1.00 0.00 ? ? ? ? ? ? ? 6 Na A Na 6 Na A Na 1
HETATM 7 Na Na . Na G 2 . 1 ? 12.660 11.110 28.940 1.00 0.00 ? ? ? ? ? ? ? 7 Na A Na 7 Na A Na 1
HETATM 8 Na Na . Na H 2 . 1 ? 9.330 18.620 34.900 1.00 0.00 ? ? ? ? ? ? ? 8 Na A Na 8 Na A Na 1
#
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