Beginnings of improvements for amber_file_parser and the CHARMM PSF file parser

to support triclinic cells.

Also extend unitcell.computePeriodicBoxVectors to take Quantity inputs rather
than *just* scalars assumed to be nanometers.

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

@@ -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)) / sin(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

@@ -820,17 +820,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 +1276,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 +1350,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 +1366,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

@@ -39,9 +39,17 @@ 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 u.is_quantity(a_length): a_length = a_length.value_in_unit(u.nanometers)
+    if u.is_quantity(b_length): a_length = a_length.value_in_unit(u.nanometers)
+    if u.is_quantity(c_length): a_length = a_length.value_in_unit(u.nanometers)
+    if u.is_quantity(alpha): alpha = alpha.value_in_unit(u.radians)
+    if u.is_quantity(beta): beta = beta.value_in_unit(u.radians)
+    if u.is_quantity(gamma): gamma = gamma.value_in_unit(u.radians)
+
     # Compute the vectors.
 
     a = [a_length, 0, 0]