Make Numpy optional in CustomGBForces

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

@@ -40,6 +40,13 @@ from simtk.openmm import CustomGBForce, Discrete2DFunction
 import simtk.unit as u
 from math import floor
 
+_have_numpy = False
+try:
+    import numpy
+    _have_numpy = True
+except ImportError:
+    pass
+
 
 def strip_unit(value, unit):
     """Strip off any units and return value in unit"""
@@ -49,8 +56,7 @@ def strip_unit(value, unit):
 
 
 # GBn and GBn2 lookup tables
-d0= numpy.array(
-   [2.26685,2.32548,2.38397,2.44235,2.50057,2.55867,2.61663,2.67444,
+d0= [2.26685,2.32548,2.38397,2.44235,2.50057,2.55867,2.61663,2.67444,
     2.73212,2.78965,2.84705,2.9043,2.96141,3.0184,3.07524,3.13196,
     3.18854,3.24498,3.30132,3.35752,3.4136,
     2.31191,2.37017,2.4283,2.48632,2.5442,2.60197,2.65961,2.71711,
@@ -112,11 +118,9 @@ d0= numpy.array(
     4.04598,4.09974,4.15347,4.20715,4.26078,
     3.22855,3.28307,3.33751,3.39188,3.4462,3.50046,3.55466,3.6088,
     3.6629,3.71694,3.77095,3.82489,3.8788,3.93265,3.98646,4.04022,
-    4.09395,4.14762,4.20126,4.25485,4.3084],
-    numpy.double)
+    4.09395,4.14762,4.20126,4.25485,4.3084]
 
-m0=numpy.array(
-   [0.0381511,0.0338587,0.0301776,0.027003,0.0242506,0.0218529,
+m0= [0.0381511,0.0338587,0.0301776,0.027003,0.0242506,0.0218529,
     0.0197547,0.0179109,0.0162844,0.0148442,0.0135647,0.0124243,
     0.0114047,0.0104906,0.00966876,0.008928,0.0082587,0.00765255,
     0.00710237,0.00660196,0.00614589,
@@ -199,11 +203,10 @@ m0=numpy.array(
     0.0609642,0.0546169,0.0491183,0.0443295,0.0401388,0.036455,
     0.0332033,0.030322,0.0277596,0.0254732,0.0234266,0.0215892,
     0.0199351,0.018442,0.0170909,0.0158654,0.0147514,0.0137365,
-    0.0128101,0.0119627,0.0111863],
-    numpy.double)
+    0.0128101,0.0119627,0.0111863]
 # Rescale to nanometers
-d0/=10
-m0*=10
+d0 = [d/10.0 for d in d0]
+m0 = [m*10.0 for m in m0]
 
 
 def _get_bonded_atom_list(topology):
@@ -250,7 +253,11 @@ def _bondi_radii(topology):
         E.sulfur:       1.8,
         E.chlorine:     1.5,
     }
+
+    if _have_numpy:
         radii = numpy.empty(topology.getNumAtoms(), numpy.double)
+    else:
+        radii = [0]*topology.getNumAtoms()
     for i, atom in enumerate(topology.atoms()):
         radii[i] = element_to_radius.get(atom.element, default_radius)
     return radii  # converted to nanometers above
@@ -268,7 +275,10 @@ def _mbondi_radii(topology, all_bonds = None):
         E.sulfur:       1.8,
         E.chlorine:     1.5,
     }
+    if _have_numpy:
         radii = numpy.empty(topology.getNumAtoms(), numpy.double)
+    else:
+        radii = [0]*topology.getNumAtoms()
     if all_bonds is None:
         all_bonds = _get_bonded_atom_list(topology)
     for i, atom in enumerate(topology.atoms()):
@@ -303,7 +313,10 @@ def _mbondi2_radii(topology, all_bonds = None):
         E.sulfur:       1.8,
         E.chlorine:     1.5,
     }
+    if _have_numpy:
         radii = numpy.empty(topology.getNumAtoms(), numpy.double)
+    else:
+        radii = [0]*topology.getNumAtoms()
     if all_bonds is None:
         all_bonds = _get_bonded_atom_list(topology)
     for i, atom in enumerate(topology.atoms()):
@@ -387,8 +400,6 @@ _SCREEN_PARAMETERS = { # normal, GBn, GBn2
         E.sulfur : (0.96, 0.602256336067, -0.703469),
         None : (0.8, 0.5, 0.5) # default
 }
-for key, val in _SCREEN_PARAMETERS.items():
-    _SCREEN_PARAMETERS[key] = numpy.array(val, numpy.double)
 _SCREEN_PARAMETERS[E.deuterium] = _SCREEN_PARAMETERS[E.hydrogen]
 
 def _screen_parameter(atom):
