Continuing triclinic box support in Python layer

wrappers/python/CMakeLists.txt

@@ -54,6 +54,7 @@ foreach(SUBDIR ${SUBDIRS})
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.sh"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.xml"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.pdb"
+        "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.pdbx"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.prmtop"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.prm"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.inpcrd"

wrappers/python/simtk/openmm/app/gromacsgrofile.py

@@ -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 Stanford University and the Authors.
+Portions copyright (c) 2012-2015 Stanford University and the Authors.
 Authors: Lee-Ping Wang, Peter Eastman
 Contributors:
 

wrappers/python/simtk/openmm/app/gromacstopfile.py

@@ -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-2014 Stanford University and the Authors.
+Portions copyright (c) 2012-2015 Stanford University and the Authors.
 Authors: Peter Eastman
 Contributors:
 
@@ -358,12 +358,14 @@ class GromacsTopFile(object):
             raise ValueError('Unsupported function type in [ cmaptypes ] line: '+line);
         self._cmapTypes[tuple(fields[:5])] = fields
 
-    def __init__(self, file, unitCellDimensions=None, includeDir=None, defines=None):
+    def __init__(self, file, periodicBoxVectors=None, unitCellDimensions=None, includeDir=None, defines=None):
         """Load a top file.
 
         Parameters:
          - file (string) the name of the file to load
-         - unitCellDimensions (Vec3=None) the dimensions of the crystallographic unit cell
+         - periodicBoxVectors (tuple of Vec3=None) the vectors defining the periodic box
+         - unitCellDimensions (Vec3=None) the dimensions of the crystallographic unit cell.  For
+           non-rectangular unit cells, specify periodicBoxVectors instead.
          - includeDir (string=None) A directory in which to look for other files
            included from the top file. If not specified, we will attempt to locate a gromacs
            installation on your system. When gromacs is installed in /usr/local, this will resolve
@@ -403,6 +405,11 @@ class GromacsTopFile(object):
         top = Topology()
         ## The Topology read from the prmtop file
         self.topology = top
+        if periodicBoxVectors is not None:
+            if unitCellDimensions is not None:
+                raise ValueError("specify either periodicBoxVectors or unitCellDimensions, but not both")
+            top.setPeriodicBoxVectors(periodicBoxVectors)
+        else:
             top.setUnitCellDimensions(unitCellDimensions)
         PDBFile._loadNameReplacementTables()
         for moleculeName, moleculeCount in self._molecules:

wrappers/python/simtk/openmm/app/internal/pdbstructure.py

@@ -8,7 +8,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: Christopher M. Bruns
 Contributors: Peter Eastman
 
@@ -41,6 +41,38 @@ import warnings
 import sys
 import math
 
+def computePeriodicBoxVectors(a_length, b_length, c_length, alpha, beta, gamma):
+    """Convert lengths and angles from a CRYST1 record to periodic box vectors."""
+
+    # Compute the vectors.
+
+    a = [a_length, 0, 0]
+    b = [b_length*math.cos(gamma), b_length*math.sin(gamma), 0]
+    cx = c_length*math.cos(beta)
+    cy = c_length*(math.cos(alpha)-math.cos(beta)*math.cos(gamma))/math.sin(gamma)
+    cz = math.sqrt(c_length*c_length-cx*cx-cy*cy)
+    c = [cx, cy, cz]
+
+    # If any elements are very close to 0, set them to exactly 0.
+
+    for i in range(3):
+        if abs(a[i]) < 1e-6:
+            a[i] = 0.0
+        if abs(b[i]) < 1e-6:
+            b[i] = 0.0
+        if abs(c[i]) < 1e-6:
+            c[i] = 0.0
+    a = Vec3(*a)
+    b = Vec3(*b)
+    c = Vec3(*c)
+
+    # Make sure they're in the reduced form required by OpenMM.
+
+    c = c - b*round(c[1]/b[1])
+    c = c - a*round(c[0]/a[0])
+    b = b - a*round(b[0]/a[0])
+    return (a, b, c)*unit.angstroms
+
 class PdbStructure(object):
     """
     PdbStructure object holds a parsed Protein Data Bank format file.
@@ -171,7 +203,13 @@ class PdbStructure(object):
                 self._current_model._current_chain._add_ter_record()
                 self._reset_residue_numbers()
             elif (pdb_line.find("CRYST1") == 0):
-                self._compute_periodic_box_vectors(pdb_line)
+                a_length = float(pdb_line[6:15])
+                b_length = float(pdb_line[15:24])
+                c_length = float(pdb_line[24:33])
+                alpha = float(pdb_line[33:40])*math.pi/180.0
+                beta = float(pdb_line[40:47])*math.pi/180.0
+                gamma = float(pdb_line[47:54])*math.pi/180.0
+                self._periodic_box_vectors = computePeriodicBoxVectors(a_length, b_length, c_length, alpha, beta, gamma)
             elif (pdb_line.find("CONECT") == 0):
                 atoms = [int(pdb_line[6:11])]
                 for pos in (11,16,21,26):
@@ -189,29 +227,6 @@ class PdbStructure(object):
                 self.modified_residues.append(ModifiedResidue(pdb_line[16], int(pdb_line[18:22]), pdb_line[12:15].strip(), pdb_line[24:27].strip()))
         self._finalize()
 
