Merge pull request #802 from swails/fix-charmm-box-vecs

Fix for unit cell vector calculation

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

@@ -39,6 +39,7 @@ import forcefield as ff
 import element as elem
 import simtk.unit as unit
 import simtk.openmm as mm
+from simtk.openmm.app.internal.unitcell import computePeriodicBoxVectors
 
 # Enumerated values for implicit solvent model
 
@@ -141,7 +142,8 @@ class AmberPrmtopFile(object):
         # Set the periodic box size.
 
         if prmtop.getIfBox():
-            top.setUnitCellDimensions(tuple(x.value_in_unit(unit.nanometer) for x in prmtop.getBoxBetaAndDimensions()[1:4])*unit.nanometer)
+            box = prmtop.getBoxBetaAndDimensions()
+            top.setPeriodicBoxVectors(computePeriodicBoxVectors(*(box[1:4] + box[0:1]*3)))
 
     def createSystem(self, nonbondedMethod=ff.NoCutoff, nonbondedCutoff=1.0*unit.nanometer,
                      constraints=None, rigidWater=True, implicitSolvent=None,

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

@@ -8,7 +8,7 @@ Structures at Stanford, funded under the NIH Roadmap for Medical Research,
 grant U54 GM072970. See https://simtk.org.  This code was originally part of
 the ParmEd program and was ported for use with OpenMM.
 
-Copyright (c) 2014 the Authors
+Copyright (c) 2014-2015 the Authors
 
 Author: Jason M. Swails
 Contributors:
@@ -44,6 +44,7 @@ from simtk.openmm.app import (forcefield as ff, Topology, element)
 from simtk.openmm.app.amberprmtopfile import HCT, OBC1, OBC2, GBn, GBn2
 from simtk.openmm.app.internal.customgbforces import (GBSAHCTForce,
                 GBSAOBC1Force, GBSAOBC2Force, GBSAGBnForce, GBSAGBn2Force, convertParameters)
+from simtk.openmm.app.internal.unitcell import computePeriodicBoxVectors
 # CHARMM imports
 from simtk.openmm.app.internal.charmm.topologyobjects import (
                 ResidueList, AtomList, TrackedList, Bond, Angle, Dihedral,
@@ -744,7 +745,7 @@ class CharmmPsfFile(object):
         Parameters:
             - a, b, c (floats) : Lengths of the periodic cell
             - alpha, beta, gamma (floats, optional) : Angles between the
-                periodic cells.
+                periodic cell vectors.
         """
         try:
             # Since we are setting the box, delete the cached box lengths if we
@@ -752,8 +753,7 @@ class CharmmPsfFile(object):
             del self._boxLengths
         except AttributeError:
             pass
-        self.box_vectors = _box_vectors_from_lengths_angles(a, b, c,
-                                                            alpha, beta, gamma)
+        self.box_vectors = computePeriodicBoxVectors(a, b, c, alpha, beta, gamma)
         # If we already have a _topology instance, then we have possibly changed
         # the existence of box information (whether or not this is a periodic
         # system), so delete any cached reference to a topology so it's
@@ -1523,58 +1523,6 @@ class CharmmPsfFile(object):
         """ Deletes the CMAP terms from the CHARMM PSF """
         self.cmap_list = TrackedList()
 
-def _box_vectors_from_lengths_angles(a, b, c, alpha, beta, gamma):
-    """
-    This method takes the lengths and angles from a unit cell and creates unit
-    cell vectors.
-
-    Parameters:
-        - a (unit, dimension length): Length of the first vector
-        - b (unit, dimension length): Length of the second vector
-        - c (unit, dimension length): Length of the third vector
-        - alpha (float): Angle between b and c in degrees
-        - beta (float): Angle between a and c in degrees
-        - gamma (float): Angle between a and b in degrees
-
-    Returns:
-        Tuple of box vectors (as Vec3 instances)
-   """
-    if not (u.is_quantity(a) and u.is_quantity(b) and u.is_quantity(c)):
-        raise TypeError('a, b, and c must be units of dimension length')
-    if u.is_quantity(alpha): alpha = alpha.value_in_unit(u.degree)
-    if u.is_quantity(beta): beta = beta.value_in_unit(u.degree)
-    if u.is_quantity(gamma): gamma = gamma.value_in_unit(u.degree)
-    a = a.value_in_unit(u.angstrom)
-    b = b.value_in_unit(u.angstrom)
-    c = c.value_in_unit(u.angstrom)
-
-    if alpha <= 2 * pi and beta <= 2 * pi and gamma <= 2 * pi:
-        raise ValueError('box angles must be given in degrees')
-
-    alpha *= pi / 180
-    beta *= pi / 180
-    gamma *= pi / 180
-
-    av = Vec3(a, 0.0, 0.0) * u.angstrom
-    bx = b * cos(gamma)
-    by = b * sin(gamma)
-    bz = 0.0
-    cx = c * cos(beta)
-    cy = c * (cos(alpha) - cos(beta) * cos(gamma))
-    cz = sqrt(c * c - cx * cx - cy * cy)
-   
-    # Make sure any components that are close to zero are set to zero exactly
-    if abs(bx) < TINY: bx = 0.0
-    if abs(by) < TINY: by = 0.0
-    if abs(cx) < TINY: cx = 0.0
-    if abs(cy) < TINY: cy = 0.0
-    if abs(cz) < TINY: cz = 0.0
-
-    bv = Vec3(bx, by, bz) * u.angstrom
-    cv = Vec3(cx, cy, cz) * u.angstrom
-
-    return (av, bv, cv)
-
 # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
 def set_molecules(atom_list):

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

@@ -52,6 +52,7 @@ except:
 import simtk.unit as units
 import simtk.openmm
 from simtk.openmm.app import element as elem
+from simtk.openmm.app.internal.unitcell import computePeriodicBoxVectors
 from simtk.openmm.vec3 import Vec3
 import customgbforces as customgb
 
@@ -691,7 +692,6 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
     >>> system = readAmberSystem(prmtop_filename)
 
     """
-
     if prmtop_filename is None and prmtop_loader is None:
         raise Exception("Must specify a filename or loader")
     if prmtop_filename is not None and prmtop_loader is not None:
@@ -709,9 +709,6 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
     if prmtop.getIfPert()>0:
         raise Exception("perturbation not currently supported")
 
-    if prmtop.getIfBox()>1:
-        raise Exception("only standard periodic boxes are currently supported")
-
     if prmtop.has_scee_scnb and (scee is not None or scnb is not None):
         warnings.warn("1-4 scaling parameters in topology file are being ignored. "
             "This is not recommended unless you know what you are doing.")
@@ -820,17 +817,7 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
         # System is periodic.
         # Set periodic box vectors for periodic system
         (boxBeta, boxX, boxY, boxZ) = prmtop.getBoxBetaAndDimensions()
-        boxBeta = boxBeta.value_in_unit(units.degrees)
-        boxX = boxX.value_in_unit(units.angstroms)
-        boxY = boxY.value_in_unit(units.angstroms)
-        boxZ = boxZ.value_in_unit(units.angstroms)
-        tmp = [[0.0, 0.0, 0.0],[0.0, 0.0, 0.0],[0.0, 0.0, 0.0]]
-        _box_vectors_from_lengths_angles([boxX, boxY, boxZ],
-                                         [boxBeta, boxBeta, boxBeta],
-                                         tmp)
-        xVec = units.Quantity(tmp[0], units.angstroms)
-        yVec = units.Quantity(tmp[1], units.angstroms)
-        zVec = units.Quantity(tmp[2], units.angstroms)
+        xVec, yVec, zVec = computePeriodicBoxVectors(boxX, boxY, boxZ, boxBeta, boxBeta, boxBeta)
         system.setDefaultPeriodicBoxVectors(xVec, yVec, zVec)
 
         # Set cutoff.
@@ -1286,10 +1273,10 @@ class AmberAsciiRestart(object):
             except (IndexError, ValueError):
                 raise ValueError('Could not parse box line in %s' %
                                  self.filename)
-            lengths = tmp[:3]
-            angles = tmp[3:]
-            _box_vectors_from_lengths_angles(lengths, angles, boxVectors)
-            self.boxVectors = [units.Quantity(Vec3(*x), units.angstrom) for x in boxVectors]
+            lengths = tmp[:3] * units.angstroms
+            angles = tmp[3:] * units.degrees
+            self.boxVectors = computePeriodicBoxVectors(lengths[0], lengths[1],
+                    lengths[2], angles[0], angles[1], angles[2])
 
 class AmberNetcdfRestart(object):
     """
@@ -1360,11 +1347,12 @@ class AmberNetcdfRestart(object):
             if ('cell_lengths' in ncfile.variables and
                 'cell_angles' in ncfile.variables):
                 self.boxVectors = np.zeros((3,3), np.float32)
-                _box_vectors_from_lengths_angles(
-                        ncfile.variables['cell_lengths'][:],
-                        ncfile.variables['cell_angles'][:],
-                        self.boxVectors,
-                )
+                leng = units.Quantity(ncfile.variables['cell_lengths'][:],
+                        units.angstroms)
+                angl = units.Quantity(ncfile.variables['cell_angles'][:],
+                        units.degrees)
+                self.boxVectors = computePeriodicBoxVectors(leng[0], leng[1],
+                        leng[2], angl[0], angl[1], angl[2])
             if 'time' in ncfile.variables:
                 self.time = ncfile.variables['time'].getValue()
         finally:
@@ -1375,52 +1363,18 @@ class AmberNetcdfRestart(object):
             self.coordinates = [Vec3(*x) for x in self.coordinates]
             if self.velocities is not None:
                 self.velocities = [Vec3(*x) for x in self.velocities]
+        else:
             if self.boxVectors is not None:
-                self.boxVectors = [Vec3(*x) for x in self.boxVectors]
+                self.boxVectors = np.asarray(self.boxVectors.value_in_unit(units.nanometers))
+                self.boxVectors = units.Quantity(self.boxVectors, units.nanometers)
 
         # Now add the units
         self.coordinates = units.Quantity(self.coordinates, units.angstroms)
         if self.velocities is not None:
             self.velocities = units.Quantity(self.velocities,
                                              units.angstroms/units.picoseconds)
-        if self.boxVectors is not None:
-            self.boxVectors = [units.Quantity(x, units.angstroms) for x in self.boxVectors]
         self.time = units.Quantity(self.time, units.picosecond)
 
-def _box_vectors_from_lengths_angles(lengths, angles, boxVectors):
-    """
-    Converts lengths and angles into a series of box vectors and modifies
-    boxVectors in-place (it must be a mutable sequence)
-
-    Parameters
-    ----------
-    lengths : 3-element array of floats
-        Lengths of the 3 periodic box vectors
-    angles : 3-element array of floats
-        Angles (in degrees) between the 3 periodic box vectors
-    boxVectors : mutable 3x3 sequence
-    """
-    alpha = angles[0] * pi / 180.0
-    beta = angles[1] * pi / 180.0
-    gamma = angles[2] * pi / 180.0
-
-    boxVectors[0][0] = lengths[0]
-
-    boxVectors[1][0] = lengths[1] * cos(gamma)
-    boxVectors[1][1] = lengths[1] * sin(gamma)
-
-    boxVectors[2][0] = cx = lengths[2] * cos(beta)
-    boxVectors[2][1] = cy = lengths[2] * (cos(alpha) - cos(beta) * cos(gamma))
-    boxVectors[2][2] = sqrt(lengths[2]*lengths[2] - cx*cx - cy*cy)
-
-    boxVectors[0][1] = boxVectors[0][2] = boxVectors[1][2] = 0.0
-
-    # Now make sure any vector close to zero is zero exactly
-    for i in range(3):
-        for j in range(3):
-            if abs(boxVectors[i][j]) < TINY:
-                boxVectors[i][j] = 0.0
-
 def readAmberCoordinates(filename, asNumpy=False):
     """
     Read atomic coordinates (and optionally, box vectors) from Amber formatted coordinate file.

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

@@ -32,16 +32,24 @@ __author__ = "Peter Eastman"
 __version__ = "1.0"
 
 from simtk.openmm import Vec3
-from simtk.unit import nanometers, is_quantity, norm, dot
+from simtk.unit import nanometers, is_quantity, norm, dot, radians
 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.
+    Lengths should be given in nanometers and angles in radians (or as Quantity
+    instances)
     """
 
+    if is_quantity(a_length): a_length = a_length.value_in_unit(nanometers)
+    if is_quantity(b_length): b_length = b_length.value_in_unit(nanometers)
+    if is_quantity(c_length): c_length = c_length.value_in_unit(nanometers)
+    if is_quantity(alpha): alpha = alpha.value_in_unit(radians)
+    if is_quantity(beta): beta = beta.value_in_unit(radians)
+    if is_quantity(gamma): gamma = gamma.value_in_unit(radians)
+
     # Compute the vectors.
 
     a = [a_length, 0, 0]
@@ -76,8 +84,10 @@ def computeLengthsAndAngles(periodicBoxVectors):
 
     Lengths are returned in nanometers and angles in radians.
     """
-
-    (a, b, c) = vectors.value_in_unit(nanometers)
+    if is_quantity(periodicBoxVectors):
+        (a, b, c) = vectors.value_in_unit(nanometers)
+    else:
+        a, b, c = vectors
     a_length = norm(a)
     b_length = norm(b)
     c_length = norm(c)

wrappers/python/tests/TestAmberInpcrdFile.py

@@ -60,5 +60,29 @@ class TestAmberInpcrdFile(unittest.TestCase):
         self.assertTrue(inpcrd.boxVectors is None)
         self.assertTrue(inpcrd.velocities is None)
 
+    def test_CrdBoxTruncoct(self):
+        # Check that the box vectors come out correct.
+        inpcrd = AmberInpcrdFile('systems/tz2.truncoct.rst7')
+        ac = Vec3(42.4388485, 0.0, 0.0) * angstroms
+        bc = Vec3(-14.146281691908937, 40.011730483685835, 0.0) * angstroms
+        cc = Vec3(-14.146281691908937, -20.0058628205162, 34.651176446201672) * angstroms
+        a, b, c = inpcrd.getBoxVectors()
+        diffa = ac - a
+        diffb = bc - b
+        diffc = cc - c
+        self.assertAlmostEqual(norm(diffa)/angstroms, 0)
+        self.assertAlmostEqual(norm(diffb)/angstroms, 0)
+        self.assertAlmostEqual(norm(diffc)/angstroms, 0)
+        # Make sure angles and lengths come out about right
+        la = norm(a).in_units_of(angstroms)
+        lb = norm(b).in_units_of(angstroms)
+        lc = norm(c).in_units_of(angstroms)
+        self.assertAlmostEqual(la/angstroms, 42.4388485, 6)
+        self.assertAlmostEqual(lb/angstroms, 42.4388485, 6)
+        self.assertAlmostEqual(lc/angstroms, 42.4388485, 6)
+        self.assertAlmostEqual(dot(a,b)/la/lb, cos(109.4712190*degrees), 6)
+        self.assertAlmostEqual(dot(a,c)/la/lc, cos(109.4712190*degrees), 6)
+        self.assertAlmostEqual(dot(b,c)/lc/lb, cos(109.4712190*degrees), 6)
+
 if __name__ == '__main__':
     unittest.main()

wrappers/python/tests/TestAmberPrmtopFile.py

@@ -9,6 +9,7 @@ prmtop1 = AmberPrmtopFile('systems/alanine-dipeptide-explicit.prmtop')
 prmtop2 = AmberPrmtopFile('systems/alanine-dipeptide-implicit.prmtop')
 prmtop3 = AmberPrmtopFile('systems/ff14ipq.parm7')
 prmtop4 = AmberPrmtopFile('systems/Mg_water.prmtop')
+prmtop5 = AmberPrmtopFile('systems/tz2.truncoct.parm7')
 inpcrd3 = AmberInpcrdFile('systems/ff14ipq.rst7')
 inpcrd4 = AmberInpcrdFile('systems/Mg_water.inpcrd')
 
@@ -229,7 +230,32 @@ class TestAmberPrmtopFile(unittest.TestCase):
         # Make sure the energy is relatively close to the value we get with
         # Amber using this force field.
         self.assertAlmostEqual(-7307.2735621/ene, 1, places=3)
-    
+
+    def test_triclinicParm(self):
+        """ Check that triclinic unit cells work correctly """
+        system = prmtop5.createSystem(nonbondedMethod=PME)
+        refa = Vec3(4.48903851, 0.0, 0.0) * nanometer
+        refb = Vec3(-1.4963460492639706, 4.232306137924705, 0.0) * nanometer
+        refc = Vec3(-1.4963460492639706, -2.116152812842565, 3.6652847799064165) * nanometer
+        a, b, c = system.getDefaultPeriodicBoxVectors()
+        la = norm(a)
+        lb = norm(b)
+        lc = norm(c)
+        diffa = a - refa
+        diffb = b - refb
+        diffc = c - refc
+        self.assertAlmostEqual(norm(diffa)/nanometers, 0)
+        self.assertAlmostEqual(norm(diffb)/nanometers, 0)
+        self.assertAlmostEqual(norm(diffc)/nanometers, 0)
+        self.assertAlmostEqual(dot(a, b)/la/lb, cos(109.4712190*degrees))
+        self.assertAlmostEqual(dot(a, c)/la/lc, cos(109.4712190*degrees))
+        self.assertAlmostEqual(dot(c, b)/lc/lb, cos(109.4712190*degrees))
+        self.assertAlmostEqual(la/nanometers, 4.48903851)
+        self.assertAlmostEqual(lb/nanometers, 4.48903851)
+        self.assertAlmostEqual(lc/nanometers, 4.48903851)
+        # Now make sure that the context builds correctly; then we can bail
+        self.assertTrue(Context(system, VerletIntegrator(1*femtoseconds)))
+
     def test_ImplicitSolventForces(self):
         """Compute forces for different implicit solvent types, and compare them to ones generated with a previous version of OpenMM to ensure they haven't changed."""