Merge pull request #1076 from swails/tests

Add a set of unit tests for unit handling in the core API.

devtools/forcefield-scripts/processTinkerForceField.py

@@ -1017,7 +1017,6 @@ if( isAmoeba ):
     torsionTorsionUnit      = 1.0
     outputString            = """ <AmoebaTorsionTorsionForce >"""
     tinkerXmlFile.write( "%s\n" % (outputString ) )
-    conversion              = 41.84/radian
     torsionTorsions         = forces['tortors']
     for (index, torsionTorsion) in enumerate(torsionTorsions):
        torInfo       = torsionTorsion[0]

plugins/amoeba/openmmapi/include/openmm/AmoebaAngleForce.h

@@ -127,7 +127,7 @@ public:
      * @param particle1     the index of the first particle connected by the angle
      * @param particle2     the index of the second particle connected by the angle
      * @param particle3     the index of the third particle connected by the angle
-     * @param length        the angle measured in degrees
+     * @param length        the equilibrium angle, measured in degrees
      * @param quadratic k   the quadratic force constant for the angle, measured in kJ/mol/radian^2
      * @return the index of the angle that was added
      */
@@ -140,7 +140,7 @@ public:
      * @param particle1     the index of the first particle connected by the angle
      * @param particle2     the index of the second particle connected by the angle
      * @param particle3     the index of the third particle connected by the angle
-     * @param length        the equilibrium angle, measured in degress
+     * @param length        the equilibrium angle, measured in degrees
      * @param quadratic k   the quadratic force constant for the angle, measured in kJ/mol/radian^2
      */
     void getAngleParameters(int index, int& particle1, int& particle2, int& particle3, double& length, double& quadraticK) const;

plugins/amoeba/openmmapi/include/openmm/AmoebaOutOfPlaneBendForce.h

@@ -74,7 +74,7 @@ public:
      * 
      * @return global cubicK term
      */
-    double getAmoebaGlobalOutOfPlaneBendCubic(void) const;
+    double getAmoebaGlobalOutOfPlaneBendCubic() const;
 
     /** 
      * Set the global cubic term
@@ -88,7 +88,7 @@ public:
      * 
      * @return global  quartic term
      */
-    double getAmoebaGlobalOutOfPlaneBendQuartic(void) const;
+    double getAmoebaGlobalOutOfPlaneBendQuartic() const;
 
     /** 
      * Set the global pentic term
@@ -102,7 +102,7 @@ public:
      * 
      * @return global penticK term
      */
-    double getAmoebaGlobalOutOfPlaneBendPentic(void) const;
+    double getAmoebaGlobalOutOfPlaneBendPentic() const;
 
     /** 
      * Set the global sextic term
@@ -116,7 +116,7 @@ public:
      * 
      * @return global sexticK term
      */
