More support for triclinic boxes

docs-source/usersguide/application.rst

@@ -572,7 +572,7 @@ with the name :file:`simulateGromacs.py`.
         from sys import stdout
 
         gro = GromacsGroFile('input.gro')
-        top = GromacsTopFile('input.top', unitCellDimensions=gro.getUnitCellDimensions(),
+        top = GromacsTopFile('input.top', periodicBoxVectors=gro.getPeriodicBoxVectors(),
                 includeDir='/usr/local/gromacs/share/gromacs/top')
         system = top.createSystem(nonbondedMethod=PME, nonbondedCutoff=1*nanometer,
                 constraints=HBonds)
@@ -593,10 +593,10 @@ This script is nearly identical to the previous one, just replacing
 :class:`AmberInpcrdFile` and :class:`AmberPrmtopFile` with :class:`GromacsGroFile` and :class:`GromacsTopFile`.
 Note that when we create the :class:`GromacsTopFile`, we specify values for two extra
 options.  First, we specify
-:code:`unitCellDimensions=gro.getUnitCellDimensions()`\ .  Unlike OpenMM and
-AMBER, which can store periodic unit cell dimensions with the topology, Gromacs
-only stores them with the coordinates.  To let :class:`GromacsTopFile` create a :class:`Topology`
-object, we therefore need to tell it the unit cell dimensions that were loaded
+:code:`periodicBoxVectors=gro.getPeriodicBoxVectors()`\ .  Unlike OpenMM and
+AMBER, which can store periodic unit cell information with the topology, Gromacs
+only stores it with the coordinates.  To let :class:`GromacsTopFile` create a :class:`Topology`
+object, we therefore need to tell it the periodic box vectors that were loaded
 from the :file:`gro` file.  You only need to do this if you are simulating a periodic
 system.  For implicit solvent simulations, it usually can be omitted.
 

examples/simulateGromacs.py

@@ -4,7 +4,7 @@ from simtk.unit import *
 from sys import stdout
 
 gro = GromacsGroFile('input.gro')
-top = GromacsTopFile('input.top', unitCellDimensions=gro.getUnitCellDimensions(), includeDir='/usr/local/gromacs/share/gromacs/top')
+top = GromacsTopFile('input.top', periodicBoxVectors=gro.getPeriodicBoxVectors(), includeDir='/usr/local/gromacs/share/gromacs/top')
 system = top.createSystem(nonbondedMethod=PME, nonbondedCutoff=1*nanometer, constraints=HBonds)
 integrator = LangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
 simulation = Simulation(top.topology, system, integrator)

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

@@ -77,14 +77,17 @@ class DCDFile(object):
         header += struct.pack('<4i', 164, 4, len(list(topology.atoms())), 4)
         file.write(header)
 
-    def writeModel(self, positions, unitCellDimensions=None):
+    def writeModel(self, positions, unitCellDimensions=None, periodicBoxVectors=None):
         """Write out a model to the DCD file.
 
+        The periodic box can be specified either by the unit cell dimensions (for a rectangular box), or the full set of box
+        vectors (for an arbitrary triclinic box).  If neither is specified, the box vectors specified in the Topology will be
+        used.  Regardless of the value specified, no dimensions will be written if the Topology does not represent a periodic system.
+
         Parameters:
          - positions (list) The list of atomic positions to write
-         - unitCellDimensions (Vec3=None) The dimensions of the crystallographic unit cell.  If None, the dimensions specified in
-           the Topology will be used.  Regardless of the value specified, no dimensions will be written if the Topology does not
-           represent a periodic system.
+         - unitCellDimensions (Vec3=None) The dimensions of the crystallographic unit cell.
+         - periodicBoxVectors (tuple of Vec3=None) The vectors defining the periodic box.
         """
         if len(list(self._topology.atoms())) != len(positions):
             raise ValueError('The number of positions must match the number of atoms')
@@ -107,12 +110,22 @@ class DCDFile(object):
         # Write the data.
 
         file.seek(0, os.SEEK_END)
-        boxSize = self._topology.getUnitCellDimensions()
-        if boxSize is not None:
-            if unitCellDimensions is not None:
-                boxSize = unitCellDimensions
-            size = boxSize.value_in_unit(angstroms)
-            file.write(struct.pack('<i6di', 48, size[0], 0, size[1], 0, 0, size[2], 48))
+        boxVectors = self._topology.getPeriodicBoxVectors()
+        if boxVectors is not None:
+            if getPeriodicBoxVectors is not None:
+                boxVectors = getPeriodicBoxVectors
+            elif unitCellDimensions is not None:
+                if is_quantity(unitCellDimensions):
+                    unitCellDimensions = unitCellDimensions.value_in_unit(nanometers)
+                boxVectors = (Vec3(unitCellDimensions[0], 0, 0), Vec3(0, unitCellDimensions[1], 0), Vec3(0, 0, unitCellDimensions[2]))*nanometers
+            (a_length, b_length, c_length, alpha, beta, gamma) = computeLengthsAndAngles(boxVectors)
+            a_length = a_length.value_in_unit(angstroms)
+            b_length = b_length.value_in_unit(angstroms)
+            c_length = c_length.value_in_unit(angstroms)
+            angle1 = math.sin(math.pi/2-gamma)
+            angle2 = math.sin(math.pi/2-beta)
+            angle3 = math.sin(math.pi/2-alpha)
+            file.write(struct.pack('<i6di', 48, a_length, angle1, b_length, angle2, angle3, c_length, 48))
         length = struct.pack('<i', 4*len(positions))
         for i in range(3):
             file.write(length)

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

@@ -37,42 +37,11 @@ __version__ = "1.0"
 from simtk.openmm.vec3 import Vec3
 import simtk.unit as unit
 from .. import element
+from unitcell import computePeriodicBoxVectors
 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.
@@ -203,9 +172,9 @@ class PdbStructure(object):
                 self._current_model._current_chain._add_ter_record()
                 self._reset_residue_numbers()
             elif (pdb_line.find("CRYST1") == 0):
-                a_length = float(pdb_line[6:15])
-                b_length = float(pdb_line[15:24])
-                c_length = float(pdb_line[24:33])
+                a_length = float(pdb_line[6:15])*0.1
+                b_length = float(pdb_line[15:24])*0.1
+                c_length = float(pdb_line[24:33])*0.1
                 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

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

+"""
+unitcell.py: Routines for converting between different representations of the periodic unit cell.
+
+This is part of the OpenMM molecular simulation toolkit originating from
+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) 2015 Stanford University and the Authors.
+Authors: Peter Eastman
+Contributors:
+
+Permission is hereby granted, free of charge, to any person obtaining a
+copy of this software and associated documentation files (the "Software"),
+to deal in the Software without restriction, including without limitation
+the rights to use, copy, modify, merge, publish, distribute, sublicense,
+and/or sell copies of the Software, and to permit persons to whom the
+Software is furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
+DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
+OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
+USE OR OTHER DEALINGS IN THE SOFTWARE.
+"""
+__author__ = "Peter Eastman"
+__version__ = "1.0"
+
+from simtk.openmm import Vec3
+from simtk.unit import nanometers, is_quantity, norm, dot
+import math
+
+
+def computePeriodicBoxVectors(a_length, b_length, c_length, alpha, beta, gamma):
+    """Convert lengths and angles to periodic box vectors.
+    
+    Lengths should be given in nanometers and angles in radians.
+    """
+
+    # 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)*nanometers
+
+def computeLengthsAndAngles(periodicBoxVectors):
+    """Convert periodic box vectors to lengths and angles.
+
+    Lengths are returned in nanometers and angles in radians.
+    """
+
+    (a, b, c) = vectors.value_in_unit(nanometers)
+    a_length = norm(a)
+    b_length = norm(b)
+    c_length = norm(c)
+    alpha = math.acos(dot(b, c)/(b_length*c_length))
+    beta = math.acos(dot(c, a)/(c_length*a_length))
+    gamma = math.acos(dot(a, b)/(a_length*b_length))
+    return (a_length, b_length, c_length, alpha, beta, gamma)

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

@@ -36,7 +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.internal.unitcell import computePeriodicBoxVectors
 from simtk.openmm.app import Topology
 from simtk.unit import nanometers, angstroms, is_quantity, norm, Quantity
 import element as elem
@@ -147,7 +147,7 @@ class PDBxFile(object):
         cell = block.getObj('cell')
         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')]
+            (a, b, c) = [float(row[cell.getAttributeIndex(attribute)])*0.1 for attribute in ('length_a', 'length_b', 'length_c')]
             (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))
 

wrappers/python/tests/TestPdbFile.py

@@ -3,6 +3,7 @@ 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"""
@@ -39,6 +40,26 @@ class TestPdbFile(unittest.TestCase):
                 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)