-    def _compute_periodic_box_vectors(self, line):
-        """Parse a CRYST1 record to compute the periodic box vectors."""
-        a_length = float(line[6:15])
-        b_length = float(line[15:24])
-        c_length = float(line[24:33])
-        alpha = float(line[33:40])*math.pi/180.0
-        beta = float(line[40:47])*math.pi/180.0
-        gamma = float(line[47:54])*math.pi/180.0
-        a = [a_length, 0, 0]
-        b = [b_length*math.cos(gamma), b_length*math.sin(gamma), 0]
-        cx = c_length*math.cos(beta)
-        cy = c_length*(math.cos(alpha)-math.cos(beta)*math.cos(gamma))
-        cz = math.sqrt(c_length*c_length-cx*cx-cy*cy)
-        c = [cx, cy, cz]
-        for i in range(3):
-            if abs(a[i]) < 1e-6:
-                a[i] = 0.0
-            if abs(b[i]) < 1e-6:
-                b[i] = 0.0
-            if abs(c[i]) < 1e-6:
-                c[i] = 0.0
-        self._periodic_box_vectors = (Vec3(*a), Vec3(*b), Vec3(*c))*unit.angstroms
-
     def _reset_atom_numbers(self):
         self._atom_numbers_are_hex = False
         self._next_atom_number = 1

wrappers/python/simtk/openmm/app/pdbfile.py

@@ -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 Stanford University and the Authors.
+Portions copyright (c) 2012-2015 Stanford University and the Authors.
 Authors: Peter Eastman
 Contributors:
 
@@ -40,7 +40,7 @@ from datetime import date
 from simtk.openmm import Vec3, Platform
 from simtk.openmm.app.internal.pdbstructure import PdbStructure
 from simtk.openmm.app import Topology
-from simtk.unit import nanometers, angstroms, is_quantity, norm, Quantity
+from simtk.unit import nanometers, angstroms, is_quantity, norm, Quantity, dot
 import element as elem
 try:
     import numpy
@@ -250,10 +250,16 @@ class PDBFile(object):
          - file (file=stdout) A file to write the file to
         """
         print >>file, "REMARK   1 CREATED WITH OPENMM %s, %s" % (Platform.getOpenMMVersion(), str(date.today()))
-        boxSize = topology.getUnitCellDimensions()
-        if boxSize is not None:
-            size = boxSize.value_in_unit(angstroms)
-            print >>file, "CRYST1%9.3f%9.3f%9.3f  90.00  90.00  90.00 P 1           1 " % size
+        vectors = topology.getPeriodicBoxVectors()
+        if vectors is not None:
+            (a, b, c) = vectors.value_in_unit(angstroms)
+            a_length = norm(a)
+            b_length = norm(b)
+            c_length = norm(c)
+            alpha = math.acos(dot(b, c)/(b_length*c_length))*180.0/math.pi
+            beta = math.acos(dot(c, a)/(c_length*a_length))*180.0/math.pi
+            gamma = math.acos(dot(a, b)/(a_length*b_length))*180.0/math.pi
+            print >>file, "CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f P 1           1 " % (a_length, b_length, c_length, alpha, beta, gamma)
 
     @staticmethod
     def writeModel(topology, positions, file=sys.stdout, modelIndex=None):

wrappers/python/simtk/openmm/app/pdbxfile.py

@@ -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) 2014 Stanford University and the Authors.
+Portions copyright (c) 2014-2015 Stanford University and the Authors.
 Authors: Peter Eastman
 Contributors:
 
@@ -36,6 +36,7 @@ import sys
 import math
 from simtk.openmm import Vec3
 from simtk.openmm.app.internal.pdbx.reader.PdbxReader import PdbxReader
+from simtk.openmm.app.internal.pdbstructure import computePeriodicBoxVectors
 from simtk.openmm.app import Topology
 from simtk.unit import nanometers, angstroms, is_quantity, norm, Quantity
 import element as elem
@@ -147,21 +148,8 @@ class PDBxFile(object):
         if cell is not None and cell.getRowCount() > 0:
             row = cell.getRow(0)
             (a, b, c) = [float(row[cell.getAttributeIndex(attribute)]) for attribute in ('length_a', 'length_b', 'length_c')]
-            (alpha, beta, gamma) = [float(row[cell.getAttributeIndex(attribute)]) for attribute in ('angle_alpha', 'angle_beta', 'angle_gamma')*math.pi/180.0]
-            a = [a_length, 0, 0]
-            b = [b_length*math.cos(gamma), b_length*math.sin(gamma), 0]
-            cx = c_length*math.cos(beta)
-            cy = c_length*(math.cos(alpha)-math.cos(beta)*math.cos(gamma))
-            cz = math.sqrt(c_length*c_length-cx*cx-cy*cy)
-            c = [cx, cy, cz]
-            for i in range(3):
-                if abs(a[i]) < 1e-6:
-                    a[i] = 0.0
-                if abs(b[i]) < 1e-6:
-                    b[i] = 0.0
-                if abs(c[i]) < 1e-6:
-                    c[i] = 0.0
-            self.topology.setPeriodicBoxVectors((Vec3(*a), Vec3(*b), Vec3(*c))*unit.angstroms)
+            (alpha, beta, gamma) = [float(row[cell.getAttributeIndex(attribute)])*math.pi/180.0 for attribute in ('angle_alpha', 'angle_beta', 'angle_gamma')]
+            self.topology.setPeriodicBoxVectors(computePeriodicBoxVectors(a, b, c, alpha, beta, gamma))
 
         # Add bonds based on struct_conn records.
 

wrappers/python/tests/TestPdbFile.py

+import unittest
+from simtk.openmm.app import *
+from simtk.openmm import *
+from simtk.unit import *
+import simtk.openmm.app.element as elem
+
+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 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

+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.pdb

+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

+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 
+#