Merge pull request #322 from choderalab/water-cluster

Add WaterCluster testsystem

docs/releasehistory.rst

@@ -1,6 +1,13 @@
 Release History
 ===============
 
+Development snapshot
+====================
+New features
+------------
+- Add a ``WaterCluster`` testsystem (`#322 <https://github.com/choderalab/openmmtools/pull/322>`_)
+
+
 0.13.4 - Barostat/External Force Bugfix, Restart Robustness
 ===========================================================
 

docs/testsystems.rst

@@ -59,6 +59,7 @@ Clusters and simple fluids
     TolueneImplicitGBn
     TolueneImplicitGBn2
     MethanolBox
+    WaterCluster
 
 Solids
 ------
@@ -91,6 +92,7 @@ Water boxes
     GiantFlexibleDischargedWaterBox
     DischargedWaterBoxHsites
     AlchemicalWaterBox
+    WaterCluster
 
 Peptide and protein systems
 ---------------------------

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

@@ -331,6 +331,17 @@ def generate_dummy_trajectory(xyz, box):
 
     return traj
 
+def construct_restraining_potential(n_particles, K):
+    """Make a CustomExternalForce that puts a spring on particles 0 through n_particles"""
+
+    # 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_particles):
+        force.addParticle(particle_index, [])
+    return force
+
 
 #=============================================================================================
 # Thermodynamic state description
@@ -1711,16 +1722,87 @@ def __init__(self, nx=3, ny=3, nz=3, K=1.0 * unit.kilojoules_per_mole / unit.nan
             topology.addAtom('Ar', element, residue)
         self.topology = topology
 
-        # Add a restrining 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(natoms):
-            force.addParticle(particle_index, [])
-        system.addForce(force)
+        # Add a restraining potential centered at the origin.
+        system.addForce(construct_restraining_potential(natoms, K))
 
         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.
+        system.addForce(construct_restraining_potential(n_atoms, K))
+
+        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)