@@ -453,14 +464,20 @@ class CustomAmberGBForceBase(CustomGBForce):
         -------
         None
         """
-        if isinstance(params, list):
+        if isinstance(params, list) and _have_numpy:
             params = numpy.array(params)
         else:
             params = copy.deepcopy(params)
 
+        if _have_numpy:
             params[:,self.RADIUS_ARG_POSITION] -= self.OFFSET
             params[:,self.SCREEN_POSITION] *= params[:,self.RADIUS_ARG_POSITION]
             self.parameters.extend(params.tolist())
+        else:
+            for p in params:
+                p[self.RADIUS_ARG_POSITION] -= self.OFFSET
+                p[self.SCREEN_POSITION] *= p[self.RADIUS_ARG_POSITION]
+            self.parameters.extend(params)
 
     @staticmethod
     def getStandardParameters(topology):
@@ -551,8 +568,11 @@ class GBSAHCTForce(CustomAmberGBForceBase):
             at the beginning of the list
 
         """
+        if _have_numpy:
             radii = numpy.empty((topology.getNumAtoms(), 2), numpy.double)
             radii[:,0] = _mbondi_radii(topology)/10
+        else:
+            radii = [[r/10, 0] for r in _mbondi_radii(topology)]
         for rad, atom in zip(radii, topology.atoms()):
             rad[1] = _screen_parameter(atom)[0]
         return radii
@@ -614,8 +634,11 @@ class GBSAOBC1Force(CustomAmberGBForceBase):
             at the beginning of the list
 
         """
+        if _have_numpy:
             radii = numpy.empty((topology.getNumAtoms(), 2), numpy.double)
             radii[:,0] = _mbondi2_radii(topology)/10
+        else:
+            radii = [[r/10, 0] for r in _mbondi2_radii(topology)]
         for rad, atom in zip(radii, topology.atoms()):
             rad[1] = _screen_parameter(atom)[0]
         return radii
@@ -704,10 +727,17 @@ class GBSAGBnForce(CustomAmberGBForceBase):
             raise ValueError('Radii must be between 1 and 2 Angstroms for neck lookup')
 
     def addParticles(self, parameters):
+        bad_radius = False
+        if _have_numpy:
             if isinstance(parameters, list):
                 parameters = numpy.array(parameters)
             radii = parameters[:,1]
             if numpy.any(numpy.logical_or(radii<0.1, radii>0.2)):
+                bad_radius = True
+        else:
+            if any(p[1]<0.1 or p[1]>0.2 for p in parameters):
+                bad_radius = True
+        if bad_radius:
             raise ValueError('Radii must be between 1 and 2 Angstroms for neck lookup')
         CustomAmberGBForceBase.addParticles(self, parameters)
 
@@ -734,8 +764,11 @@ class GBSAGBnForce(CustomAmberGBForceBase):
             at the beginning of the list
 
         """
+        if _have_numpy:
             radii = numpy.empty((topology.getNumAtoms(), 2), numpy.double)
             radii[:,0] = _bondi_radii(topology)/10
+        else:
+            radii = [[r/10, 0] for r in _bondi_radii(topology)]
         for rad, atom in zip(radii, topology.atoms()):
             rad[1] = _screen_parameter(atom)[1]
         return radii
@@ -834,14 +867,14 @@ class GBSAGBn2Force(GBSAGBnForce):
         GBSAGBnForce.__init__(self, solventDielectric, soluteDielectric, SA, cutoff, kappa)
 
     _atom_params = {
-        E.hydrogen:     numpy.array([0.788440, 0.798699, 0.437334]),
-        E.deuterium:    numpy.array([0.788440, 0.798699, 0.437334]),
-        E.carbon:       numpy.array([0.733756, 0.506378, 0.205844]),
-        E.nitrogen:     numpy.array([0.503364, 0.316828, 0.192915]),
-        E.oxygen:       numpy.array([0.867814, 0.876635, 0.387882]),
-        E.sulfur:       numpy.array([0.867814, 0.876635, 0.387882]),
+        E.hydrogen:     [0.788440, 0.798699, 0.437334],
+        E.deuterium:    [0.788440, 0.798699, 0.437334],
+        E.carbon:       [0.733756, 0.506378, 0.205844],
+        E.nitrogen:     [0.503364, 0.316828, 0.192915],
+        E.oxygen:       [0.867814, 0.876635, 0.387882],
+        E.sulfur:       [0.867814, 0.876635, 0.387882],
     }
-    _default_atom_params = numpy.array([0.8, 4.85, 0.5]);
+    _default_atom_params = [0.8, 4.85, 0.5]
 
     @classmethod
     def getStandardParameters(cls, topology):
@@ -861,11 +894,18 @@ class GBSAGBn2Force(GBSAGBnForce):
 
         """
         natoms = topology.getNumAtoms()
+        if _have_numpy:
             radii = numpy.empty([natoms,5], numpy.double)
             radii[:,0] = _mbondi3_radii(topology)/10
             for atom, rad in zip(topology.atoms(), radii):
                 rad[1] = _screen_parameter(atom)[2]
                 rad[2:] = cls._atom_params.get(atom.element, cls._default_atom_params)
+        else:
+            radii = [[r/10, 0, 0, 0, 0] for r in _mbondi3_radii(topology)]
+            for atom, rad in zip(topology.atoms(), radii):
+                rad[1] = _screen_parameter(atom)[2]
+                for i, p in enumerate(cls._atom_params.get(atom.element, cls._default_atom_params)):
+                    rad[2+i] = p
         return radii
 
     def _addEnergyTerms(self):