Add WaterCluster testsystem

openmmtools/tests/test_utils.py

@@ -79,7 +79,8 @@ def test_is_quantity_close():
                   (1.01325*unit.bar, 1.01325000006*unit.bar, True),
                   (1.01325*unit.bar, 1.0132500006*unit.bar, False)]
     for quantity1, quantity2, test_result in test_cases:
-        assert is_quantity_close(quantity1, quantity2) is test_result
+        msg = "Test failed: ({}, {}, {})".format(quantity1, quantity2, test_result)
+        assert is_quantity_close(quantity1, quantity2) == test_result, msg
 
     # Passing quantities with different units raise an exception.
     with nose.tools.assert_raises(TypeError):

openmmtools/testsystems.py

@@ -1721,6 +1721,87 @@ def __init__(self, nx=3, ny=3, nz=3, K=1.0 * unit.kilojoules_per_mole / unit.nan
 
         self.system, self.positions = system, positions
 
+
+class WaterCluster(TestSystem):
+    """Create a few water molecules in a harmonic restraining potential"""
+
+    def __init__(self,
+                 n_waters=20,
+                 K=1.0 * unit.kilojoules_per_mole / unit.nanometer ** 2,
+                 model='tip3p',
+                 constrained=True,
+                 **kwargs):
+        """
+
+        Parameters
+        ----------
+        n_waters : int
+            Number of water molecules in the cluster
+        K : simtk.unit.Quantity (energy / distance^2)
+            spring constant for restraining potential
+        model : string
+            Must be one of ['tip3p', 'tip4pew', 'tip5p', 'spce']
+        constrained: bool
+            Whether to use rigid water or not
+
+        Examples
+        --------
+
+        Create water cluster with default settings
+
+        >>> cluster = WaterCluster()
+        >>> system, positions = cluster.system, cluster.positions
+        """
+
+        TestSystem.__init__(self, **kwargs)
+
+        supported_models = ['tip3p', 'tip4pew', 'tip5p', 'spce']
+        if model not in supported_models:
+            raise Exception(
+                "Specified water model '%s' is not in list of supported models: %s" % (model, str(supported_models)))
+
+        # Load forcefield for solvent model and ions.
+        ff = app.ForceField(model + '.xml')
+
+        # Create empty topology and coordinates.
+        top = app.Topology()
+        pos = unit.Quantity((), unit.angstroms)
+
+        # Create new Modeller instance.
+        modeller = app.Modeller(top, pos)
+
+        # Add solvent
+        modeller.addSolvent(ff, model=model, numAdded=n_waters)
+
+        # Get new topology and coordinates.
+        new_top = modeller.getTopology()
+        new_pos = modeller.getPositions()
+
+        # Convert positions to numpy.
+        positions = unit.Quantity(numpy.array(new_pos / new_pos.unit), new_pos.unit)
+
+        # Create OpenMM System.
+        system = ff.createSystem(new_top,
+                                 nonbondedCutoff=openmm.NonbondedForce.NoCutoff,
+                                 constraints=None,
+                                 rigidWater=constrained,
+                                 removeCMMotion=False)
+
+        n_atoms = system.getNumParticles()
+        self.ndof = 3 * n_atoms - 3 * constrained
+
+        # Add a restraining potential centered at the origin.
+        energy_expression = '(K/2.0) * (x^2 + y^2 + z^2);'
+        energy_expression += 'K = %f;' % (K / (unit.kilojoules_per_mole / unit.nanometers ** 2))  # in OpenMM units
+        force = openmm.CustomExternalForce(energy_expression)
+        for particle_index in range(n_atoms):
+            force.addParticle(particle_index, [])
+        system.addForce(force)
+
+        self.topology = modeller.getTopology()
+        self.system = system
+        self.positions = positions
+
 #=============================================================================================
 # Lennard-Jones fluid
 #=============================================================================================

openmmtools/utils.py

@@ -363,7 +363,7 @@ def is_quantity_close(quantity1, quantity2):
     value2 = quantity2.value_in_unit_system(unit.md_unit_system)
 
     # np.isclose is not symmetric, so we make it so.
-    if value2 >= value1:
+    if abs(value2) >= abs(value1):
         return np.isclose(value1, value2, rtol=1e-10, atol=0.0)
     else:
         return np.isclose(value2, value1, rtol=1e-10, atol=0.0)