Adding triclinic box support to Python layer

wrappers/python/CMakeLists.txt

@@ -58,6 +58,7 @@ foreach(SUBDIR ${SUBDIRS})
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.prm"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.inpcrd"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.crd"
+        "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.gro"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.parm7"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.rst7"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.ncrst"

wrappers/python/simtk/openmm/app/forcefield.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, Mark Friedrichs
 Contributors:
 
@@ -403,9 +403,9 @@ class ForceField(object):
 
         # Set periodic boundary conditions.
 
-        boxSize = topology.getUnitCellDimensions()
-        if boxSize is not None:
-            sys.setDefaultPeriodicBoxVectors((boxSize[0], 0, 0), (0, boxSize[1], 0), (0, 0, boxSize[2]))
+        boxVectors = topology.getPeriodicBoxVectors()
+        if boxVectors is not None:
+            sys.setDefaultPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2])
         elif nonbondedMethod is not NoCutoff and nonbondedMethod is not CutoffNonPeriodic:
             raise ValueError('Requested periodic boundary conditions for a Topology that does not specify periodic box dimensions')
 

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

@@ -89,6 +89,18 @@ def _is_gro_box(line):
     else:
         return 0
 
+def _construct_box_vectors(line):
+    """Create the periodic box vectors based on the values stored in the file.
+
+    @param[in] line The line containing the description
+    """
+
+    sline = line.split()
+    values = [float(i) for i in sline]
+    if len(sline) == 3:
+        return (Vec3(values[0], 0, 0), Vec3(0, values[1], 0), Vec3(0, 0, values[2]))*nanometers
+    return (Vec3(values[0], values[3], values[4]), Vec3(values[5], values[1], values[6]), Vec3(values[7], values[8], values[2]))*nanometers
+
 class GromacsGroFile(object):
     """GromacsGroFile parses a Gromacs .gro file and constructs a set of atom positions from it.
 
@@ -140,7 +152,7 @@ class GromacsGroFile(object):
                 xyz.append(Vec3(pos[0], pos[1], pos[2]))
             elif _is_gro_box(line) and ln == na + 2:
                 sline = line.split()
-                boxes.append(tuple([float(i) for i in sline])*nanometers)
+                boxes.append(_construct_box_vectors(line))
                 xyzs.append(xyz*nanometers)
                 xyz = []
                 ln = -1
@@ -160,7 +172,7 @@ class GromacsGroFile(object):
         ## A list containing the name of the residue that each atom belongs to
         self.residueNames = resname
         self._positions = xyzs
-        self._unitCellDimensions = boxes
+        self._periodicBoxVectors = boxes
         self._numpyPositions = None
 
     def getNumFrames(self):
@@ -182,10 +194,21 @@ class GromacsGroFile(object):
             return self._numpyPositions[frame]
         return self._positions[frame]
 
+    def getPeriodicBoxVectors(self, frame=0):
+        """Get the vectors defining the periodic box.
+
+        Parameters:
+         - frame (int=0) the index of the frame for which to get the box vectors
+        """
+        return self._periodicBoxVectors[frame]
+
     def getUnitCellDimensions(self, frame=0):
         """Get the dimensions of the crystallographic unit cell.
 
         Parameters:
          - frame (int=0) the index of the frame for which to get the unit cell dimensions
         """
-        return self._unitCellDimensions[frame]
+        xsize = self._periodicBoxVectors[frame][0][0].value_in_unit(nanometers)
+        ysize = self._periodicBoxVectors[frame][1][1].value_in_unit(nanometers)
+        zsize = self._periodicBoxVectors[frame][2][2].value_in_unit(nanometers)
+        return Vec3(xsize, ysize, zsize)*nanometers

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

@@ -39,6 +39,7 @@ import simtk.unit as unit
 from .. import element
 import warnings
 import sys
+import math
 
 class PdbStructure(object):
     """
@@ -137,7 +138,7 @@ class PdbStructure(object):
         self._current_model = None
         self.default_model = None
         self.models_by_number = {}
-        self._unit_cell_dimensions = None
+        self._periodic_box_vectors = None
         self.sequences = []
         self.modified_residues = []
         # read file
@@ -170,7 +171,7 @@ class PdbStructure(object):
                 self._current_model._current_chain._add_ter_record()
                 self._reset_residue_numbers()
             elif (pdb_line.find("CRYST1") == 0):
-                self._unit_cell_dimensions = Vec3(float(pdb_line[6:15]), float(pdb_line[15:24]), float(pdb_line[24:33]))*unit.angstroms
+                self._compute_periodic_box_vectors(pdb_line)
             elif (pdb_line.find("CONECT") == 0):
                 atoms = [int(pdb_line[6:11])]
                 for pos in (11,16,21,26):