-    double getAmoebaGlobalOutOfPlaneBendSextic(void) const;
+    double getAmoebaGlobalOutOfPlaneBendSextic() const;
 
     /**
      * Add an out-of-plane bend term to the force field.

wrappers/python/src/swig_doxygen/swigInputBuilder.py

@@ -10,6 +10,7 @@ import getopt
 import re
 import xml.etree.ElementTree as etree
 from distutils.version import LooseVersion
+import copy
 
 try:
     from html.parser import HTMLParser
@@ -20,7 +21,6 @@ except ImportError:
 INDENT = "   "
 docTags = {'emphasis':'i', 'bold':'b', 'itemizedlist':'ul', 'listitem':'li', 'preformatted':'pre', 'computeroutput':'tt', 'subscript':'sub'}
 
-
 def striphtmltags(s):
     """Strip a couple html tags used inside docstrings in the C++ source
     to produce something more easily read as plain text.
@@ -157,7 +157,7 @@ class SwigInputBuilder:
                  skipAdditionalMethods=[],
                  SWIG_VERSION='3.0.2'):
         self.nodeByID={}
-        self.SWIG_COMPACT_ARGUMENTS = LooseVersion(SWIG_VERSION) < LooseVersion('3.0.6')
+        self.SWIG_COMPACT_ARGUMENTS = LooseVersion(SWIG_VERSION) < LooseVersion('3.0.5')
 
         self.configModule = __import__(os.path.splitext(configFilename)[0])
 
@@ -477,7 +477,7 @@ class SwigInputBuilder:
                     valueUnits=[None, ()]
 
                 index=0
-                if valueUnits[0]:
+                if valueUnits[0] is not None:
                     sys.stdout.write("%s.%s() returns %s\n" %
                                      (shortClassName, methName, valueUnits[0]))
                     if len(valueUnits[1])>0:
@@ -491,10 +491,10 @@ class SwigInputBuilder:
                                  % (addText, INDENT, valueUnits[0])
 
                 for vUnit in valueUnits[1]:
-                        if vUnit!=None:
-                            addText = "%s%sval[%s]=unit.Quantity(val[%s], %s)\n" \
+                    if vUnit is not None:
+                        addText = "%s%sval[%s]=unit.Quantity(val[%s], %s)\n" \
                                      % (addText, INDENT, index, index, vUnit)
-                        index+=1
+                    index+=1
 
                 if key in self.configModule.STEAL_OWNERSHIP:
                     for argNum in self.configModule.STEAL_OWNERSHIP[key]:

wrappers/python/src/swig_doxygen/swigInputConfig.py

@@ -216,22 +216,22 @@ UNITS = {
 ("AmoebaGeneralizedKirkwoodForce",       "getDielectricOffset")                           :  ( 'unit.nanometer', ()),
 ("AmoebaGeneralizedKirkwoodForce",       "getIncludeCavityTerm")                          :  ( None,()),
 ("AmoebaGeneralizedKirkwoodForce",       "getProbeRadius")                                :  ( 'unit.nanometer', ()),
-("AmoebaGeneralizedKirkwoodForce",       "getSurfaceAreaFactor")                          :  ( '(unit.nanometer*unit.nanometer)/unit.kilojoule_per_mole',()),
+("AmoebaGeneralizedKirkwoodForce",       "getSurfaceAreaFactor")                          :  ( 'unit.kilojoule_per_mole/(unit.nanometer*unit.nanometer)',()),
 
-("AmoebaAngleForce",             "getAmoebaGlobalAngleCubic")             :  ( None,()),
-("AmoebaAngleForce",             "getAmoebaGlobalAngleQuartic")           :  ( None,()),
-("AmoebaAngleForce",             "getAmoebaGlobalAnglePentic")            :  ( None,()),
-("AmoebaAngleForce",             "getAmoebaGlobalAngleSextic")            :  ( None,()),
-("AmoebaAngleForce",             "getAngleParameters")                            :  ( None, (None, None, None, 'unit.radian', 'unit.kilojoule_per_mole/(unit.radian*unit.radian)')),
+("AmoebaAngleForce",             "getAmoebaGlobalAngleCubic")             :  ( '1/unit.radian',()),
+("AmoebaAngleForce",             "getAmoebaGlobalAngleQuartic")           :  ( '1/unit.radian**2',()),
+("AmoebaAngleForce",             "getAmoebaGlobalAnglePentic")            :  ( '1/unit.radian**3',()),
+("AmoebaAngleForce",             "getAmoebaGlobalAngleSextic")            :  ( '1/unit.radian**4',()),
+("AmoebaAngleForce",             "getAngleParameters")                            :  ( None, (None, None, None, 'unit.degree', 'unit.kilojoule_per_mole/(unit.radian*unit.radian)')),
 
-("AmoebaBondForce",              "getAmoebaGlobalBondCubic")              :  ( None,()),
-("AmoebaBondForce",              "getAmoebaGlobalBondQuartic")            :  ( None,()),
+("AmoebaBondForce",              "getAmoebaGlobalBondCubic")              :  ( '1/unit.nanometer',()),
+("AmoebaBondForce",              "getAmoebaGlobalBondQuartic")            :  ( '1/unit.nanometer**2',()),
 ("AmoebaBondForce",              "getBondParameters")                             :  ( None, (None, None, 'unit.nanometer', 'unit.kilojoule_per_mole/(unit.nanometer*unit.nanometer)')),
 
-("AmoebaInPlaneAngleForce",      "getAmoebaGlobalInPlaneAngleCubic")      :  ( None,()),
-("AmoebaInPlaneAngleForce",      "getAmoebaGlobalInPlaneAngleQuartic")    :  ( None,()),
-("AmoebaInPlaneAngleForce",      "getAmoebaGlobalInPlaneAnglePentic")     :  ( None,()),
-("AmoebaInPlaneAngleForce",      "getAmoebaGlobalInPlaneAngleSextic")     :  ( None,()),
+("AmoebaInPlaneAngleForce",      "getAmoebaGlobalInPlaneAngleCubic")      :  ( '1/unit.radian',()),
+("AmoebaInPlaneAngleForce",      "getAmoebaGlobalInPlaneAngleQuartic")    :  ( '1/unit.radian**2',()),
+("AmoebaInPlaneAngleForce",      "getAmoebaGlobalInPlaneAnglePentic")     :  ( '1/unit.radian**3',()),
+("AmoebaInPlaneAngleForce",      "getAmoebaGlobalInPlaneAngleSextic")     :  ( '1/unit.radian**4',()),
 ("AmoebaInPlaneAngleForce",      "getAngleParameters")                            :  ( None, (None, None, None, None, 'unit.radian', 'unit.kilojoule_per_mole/(unit.radian*unit.radian)')),
 
 ("AmoebaMultipoleForce",                 "getNumMultipoles")                              :  ( None,()),
@@ -256,7 +256,9 @@ UNITS = {
 #    void getCovalentMap(int index, CovalentType typeId, std::vector<int>& covalentAtoms )
 #    void getCovalentMaps(int index, std::vector < std::vector<int> >& covalentLists )
 
-("AmoebaMultipoleForce",                 "getMultipoleParameters")                        :  ( None, ()),
+("AmoebaMultipoleForce",                 "getMultipoleParameters")                        :  ( None, ('unit.elementary_charge', 'unit.elementary_charge/unit.nanometer',
+                                                                                                      'unit.elementary_charge/unit.nanometer**2', None, None, None, None, None, None,
+                                                                                                      'unit.nanometer**3')),
 ("AmoebaMultipoleForce",                 "getCovalentMap")                                :  ( None, ()),
 ("AmoebaMultipoleForce",                 "getCovalentMaps")                               :  ( None, ()),
 ("AmoebaMultipoleForce",                 "getScalingDistanceCutoff")                      :  ( 'unit.nanometer', ()),
@@ -268,17 +270,17 @@ UNITS = {
 ("AmoebaMultipoleForce",                 "getSystemMultipoleMoments")                     :  ( None, ()),
 
 ("AmoebaOutOfPlaneBendForce",            "getNumOutOfPlaneBends")                         :  ( None, ()),
-("AmoebaOutOfPlaneBendForce",            "getAmoebaGlobalOutOfPlaneBendCubic")            :  ( None,()),
-("AmoebaOutOfPlaneBendForce",            "getAmoebaGlobalOutOfPlaneBendQuartic")          :  ( None,()),
-("AmoebaOutOfPlaneBendForce",            "getAmoebaGlobalOutOfPlaneBendPentic")           :  ( None,()),
-("AmoebaOutOfPlaneBendForce",            "getAmoebaGlobalOutOfPlaneBendSextic")           :  ( None,()),
-("AmoebaOutOfPlaneBendForce",            "getOutOfPlaneBendParameters")                   :  ( None, (None, None, None, None, 'unit.kilojoule_per_mole')),
+("AmoebaOutOfPlaneBendForce",            "getAmoebaGlobalOutOfPlaneBendCubic")            :  ( '1/unit.radian',()),
+("AmoebaOutOfPlaneBendForce",            "getAmoebaGlobalOutOfPlaneBendQuartic")          :  ( '1/unit.radian**2',()),
+("AmoebaOutOfPlaneBendForce",            "getAmoebaGlobalOutOfPlaneBendPentic")           :  ( '1/unit.radian**3',()),
+("AmoebaOutOfPlaneBendForce",            "getAmoebaGlobalOutOfPlaneBendSextic")           :  ( '1/unit.radian**4',()),
+("AmoebaOutOfPlaneBendForce",            "getOutOfPlaneBendParameters")                   :  ( None, (None, None, None, None, 'unit.kilojoule_per_mole/unit.radians**2')),
 
 ("AmoebaPiTorsionForce",                  "getNumPiTorsions")                              :  ( None, ()),
 ("AmoebaPiTorsionForce",                  "getPiTorsionParameters")                        :  ( None, (None, None, None, None, None,  None, 'unit.kilojoule_per_mole')),
 
 ("AmoebaStretchBendForce",                "getNumStretchBends")                            :  ( None, ()),
-("AmoebaStretchBendForce",                "getStretchBendParameters")                      :  ( None, (None, None, None, 'unit.nanometer', 'unit.nanometer', 'unit.radian', 'unit.kilojoule_per_mole/unit.nanometer/unit.degree', 'unit.kilojoule_per_mole/unit.nanometer/unit.degree')),
+("AmoebaStretchBendForce",                "getStretchBendParameters")                      :  ( None, (None, None, None, 'unit.nanometer', 'unit.nanometer', 'unit.radian', 'unit.kilojoule_per_mole/unit.nanometer/unit.radian', 'unit.kilojoule_per_mole/unit.nanometer/unit.radian')),
 
 ("AmoebaTorsionTorsionForce",             "getNumTorsionTorsions")                         :  ( None, ()),
 ("AmoebaTorsionTorsionForce",             "getNumTorsionTorsionGrids")                     :  ( None, ()),
@@ -376,6 +378,9 @@ UNITS = {
 ("CustomTorsionForce", "getPerTorsionParameterName") : (None, ()),
 ("CustomTorsionForce", "getGlobalParameterName") : (None, ()),
 ("CustomTorsionForce", "getTorsionParameters") : (None, ()),
+("DrudeForce", "getParticleParameters") : (None, (None, None, None, None, None, 'unit.elementary_charge', 'unit.nanometer**3', None, None)),
+("DrudeForce", "getNumScreenedPairs") : (None, ()),
+("DrudeForce", "getScreenedPairParameters") : (None, ()),
 ("GBSAOBCForce", "getParticleParameters")
  : (None, ('unit.elementary_charge',
            'unit.nanometer', None)),

wrappers/python/src/swig_doxygen/swig_lib/python/pythoncode.i

@@ -2,9 +2,10 @@
 
 try:
     import numpy
-except:
-    pass
+except ImportError:
+    numpy = None
 
+import copy
 import sys
 import math
 RMIN_PER_SIGMA=math.pow(2, 1/6.0)
@@ -224,3 +225,70 @@ class State(_object):
     self._system = args[0]
     self._integrator = args[1]
 %}
+
+%pythonprepend OpenMM::AmoebaAngleForce::addAngle %{
+    try:
+        length = args[3]
+        if isinstance(args, tuple):
+            args = list(args)
+    except (NameError, UnboundLocalError):
+        if unit.is_quantity(length):
+            length = length.value_in_unit(unit.degree)
+    else:
+        if unit.is_quantity(length):
+            args[3] = length.value_in_unit(unit.degree)
+%}
+
+%pythonprepend OpenMM::AmoebaAngleForce::setAngleParameters %{
+    try:
+        length = args[4]
+        if isinstance(args, tuple):
+            args = list(args)
+    except (NameError, UnboundLocalError):
+        if unit.is_quantity(length):
+            length = length.value_in_unit(unit.degree)
+    else:
+        if unit.is_quantity(length):
+            args[4] = length.value_in_unit(unit.degree)
+%}
+
+%pythonprepend OpenMM::AmoebaTorsionTorsionForce::setTorsionTorsionGrid %{
+    def deunitize_grid(grid):
+        if isinstance(grid, tuple):
+            grid = list(grid)
+        for i, row in enumerate(grid):
+            if isinstance(row, tuple):
+                row = list(row)
+                grid[i] = row
+            for i, column in enumerate(row):
+                if isinstance(column, tuple):
+                    column = list(column)
+                    row[i] = column
+                # Data is angle, angle, energy, de/dang1, de/dang2, d^2e/dang1dang2
+                if unit.is_quantity(column[0]):
+                    column[0] = column[0].value_in_unit(unit.degree)
+                if unit.is_quantity(column[1]):
+                    column[1] = column[1].value_in_unit(unit.degree)
+                if unit.is_quantity(column[2]):
+                    column[2] = column[2].value_in_unit(unit.kilojoule_per_mole)
+                if len(column) > 3 and unit.is_quantity(column[3]):
+                    column[3] = column[3].value_in_unit(unit.kilojoule_per_mole/unit.radians)
+                if len(column) > 4 and unit.is_quantity(column[4]):
+                    column[4] = column[4].value_in_unit(unit.kilojoule_per_mole/unit.radians)
+                if len(column) > 5 and unit.is_quantity(column[5]):
+                    column[5] = column[5].value_in_unit(unit.kilojoule_per_mole/unit.radians**2)
+        return grid
+    try:
+        grid = copy.deepcopy(args[1])
+        if isinstance(args, tuple):
+            args = list(args)
+    except (NameError, UnboundLocalError):
+        try:
+            # Support numpy arrays
+            grid = grid.tolist()
+        except AttributeError:
+            grid = copy.deepcopy(grid)
+        grid = deunitize_grid(grid)
+    else:
+        args[1] = deunitize_grid(grid)
+%}