@@ -188,6 +189,29 @@ 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
@@ -283,9 +307,9 @@ class PdbStructure(object):
         for model in self.models:
             model._finalize()
 
-    def get_unit_cell_dimensions(self):
-        """Get the dimensions of the crystallographic unit cell (may be None)."""
-        return self._unit_cell_dimensions
+    def get_periodic_box_vectors(self):
+        """Get the vectors defining the crystallographic unit cell (may be None)."""
+        return self._periodic_box_vectors
 
 
 class Sequence(object):

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

@@ -94,7 +94,7 @@ class Modeller(object):
         # Copy over the existing model.
 
         newTopology = Topology()
-        newTopology.setUnitCellDimensions(deepcopy(self.topology.getUnitCellDimensions()))
+        newTopology.setUnitCellDimensions(self.topology.getUnitCellDimensions())
         newAtoms = {}
         newPositions = []*nanometer
         for chain in self.topology.chains():
@@ -140,7 +140,7 @@ class Modeller(object):
          - toDelete (list) a list of Atoms, Residues, Chains, and bonds (specified as tuples of Atoms) to delete
         """
         newTopology = Topology()
-        newTopology.setUnitCellDimensions(deepcopy(self.topology.getUnitCellDimensions()))
+        newTopology.setUnitCellDimensions(self.topology.getUnitCellDimensions())
         newAtoms = {}
         newPositions = []*nanometer
         deleteSet = set(toDelete)
@@ -189,7 +189,7 @@ class Modeller(object):
         else:
             raise ValueError('Unknown water model: %s' % model)
         newTopology = Topology()
-        newTopology.setUnitCellDimensions(deepcopy(self.topology.getUnitCellDimensions()))
+        newTopology.setUnitCellDimensions(self.topology.getUnitCellDimensions())
         newAtoms = {}
         newPositions = []*nanometer
         for chain in self.topology.chains():
@@ -472,7 +472,7 @@ class Modeller(object):
                     for atom2 in molAtoms:
                         if atom2.element == elem.hydrogen:
                             newTopology.addBond(atom1, atom2)
-        newTopology.setUnitCellDimensions(deepcopy(box)*nanometer)
+        newTopology.setUnitCellDimensions(box*nanometer)
         self.topology = newTopology
         self.positions = newPositions
 
@@ -610,7 +610,7 @@ class Modeller(object):
         # Loop over residues.
 
         newTopology = Topology()
-        newTopology.setUnitCellDimensions(deepcopy(self.topology.getUnitCellDimensions()))
+        newTopology.setUnitCellDimensions(self.topology.getUnitCellDimensions())
         newAtoms = {}
         newPositions = []*nanometer
         newIndices = []
@@ -883,7 +883,7 @@ class Modeller(object):
         # Create the new Topology.
 
         newTopology = Topology()
-        newTopology.setUnitCellDimensions(deepcopy(self.topology.getUnitCellDimensions()))
+        newTopology.setUnitCellDimensions(self.topology.getUnitCellDimensions())
         newAtoms = {}
         newPositions = []*nanometer
         for chain in self.topology.chains():

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

@@ -142,7 +142,7 @@ class PDBFile(object):
             self._positions.append(coords*nanometers)
         ## The atom positions read from the PDB file.  If the file contains multiple frames, these are the positions in the first frame.
         self.positions = self._positions[0]
-        self.topology.setUnitCellDimensions(pdb.get_unit_cell_dimensions())
+        self.topology.setPeriodicBoxVectors(pdb.get_periodic_box_vectors())
         self.topology.createStandardBonds()
         self.topology.createDisulfideBonds(self.positions)
         self._numpyPositions = None

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

@@ -33,6 +33,7 @@ __version__ = "1.0"
 
 import os
 import sys
+import math
 from simtk.openmm import Vec3
 from simtk.openmm.app.internal.pdbx.reader.PdbxReader import PdbxReader
 from simtk.openmm.app import Topology
@@ -145,8 +146,22 @@ class PDBxFile(object):
         cell = block.getObj('cell')
         if cell is not None and cell.getRowCount() > 0:
             row = cell.getRow(0)
-            cellSize = [float(row[cell.getAttributeIndex(attribute)]) for attribute in ('length_a', 'length_b', 'length_c')]*angstroms
-            self.topology.setUnitCellDimensions(cellSize)
+            (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)
 
         # Add bonds based on struct_conn records.
 

wrappers/python/simtk/openmm/app/topology.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-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

@@ -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

+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/systems/triclinic.gro

+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