Repository: choderalab/openmmtools

Issue:
Where to put replace_reaction_field?
Currently it's an instance method of `AbsoluteAlchemicalFactory`. The relative factory will surely want to use it too, so we should probably move it as a function at the `alchemy` module level at least. Are there any other plans for that feature?

Discussion:

Patch:
diff --git a/docs/forcefactories.rst b/docs/forcefactories.rst
new file mode 100644
index 0000000..cd24f18
--- /dev/null
+++ b/docs/forcefactories.rst
@@ -0,0 +1,18 @@
+.. _forcefactories::
+
+Cache
+=====
+
+The module :mod:`openmmtools.forcefactories` implements utility methods and factories to configure system forces.
+
+ - :func:`replace_reaction_field`: Configure a system to model the electrostatics with an :class:`UnshiftedReactionField` force.
+
+Cache objects
+-------------
+
+.. currentmodule:: openmmtools.forcefactories
+.. autosummary::
+    :nosignatures:
+    :toctree: api/generated/
+
+    replace_reaction_field
diff --git a/docs/forces.rst b/docs/forces.rst
new file mode 100644
index 0000000..49d9bcb
--- /dev/null
+++ b/docs/forces.rst
@@ -0,0 +1,22 @@
+.. _forces::
+
+Forces
+======
+
+The module :mod:`openmmtools.forces` implements custom forces that are not natively found in OpenMM.
+
+ - :class:`UnshiftedReactionFieldForce`: A `CustomNonbondedForce <http://docs.openmm.org/7.1.0/api-python/generated/simtk.openmm.openmm.CustomNonbondedForce.html>`_ implementing a reaction field variant with `c_rf` term set to zero and a switching function. Using the native OpenMM reaction field implementation with `c_rf != 0` can cause issues with hydration free energy calculations.
+ - :func:`find_nonbonded_force`: Find the first ``NonbondedForce`` in an OpenMM ``System``.
+ - :func:`iterate_nonbonded_forces`: Iterate over all the ``NonbondedForce``s in an OpenMM ``System``.
+
+Cache objects
+-------------
+
+.. currentmodule:: openmmtools.forces
+.. autosummary::
+    :nosignatures:
+    :toctree: api/generated/
+
+    UnshiftedReactionFieldForce
+    find_nonbonded_force
+    iterate_nonbonded_forces
diff --git a/docs/index.rst b/docs/index.rst
index 108ac4c..c561815 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -53,6 +53,8 @@ Modules
   mcmc
   sampling
   alchemy
+  forces
+  forcefactories
   utils
   scripts
 
diff --git a/openmmtools/alchemy.py b/openmmtools/alchemy.py
index 972ef8e..81e7650 100644
--- a/openmmtools/alchemy.py
+++ b/openmmtools/alchemy.py
@@ -36,14 +36,13 @@ usable for the calculation of free energy differences of hydration or ligand bin
 # =============================================================================
 
 import copy
-import inspect
 import logging
 import collections
 
 import numpy as np
 from simtk import openmm, unit
 
-from openmmtools import states, utils
+from openmmtools import states, forcefactories, utils
 from openmmtools.constants import ONE_4PI_EPS0
 
 logger = logging.getLogger(__name__)
@@ -843,9 +842,10 @@ class AbsoluteAlchemicalFactory(object):
 
         # If the System uses a NonbondedForce, replace its NonbondedForce implementation of reaction field
         # with a Custom*Force implementation that uses c_rf = 0.
-        # NOTE: This adds an additional CustomNonbondedForce and CustomBondForce
-        if (self.alchemical_rf_treatment == 'switched'):
-            alchemical_system = self.replace_reaction_field(alchemical_system)
+        # NOTE: This adds an additional CustomNonbondedForce
+        if self.alchemical_rf_treatment == 'switched':
+            forcefactories.replace_reaction_field(alchemical_system, return_copy=False,
+                                                  switch_width=self.switch_width)
 
         return alchemical_system
 
@@ -905,73 +905,6 @@ class AbsoluteAlchemicalFactory(object):
 
         return energy_components
 
-    def replace_reaction_field(self, reference_system):
-        """Replace reaction-field electrostatics with Custom*Force terms to ensure c_rf = 0.
-
-        .. warning:: Unstable API.
-            This method is still experimental. It could be moved to some
-            other module or have its signature changed in the near future.
-
-        A deep copy of the system is made.
-
-        If reaction field electrostatics is in use, this will add a CustomNonbondedForce and CustomBondForce to the System
-        for each NonbondedForce that utilizes CutoffPeriodic.
-
-        Note that the resulting System object can NOT be fed to `create_alchemical_system` since the CustomNonbondedForce
-        will not be recognized and re-coded.
-
-        Parameters
-        ----------
-        reference_system : simtk.openmm.System
-            The system to use as a reference for the creation of the
-            alchemical system. This will not be modified.
-
-        Returns
-        -------
-        system : simtk.openmm.System
-            System with reaction-field converted to c_rf = 0
-
-        """
-        system = copy.deepcopy(reference_system)
-        for force_index, reference_force in enumerate(system.getForces()):
-            reference_force_name = reference_force.__class__.__name__
-            if (reference_force_name == 'NonbondedForce' and
-                        reference_force.getNonbondedMethod() == openmm.NonbondedForce.CutoffPeriodic):
-                # Create CustomNonbondedForce to handle switched reaction field
-                epsilon_solvent = reference_force.getReactionFieldDielectric()
-                r_cutoff = reference_force.getCutoffDistance()
-                energy_expression = "ONE_4PI_EPS0*chargeprod*(r^(-1) + k_rf*r^2);"  # Omit c_rf constant term.
-                k_rf = r_cutoff**(-3) * ((epsilon_solvent - 1.0) / (2.0*epsilon_solvent + 1.0))
-                energy_expression += "chargeprod = charge1*charge2;"
-                energy_expression += "k_rf = %f;" % (k_rf.value_in_unit_system(unit.md_unit_system))
-                energy_expression += "ONE_4PI_EPS0 = %f;" % ONE_4PI_EPS0 # already in OpenMM units
-                custom_nonbonded_force = openmm.CustomNonbondedForce(energy_expression)
-                custom_nonbonded_force.addPerParticleParameter("charge")
-                custom_nonbonded_force.setNonbondedMethod(openmm.CustomNonbondedForce.CutoffPeriodic)
-                custom_nonbonded_force.setCutoffDistance(reference_force.getCutoffDistance())
-                custom_nonbonded_force.setUseLongRangeCorrection(False)
-                system.addForce(custom_nonbonded_force)
-
-                # Add switch
-                if self.switch_width is not None:
-                    custom_nonbonded_force.setUseSwitchingFunction(True)
-                    custom_nonbonded_force.setSwitchingDistance(reference_force.getCutoffDistance() - self.switch_width)
-                else:
-                    custom_nonbonded_force.setUseSwitchingFunction(False)
-
-                # Rewrite particle charges
-                for particle_index in range(reference_force.getNumParticles()):
-                    [charge, sigma, epsilon] = reference_force.getParticleParameters(particle_index)
-                    reference_force.setParticleParameters(particle_index, abs(0.0*charge), sigma, epsilon)
-                    custom_nonbonded_force.addParticle([charge])
-
-                # Add exclusions to CustomNonbondedForce.
-                for exception_index in range(reference_force.getNumExceptions()):
-                    iatom, jatom, chargeprod, sigma, epsilon = reference_force.getExceptionParameters(exception_index)
-                    custom_nonbonded_force.addExclusion(iatom, jatom)
-
-        return system
-
     # -------------------------------------------------------------------------
     # Internal usage: AlchemicalRegion
     # -------------------------------------------------------------------------
@@ -2091,20 +2024,20 @@ class AbsoluteAlchemicalFactory(object):
                     add_label(label.format('non-'), force_index2)
 
             # If they are both empty they are identical and any label works.
-            elif force1.getNumBonds() == 0:
+            elif force1.getNumBonds() == 0 and force2.getNumBonds() == 0:
                 add_label(label.format(''), force_index1)
                 add_label(label.format('non-'), force_index2)
 
             # We check that the bond atoms are both alchemical or not.
             else:
-                atom_i, atom_j, _ = force1.getBondParameters(0)
+                atom_i, atom_j, _ = force2.getBondParameters(0)
                 both_alchemical = atom_i in alchemical_atoms and atom_j in alchemical_atoms
                 if both_alchemical:
-                    add_label(label.format(''), force_index1)
-                    add_label(label.format('non-'), force_index2)
-                else:
-                    add_label(label.format('non-'), force_index1)
                     add_label(label.format(''), force_index2)
+                    add_label(label.format('non-'), force_index1)
+                else:
+                    add_label(label.format('non-'), force_index2)
+                    add_label(label.format(''), force_index1)
 
         return force_labels
 
diff --git a/openmmtools/cache.py b/openmmtools/cache.py
index 11ac9b5..87ec8d6 100644
--- a/openmmtools/cache.py
+++ b/openmmtools/cache.py
@@ -347,8 +347,7 @@ class ContextCache(object):
 
         If the integrator is not provided, this will search the cache for
         any Context in the given ThermodynamicState, regardless of its
-        integrator. In this case, the method guarantees that two consecutive
-        calls with the same thermodynamic state will retrieve the same context.
+        integrator.
 
         This creates a new Context if no compatible one has been cached.
         If a compatible Context exists, the ThermodynamicState is applied
@@ -388,10 +387,8 @@ class ContextCache(object):
                 # Only one match.
                 context = self._lru[matching_context_ids[0]]
             else:
-                # Multiple matches, prefer the non-default Integrator.
-                # Always pick the least recently used to make two consective
-                # calls retrieving the same integrator.
-                for context_id in reversed(matching_context_ids):
+                # Multiple matches, prefer non-default Integrator.
+                for context_id in matching_context_ids:
                     if context_id[1] != self._default_integrator_id():
                         context = self._lru[context_id]
                         break
diff --git a/openmmtools/forcefactories.py b/openmmtools/forcefactories.py
new file mode 100644
index 0000000..777d110
--- /dev/null
+++ b/openmmtools/forcefactories.py
@@ -0,0 +1,79 @@
+#!/usr/bin/env python
+
+# =============================================================================
+# MODULE DOCSTRING
+# =============================================================================
+
+"""
+Factories to manipulate OpenMM System forces.
+
+"""
+
+
+# =============================================================================
+# GLOBAL IMPORTS
+# =============================================================================
+
+import copy
+
+from simtk import openmm, unit
+
+from openmmtools import forces
+
+
+# =============================================================================
+# FACTORY FUNCTIONS
+# =============================================================================
+
+def replace_reaction_field(reference_system, switch_width=1.0*unit.angstrom,
+                           return_copy=True):
+    """Return a system converted to use a switched reaction-field electrostatics.
+
+    This will add an `UnshiftedReactionFieldForce` for each `NonbondedForce`
+    that utilizes `CutoffPeriodic`.
+
+    Note that `AbsoluteAlchemicalFactory.create_alchemical_system()` can NOT
+    handle the resulting `System` object yet since the `CustomNonbondedForce`
+    are not recognized and re-coded.
+
+    Parameters
+    ----------
+    reference_system : simtk.openmm.System
+        The system to use as a reference. This will be modified only if
+        `return_copy` is `False`.
+    switch_width : simtk.unit.Quantity, default 1.0*angstrom
+        Switch width for electrostatics (units of distance).
+    return_copy : bool
+        If `True`, `reference_system` is not modified, and a copy is returned.
+        Setting it to `False` speeds up the function execution but modifies
+        the `reference_system` object.
+
+    Returns
+    -------
+    system : simtk.openmm.System
+        System with reaction-field converted to c_rf = 0
+
+    """
+    if return_copy:
+        system = copy.deepcopy(reference_system)
+    else:
+        system = reference_system
+
+    # Add an UnshiftedReactionFieldForce for each CutoffPeriodic NonbondedForce.
+    for reference_force in forces.iterate_nonbonded_forces(system):
+        if reference_force.getNonbondedMethod() == openmm.NonbondedForce.CutoffPeriodic:
+            reaction_field_force = forces.UnshiftedReactionFieldForce.from_nonbonded_force(reference_force,
+                                                                                           switch_width=switch_width)
+            system.addForce(reaction_field_force)
+
+            # Remove particle electrostatics from reference force, but leave exceptions.
+            for particle_index in range(reference_force.getNumParticles()):
+                charge, sigma, epsilon = reference_force.getParticleParameters(particle_index)
+                reference_force.setParticleParameters(particle_index, abs(0.0*charge), sigma, epsilon)
+
+    return system
+
+
+if __name__ == '__main__':
+    import doctest
+    doctest.testmod()
diff --git a/openmmtools/forces.py b/openmmtools/forces.py
new file mode 100644
index 0000000..503301f
--- /dev/null
+++ b/openmmtools/forces.py
@@ -0,0 +1,205 @@
+#!/usr/bin/env python
+
+# =============================================================================
+# MODULE DOCSTRING
+# =============================================================================
+
+"""
+Custom OpenMM Forces classes and utilities.
+
+"""
+
+
+# =============================================================================
+# GLOBAL IMPORTS
+# =============================================================================
+
+from simtk import openmm, unit
+
+from openmmtools.constants import ONE_4PI_EPS0
+
+
+# =============================================================================
+# UTILITY FUNCTIONS
+# =============================================================================
+
+def find_nonbonded_force(system):
+    """Find the first OpenMM `NonbondedForce` in the system.
+
+    Parameters
+    ----------
+    system : simtk.openmm.System
+        The system to search.
+
+    Returns
+    -------
+    nonbonded_force : simtk.openmm.NonbondedForce
+        The first `NonbondedForce` object in `system`.
+
+    Raises
+    ------
+    ValueError
+        If the system contains multiple `NonbondedForce`s
+
+    """
+    nonbonded_force = None
+    for force in system.getForces():
+        if isinstance(force, openmm.NonbondedForce):
+            if nonbonded_force is not None:
+                raise ValueError('The System has multiple NonbondedForces')
+            nonbonded_force = force
+    return nonbonded_force
+
+
+def iterate_nonbonded_forces(system):
+    """Iterate over all OpenMM `NonbondedForce`s in `system`.
+
+    Parameters
+    ----------
+    system : simtk.openmm.System
+        The system to search.
+
+    Yields
+    ------
+    force : simtk.openmm.NonbondedForce
+        A `NonbondedForce` object in `system`.
+
+    """
+    for force in system.getForces():
+        if isinstance(force, openmm.NonbondedForce):
+            yield force
+
+
+# =============================================================================
+# REACTION FIELD
+# =============================================================================
+
+class UnshiftedReactionFieldForce(openmm.CustomNonbondedForce):
+    """A force modelling switched reaction-field electrostatics.
+
+    Contrarily to a normal `NonbondedForce` with `CutoffPeriodic` nonbonded
+    method, this force sets the `c_rf` to 0.0 and uses a switching function
+    to avoid forces discontinuities at the cutoff distance.
+
+    Parameters
+    ----------
+    cutoff_distance : simtk.unit.Quantity, default 15*angstroms
+        The cutoff distance (units of distance).
+    switch_width : simtk.unit.Quantity, default 1.0*angstrom
+        Switch width for electrostatics (units of distance).
+    reaction_field_dielectric : float
+        The dielectric constant used for the solvent.
+
+    """
+
+    def __init__(self, cutoff_distance=15*unit.angstroms, switch_width=1.0*unit.angstrom,
+                 reaction_field_dielectric=78.3):
+        k_rf = cutoff_distance**(-3) * (reaction_field_dielectric - 1.0) / (2.0*reaction_field_dielectric + 1.0)
+
+        # Energy expression omits c_rf constant term.
+        energy_expression = "ONE_4PI_EPS0*chargeprod*(r^(-1) + k_rf*r^2);"
+        energy_expression += "chargeprod = charge1*charge2;"
+        energy_expression += "k_rf = {:f};".format(k_rf.value_in_unit_system(unit.md_unit_system))
+        energy_expression += "ONE_4PI_EPS0 = {:f};".format(ONE_4PI_EPS0)  # already in OpenMM units
+
+        # Create CustomNonbondedForce.
+        super(UnshiftedReactionFieldForce, self).__init__(energy_expression)
+
+        # Add parameters.
+        self.addPerParticleParameter("charge")
+
+        # Configure force.
+        self.setNonbondedMethod(openmm.CustomNonbondedForce.CutoffPeriodic)
+        self.setCutoffDistance(cutoff_distance)
+        self.setUseLongRangeCorrection(False)
+        if switch_width is not None:
+            self.setUseSwitchingFunction(True)
+            self.setSwitchingDistance(cutoff_distance - switch_width)
+        else:  # Truncated
+            self.setUseSwitchingFunction(False)
+
+    @classmethod
+    def from_nonbonded_force(cls, nonbonded_force, switch_width=1.0*unit.angstrom):
+        """Copy constructor from an OpenMM `NonbondedForce`.
+
+        The returned force has same cutoff distance and dielectric, and
+        its particles have the same charges. Exclusions corresponding to
+        `nonbonded_force` exceptions are also added.
+
+        .. warning
+            This only creates the force object. The electrostatics in
+            `nonbonded_force` remains unmodified. Use the function
+            `replace_reaction_field` to correctly convert a system to
+            use an unshifted reaction field potential.
+
+        Parameters
+        ----------
+        nonbonded_force : simtk.openmm.NonbondedForce
+            The nonbonded force to copy.
+        switch_width : simtk.unit.Quantity
+            Switch width for electrostatics (units of distance).
+
+        Returns
+        -------
+        reaction_field_force : UnshiftedReactionFieldForce
+            The reaction field force with copied particles.
+
+        """
+        # OpenMM gives unitless values.
+        cutoff_distance = nonbonded_force.getCutoffDistance()
+        reaction_field_dielectric = nonbonded_force.getReactionFieldDielectric()
+        reaction_field_force = cls(cutoff_distance, switch_width, reaction_field_dielectric)
+
+        # Set particle charges.
+        for particle_index in range(nonbonded_force.getNumParticles()):
+            charge, sigma, epsilon = nonbonded_force.getParticleParameters(particle_index)
+            reaction_field_force.addParticle([charge])
+
+        # Add exclusions to CustomNonbondedForce.
+        for exception_index in range(nonbonded_force.getNumExceptions()):
+            iatom, jatom, chargeprod, sigma, epsilon = nonbonded_force.getExceptionParameters(exception_index)
+            reaction_field_force.addExclusion(iatom, jatom)
+
+        return reaction_field_force
+
+    @classmethod
+    def from_system(cls, system, switch_width=1.0*unit.angstrom):
+        """Copy constructor from the first OpenMM `NonbondedForce` in `system`.
+
+        If multiple `NonbondedForce`s are found, an exception is raised.
+
+        .. warning
+            This only creates the force object. The electrostatics in
+            `nonbonded_force` remains unmodified. Use the function
+            `replace_reaction_field` to correctly convert a system to
+            use an unshifted reaction field potential.
+
+        Parameters
+        ----------
+        system : simtk.openmm.System
+            The system containing the nonbonded force to copy.
+        switch_width : simtk.unit.Quantity
+            Switch width for electrostatics (units of distance).
+
+        Returns
+        -------
+        reaction_field_force : UnshiftedReactionFieldForce
+            The reaction field force.
+
+        Raises
+        ------
+        ValueError
+            If multiple `NonbondedForce`s are found in the system.
+
+        See Also
+        --------
+        UnshiftedReactionField.from_nonbonded_force
+
+        """
+        nonbonded_force = find_nonbonded_force(system)
+        return cls.from_nonbonded_force(nonbonded_force, switch_width)
+
+
+if __name__ == '__main__':
+    import doctest
+    doctest.testmod()


Test patch:
diff --git a/openmmtools/tests/test_alchemy.py b/openmmtools/tests/test_alchemy.py
index aac29c1..c6581c1 100644
--- a/openmmtools/tests/test_alchemy.py
+++ b/openmmtools/tests/test_alchemy.py
@@ -27,7 +27,7 @@ import nose
 import scipy
 from nose.plugins.attrib import attr
 
-from openmmtools import testsystems
+from openmmtools import testsystems, forces
 from openmmtools.constants import kB
 from openmmtools.alchemy import *
 
@@ -41,9 +41,7 @@ logger = logging.getLogger(__name__)
 temperature = 300.0 * unit.kelvin  # reference temperature
 # MAX_DELTA = 0.01 * kB * temperature # maximum allowable deviation
 MAX_DELTA = 1.0 * kB * temperature  # maximum allowable deviation
-MAX_FORCE_RELATIVE_ERROR = 1.0e-6 # maximum allowable relative force error
 GLOBAL_ENERGY_UNIT = unit.kilojoules_per_mole  # controls printed units
-GLOBAL_FORCE_UNIT = unit.kilojoules_per_mole / unit.nanometers # controls printed units
 GLOBAL_ALCHEMY_PLATFORM = None  # This is used in every energy calculation.
 # GLOBAL_ALCHEMY_PLATFORM = openmm.Platform.getPlatformByName('OpenCL') # DEBUG: Use OpenCL over CPU platform for testing since OpenCL is deterministic, while CPU is not
 
@@ -88,52 +86,6 @@ def compute_energy(system, positions, platform=None, force_group=-1):
     return potential
 
 
-def compute_forces(system, positions, platform=None, force_group=-1):
-    """Compute forces of the system in the given positions.
-
-    Parameters
-    ----------
-    platform : simtk.openmm.Platform or None, optional
-        If None, the global GLOBAL_ALCHEMY_PLATFORM will be used.
-    force_group : int flag or set of int, optional
-        Passed to the groups argument of Context.getState().
-
-    """
-    timestep = 1.0 * unit.femtoseconds
-    integrator = openmm.VerletIntegrator(timestep)
-    context = create_context(system, integrator, platform)
-    context.setPositions(positions)
-    state = context.getState(getForces=True, groups=force_group)
-    forces = state.getForces(asNumpy=True)
-    del context, integrator, state
-    return forces
-
-
-def generate_new_positions(system, positions, platform=None, nsteps=50):
-    """Generate new positions by taking a few steps from the old positions.
-    Parameters
-    ----------
-    platform : simtk.openmm.Platform or None, optional
-        If None, the global GLOBAL_ALCHEMY_PLATFORM will be used.
-    nsteps : int, optional, default=50
-        Number of steps of dynamics to take.
-    Returns
-    -------
-    new_positions : simtk.unit.Quantity of shape [nparticles,3] with units compatible with distance
-        New positions
-    """
-    temperature = 300 * unit.kelvin
-    collision_rate = 90 / unit.picoseconds
-    timestep = 1.0 * unit.femtoseconds
-    integrator = openmm.LangevinIntegrator(temperature, collision_rate, timestep)
-    context = create_context(system, integrator, platform)
-    context.setPositions(positions)
-    integrator.step(nsteps)
-    new_positions = context.getState(getPositions=True).getPositions(asNumpy=True)
-    del context, integrator
-    return new_positions
-
-
 def minimize(system, positions, platform=None, tolerance=1.0*unit.kilocalories_per_mole/unit.angstroms, maxIterations=50):
     """Minimize the energy of the given system.
 
@@ -162,7 +114,7 @@ def minimize(system, positions, platform=None, tolerance=1.0*unit.kilocalories_p
     return minimized_positions
 
 
-def compute_energy_force(system, positions, force_name):
+def compute_force_energy(system, positions, force_name):
     """Compute the energy of the force with the given name."""
     system = copy.deepcopy(system)  # Copy to avoid modifications
     force_name_index = 1
@@ -192,11 +144,11 @@ def assert_almost_equal(energy1, energy2, err_msg):
 
 
 def dissect_nonbonded_energy(reference_system, positions, alchemical_atoms):
-    """Dissect the contributions to NonbondedForce of the reference system by atom group
-    and sterics/electrostatics.
+    """Dissect the nonbonded energy contributions of the reference system
+    by atom group and sterics/electrostatics.
 
-    Note that this can only work on reference_system objects whose CutoffPeriodic forces
-    have not been replaced by Custom*Force objects to set c_rf = 0.
+    This works also for systems objects whose CutoffPeriodic force
+    has been replaced by a CustomNonbondedForce to set c_rf = 0.
 
     Parameters
     ----------
@@ -229,60 +181,76 @@ def dissect_nonbonded_energy(reference_system, positions, alchemical_atoms):
 
     """
 
-    def turn_off(force, sterics=False, electrostatics=False,
+    def turn_off(system, sterics=False, electrostatics=False,
                  exceptions=False, only_atoms=frozenset()):
+        """Turn off sterics and/or electrostatics interactions.
+
+        If `exceptions` is True, only the exceptions are turned off.
+        Support also system that have gone through replace_reaction_field.
+        The `system` must have only nonbonded forces.
+        If `only_atoms` is specified, only the those atoms will be turned off.
+
+        """
         if len(only_atoms) == 0:  # if empty, turn off all particles
-            only_atoms = set(range(force.getNumParticles()))
-        e_coeff = 0.0 if sterics else 1.0
-        c_coeff = 0.0 if electrostatics else 1.0
+            only_atoms = set(range(system.getNumParticles()))
+        epsilon_coeff = 0.0 if sterics else 1.0
+        charge_coeff = 0.0 if electrostatics else 1.0
+
+        # Only a Nonbonded and a CustomNonbonded (for RF) force should be here.
         if exceptions:  # Turn off exceptions
-            for exception_index in range(force.getNumExceptions()):
-                [iatom, jatom, charge, sigma, epsilon] = force.getExceptionParameters(exception_index)
+            nonbonded_force = system.getForces()[0]  # NonbondedForce
+            for exception_index in range(nonbonded_force.getNumExceptions()):
+                iatom, jatom, charge, sigma, epsilon = nonbonded_force.getExceptionParameters(exception_index)
                 if iatom in only_atoms or jatom in only_atoms:
-                    force.setExceptionParameters(exception_index, iatom, jatom, c_coeff*charge,
-                                                 sigma, e_coeff*epsilon)
+                    nonbonded_force.setExceptionParameters(exception_index, iatom, jatom,
+                                                           charge_coeff*charge, sigma, epsilon_coeff*epsilon)
         else:  # Turn off particle interactions
-            for particle_index in range(force.getNumParticles()):
-                if particle_index in only_atoms:
-                    [charge, sigma, epsilon] = force.getParticleParameters(particle_index)
-                    force.setParticleParameters(particle_index, c_coeff*charge, sigma, e_coeff*epsilon)
-
-    def restore_system(reference_system):
-        system = copy.deepcopy(reference_system)
-        nonbonded_force = system.getForces()[0]
-        return system, nonbonded_force
+            for force in system.getForces():
+                for particle_index in range(force.getNumParticles()):
+                    if particle_index in only_atoms:
+                        # Convert tuple parameters to list to allow changes.
+                        parameters = list(force.getParticleParameters(particle_index))
+                        parameters[0] *= charge_coeff  # charge
+                        try:  # CustomNonbondedForce
+                            force.setParticleParameters(particle_index, parameters)
+                        except TypeError:  # NonbondedForce
+                            parameters[2] *= epsilon_coeff  # epsilon
+                            force.setParticleParameters(particle_index, *parameters)
 
     nonalchemical_atoms = set(range(reference_system.getNumParticles())).difference(alchemical_atoms)
 
-    # Remove all forces but NonbondedForce and, if CutoffPeriodic is in use,
-    # the CustomNonbondedForce and CustomBondForce used to replace it
+    # Remove all forces but NonbondedForce and eventually the
+    # CustomNonbondedForce used to model reaction field.
     reference_system = copy.deepcopy(reference_system)  # don't modify original system
     forces_to_remove = list()
     for force_index, force in enumerate(reference_system.getForces()):
-        if force.__class__.__name__ != 'NonbondedForce':
-            forces_to_remove.append(force_index)
-        else:
-            force.setForceGroup(0)
+        force.setForceGroup(0)
+        if isinstance(force, openmm.NonbondedForce):
             force.setReciprocalSpaceForceGroup(30)  # separate PME reciprocal from direct space
+        # We keep only CustomNonbondedForces that are not alchemically modified.
+        elif not (isinstance(force, openmm.CustomNonbondedForce) and
+                          'lambda' not in force.getEnergyFunction()):
+            forces_to_remove.append(force_index)
+
     for force_index in reversed(forces_to_remove):
         reference_system.removeForce(force_index)
-    assert len(reference_system.getForces()) == 1
+    assert len(reference_system.getForces()) <= 2
 
     # Compute particle interactions between different groups of atoms
     # ----------------------------------------------------------------
-    system, nonbonded_force = restore_system(reference_system)
+    system = copy.deepcopy(reference_system)
 
     # Compute total energy from nonbonded interactions
     tot_energy = compute_energy(system, positions)
     tot_reciprocal_energy = compute_energy(system, positions, force_group={30})
 
     # Compute contributions from particle sterics
-    turn_off(nonbonded_force, sterics=True, only_atoms=alchemical_atoms)
+    turn_off(system, sterics=True, only_atoms=alchemical_atoms)
     tot_energy_no_alchem_particle_sterics = compute_energy(system, positions)
-    system, nonbonded_force = restore_system(reference_system)  # Restore alchemical sterics
-    turn_off(nonbonded_force, sterics=True, only_atoms=nonalchemical_atoms)
+    system = copy.deepcopy(reference_system)  # Restore alchemical sterics
+    turn_off(system, sterics=True, only_atoms=nonalchemical_atoms)
     tot_energy_no_nonalchem_particle_sterics = compute_energy(system, positions)
-    turn_off(nonbonded_force, sterics=True)
+    turn_off(system, sterics=True)
     tot_energy_no_particle_sterics = compute_energy(system, positions)
 
     tot_particle_sterics = tot_energy - tot_energy_no_particle_sterics
@@ -291,15 +259,15 @@ def dissect_nonbonded_energy(reference_system, positions, alchemical_atoms):
     na_particle_sterics = tot_particle_sterics - nn_particle_sterics - aa_particle_sterics
 
     # Compute contributions from particle electrostatics
-    system, nonbonded_force = restore_system(reference_system)  # Restore sterics
-    turn_off(nonbonded_force, electrostatics=True, only_atoms=alchemical_atoms)
+    system = copy.deepcopy(reference_system)  # Restore sterics
+    turn_off(system, electrostatics=True, only_atoms=alchemical_atoms)
     tot_energy_no_alchem_particle_electro = compute_energy(system, positions)
     nn_reciprocal_energy = compute_energy(system, positions, force_group={30})
-    system, nonbonded_force = restore_system(reference_system)  # Restore alchemical electrostatics
-    turn_off(nonbonded_force, electrostatics=True, only_atoms=nonalchemical_atoms)
+    system = copy.deepcopy(reference_system)  # Restore alchemical electrostatics
+    turn_off(system, electrostatics=True, only_atoms=nonalchemical_atoms)
     tot_energy_no_nonalchem_particle_electro = compute_energy(system, positions)
     aa_reciprocal_energy = compute_energy(system, positions, force_group={30})
-    turn_off(nonbonded_force, electrostatics=True)
+    turn_off(system, electrostatics=True)
     tot_energy_no_particle_electro = compute_energy(system, positions)
 
     na_reciprocal_energy = tot_reciprocal_energy - nn_reciprocal_energy - aa_reciprocal_energy
@@ -316,13 +284,13 @@ def dissect_nonbonded_energy(reference_system, positions, alchemical_atoms):
     # -----------------------------------------------------
 
     # Compute contributions from exceptions sterics
-    system, nonbonded_force = restore_system(reference_system)  # Restore particle interactions
-    turn_off(nonbonded_force, sterics=True, exceptions=True, only_atoms=alchemical_atoms)
+    system = copy.deepcopy(reference_system)  # Restore particle interactions
+    turn_off(system, sterics=True, exceptions=True, only_atoms=alchemical_atoms)
     tot_energy_no_alchem_exception_sterics = compute_energy(system, positions)
-    system, nonbonded_force = restore_system(reference_system)  # Restore alchemical sterics
-    turn_off(nonbonded_force, sterics=True, exceptions=True, only_atoms=nonalchemical_atoms)
+    system = copy.deepcopy(reference_system)  # Restore alchemical sterics
+    turn_off(system, sterics=True, exceptions=True, only_atoms=nonalchemical_atoms)
     tot_energy_no_nonalchem_exception_sterics = compute_energy(system, positions)
-    turn_off(nonbonded_force, sterics=True, exceptions=True)
+    turn_off(system, sterics=True, exceptions=True)
     tot_energy_no_exception_sterics = compute_energy(system, positions)
 
     tot_exception_sterics = tot_energy - tot_energy_no_exception_sterics
@@ -331,13 +299,13 @@ def dissect_nonbonded_energy(reference_system, positions, alchemical_atoms):
     na_exception_sterics = tot_exception_sterics - nn_exception_sterics - aa_exception_sterics
 
     # Compute contributions from exceptions electrostatics
-    system, nonbonded_force = restore_system(reference_system)  # Restore exceptions sterics
-    turn_off(nonbonded_force, electrostatics=True, exceptions=True, only_atoms=alchemical_atoms)
+    system = copy.deepcopy(reference_system)  # Restore exceptions sterics
+    turn_off(system, electrostatics=True, exceptions=True, only_atoms=alchemical_atoms)
     tot_energy_no_alchem_exception_electro = compute_energy(system, positions)
-    system, nonbonded_force = restore_system(reference_system)  # Restore alchemical electrostatics
-    turn_off(nonbonded_force, electrostatics=True, exceptions=True, only_atoms=nonalchemical_atoms)
+    system = copy.deepcopy(reference_system)  # Restore alchemical electrostatics
+    turn_off(system, electrostatics=True, exceptions=True, only_atoms=nonalchemical_atoms)
     tot_energy_no_nonalchem_exception_electro = compute_energy(system, positions)
-    turn_off(nonbonded_force, electrostatics=True, exceptions=True)
+    turn_off(system, electrostatics=True, exceptions=True)
     tot_energy_no_exception_electro = compute_energy(system, positions)
 
     tot_exception_electro = tot_energy - tot_energy_no_exception_electro
@@ -490,40 +458,6 @@ def compare_system_energies(reference_system, alchemical_system, alchemical_regi
         raise Exception(err_msg.format(delta / unit.kilocalories_per_mole, MAX_DELTA / unit.kilocalories_per_mole))
 
 
-def compare_system_forces(reference_system, alchemical_system, positions, name="", platform=None):
-    """Check that the forces of reference and modified systems are close.
-
-    Parameters
-    ---------
-    reference_system : simtk.openmm.System
-        Reference System
-    alchemical_system : simtk.openmm.System
-        System to compare to reference
-    positions : simtk.unit.Quantity of shape [nparticles,3] with units of distance
-        The particle positions to use
-    name : str, optional, default=""
-        System name to use for debugging.
-    platform : simtk.openmm.Platform, optional, default=None
-        If specified, use this platform
-
-    """
-    # Compute forces
-    reference_force = compute_forces(reference_system, positions, platform=platform) / GLOBAL_FORCE_UNIT
-    alchemical_force = compute_forces(alchemical_system, positions, platform=platform) / GLOBAL_FORCE_UNIT
-
-    # Check that error is small.
-    def magnitude(vec):
-        return np.sqrt(np.mean(np.sum(vec**2, axis=1)))
-
-    relative_error = magnitude(alchemical_force - reference_force) / magnitude(reference_force)
-    if np.any(np.abs(relative_error) > MAX_FORCE_RELATIVE_ERROR):
-        print("========")
-        err_msg = ("Maximum allowable relative force error exceeded (was {:.8f}; allowed {:.8f}).\n"
-                   "alchemical_force = {:.8f}, reference_force = {:.8f}, difference = {:.8f}")
-        raise Exception(err_msg.format(relative_error, MAX_FORCE_RELATIVE_ERROR, magnitude(alchemical_force),
-                                       magnitude(reference_force), magnitude(alchemical_force-reference_force)))
-
-
 def check_interacting_energy_components(reference_system, alchemical_system, alchemical_regions, positions):
     """Compare full and alchemically-modified system energies by energy component.
 
@@ -660,7 +594,7 @@ def check_interacting_energy_components(reference_system, alchemical_system, alc
     for force_name in ['HarmonicBondForce', 'HarmonicAngleForce', 'PeriodicTorsionForce',
                        'GBSAOBCForce', 'CustomGBForce']:
         alchemical_forces_energies = [energy for label, energy in energy_components.items() if force_name in label]
-        reference_force_energy = compute_energy_force(reference_system, positions, force_name)
+        reference_force_energy = compute_force_energy(reference_system, positions, force_name)
 
         # There should be no force in the alchemical system if force_name is missing from the reference
         if reference_force_energy is None:
@@ -795,7 +729,7 @@ def check_noninteracting_energy_components(reference_system, alchemical_system,
                                     if i not in alchemical_regions.alchemical_atoms]
 
         # Compute reference force energy.
-        reference_force_energy = compute_energy_force(system, non_alchemical_positions, force_name)
+        reference_force_energy = compute_force_energy(system, non_alchemical_positions, force_name)
         assert_almost_equal(reference_force_energy, alchemical_energy,
                             'reference {}, alchemical {}'.format(reference_force_energy, alchemical_energy))
 
@@ -1224,12 +1158,7 @@ class TestAbsoluteAlchemicalFactory(object):
     def generate_cases(cls):
         """Generate all test cases in cls.test_cases combinatorially."""
         cls.test_cases = dict()
-        switched_rf_factory = AbsoluteAlchemicalFactory(alchemical_rf_treatment='switched')
-        shifted_rf_factory = AbsoluteAlchemicalFactory(alchemical_rf_treatment='shifted')
-
-        # Create reference versions of all rf-containing systems with their switched counterparts with c_rf = 0
-        for (name, testsystem) in cls.test_systems.items():
-            setattr(testsystem, 'modified_rf_system', switched_rf_factory.replace_reaction_field(testsystem.system))
+        factory = AbsoluteAlchemicalFactory(alchemical_rf_treatment='switched')
 
         # We generate all possible combinations of annihilate_sterics/electrostatics
         # for each test system. We also annihilate bonds, angles and torsions every
@@ -1268,20 +1197,24 @@ class TestAbsoluteAlchemicalFactory(object):
                                              softcore_c=1.0, softcore_d=1.0, softcore_e=1.0, softcore_f=1.0)
                     test_case_name += ', modified softcore parameters'
 
-                # Also store shifted rf alchemical system for energy component comparisons
-                shifted_rf_alchemical_system = shifted_rf_factory.create_alchemical_system(test_system.system, region)
-                setattr(test_system, 'shifted_rf_alchemical_system', shifted_rf_alchemical_system)
+                # Pre-generate alchemical system.
+                alchemical_system = factory.create_alchemical_system(test_system.system, region)
 
-                # Pre-generate alchemical system with switched rf
-                alchemical_system = switched_rf_factory.create_alchemical_system(test_system.system, region)
+                # Add test case.
                 cls.test_cases[test_case_name] = (test_system, alchemical_system, region)
-
                 n_test_cases += 1
 
+            # If the test system uses reaction field replace reaction field
+            # of the reference system to allow comparisons.
+            nonbonded_force = forces.find_nonbonded_force(test_system.system)
+            if nonbonded_force.getNonbondedMethod() == openmm.NonbondedForce.CutoffPeriodic:
+                forcefactories.replace_reaction_field(test_system.system, return_copy=False,
+                                                      switch_width=factory.switch_width)
+
     def test_fully_interacting_energy(self):
         """Compare the energies of reference and fully interacting alchemical system."""
         for test_name, (test_system, alchemical_system, alchemical_region) in self.test_cases.items():
-            f = partial(compare_system_energies, test_system.modified_rf_system,
+            f = partial(compare_system_energies, test_system.system,
                         alchemical_system, alchemical_region, test_system.positions)
             f.description = "Testing fully interacting energy of {}".format(test_name)
             yield f
@@ -1294,25 +1227,6 @@ class TestAbsoluteAlchemicalFactory(object):
             f.description = "Testing non-interacting energy of {}".format(test_name)
             yield f
 
-    def test_replace_reaction_field(self):
-        """Check that replacing reaction-field electrostatics with Custom*Force
-        yields minimal force differences with original system.
-
-        Note that we cannot test for energy consistency or energy overlap because
-        which atoms are within the cutoff will cause energy difference to vary wildly.
-
-        """
-        platform = openmm.Platform.getPlatformByName('Reference')
-        factory = AbsoluteAlchemicalFactory(alchemical_rf_treatment='switched', switch_width=None)
-        for test_name, (test_system, alchemical_system, alchemical_region) in self.test_cases.items():
-            if (test_system.system.getNumForces() != test_system.modified_rf_system.getNumForces()):
-                modified_rf_system = factory.replace_reaction_field(test_system.system)
-                # Make sure positions are not at minimum
-                positions = generate_new_positions(test_system.system, test_system.positions)
-                f = partial(compare_system_forces, test_system.system, modified_rf_system, positions, name=test_name, platform=platform)
-                f.description = "Testing replace_reaction_field on system {}".format(test_name)
-                yield f
-
     @attr('slow')
     def test_fully_interacting_energy_components(self):
         """Test interacting state energy by force component."""
@@ -1322,10 +1236,8 @@ class TestAbsoluteAlchemicalFactory(object):
                             if 'Explicit' in test_name]
         for test_name in test_cases_names:
             test_system, alchemical_system, alchemical_region = self.test_cases[test_name]
-            # We have to compare shifted rf system with original system because test cannot handle
-            # re-coded reaction-field forces with c_rf = 0
-            f = partial(check_interacting_energy_components, test_system.system, test_system.shifted_rf_alchemical_system,
-                    alchemical_region, test_system.positions)
+            f = partial(check_interacting_energy_components, test_system.system, alchemical_system,
+                        alchemical_region, test_system.positions)
             f.description = "Testing energy components of %s..." % test_name
             yield f
 
@@ -1347,10 +1259,12 @@ class TestAbsoluteAlchemicalFactory(object):
         for platform in platforms:
             GLOBAL_ALCHEMY_PLATFORM = platform
             for test_name, (test_system, alchemical_system, alchemical_region) in self.test_cases.items():
-                f = partial(compare_system_energies, test_system.modified_rf_system, alchemical_system, alchemical_region, test_system.positions)
+                f = partial(compare_system_energies, test_system.system, alchemical_system,
+                            alchemical_region, test_system.positions)
                 f.description = "Test fully interacting energy of {} on {}".format(test_name, platform.getName())
                 yield f
-                f = partial(check_noninteracting_energy_components, alchemical_system, alchemical_region, test_system.positions)
+                f = partial(check_noninteracting_energy_components, test_system.system, alchemical_system,
+                            alchemical_region, test_system.positions)
                 f.description = "Test non-interacting energy of {} on {}".format(test_name, platform.getName())
                 yield f
 
@@ -1364,7 +1278,7 @@ class TestAbsoluteAlchemicalFactory(object):
             #cached_trajectory_filename = os.path.join(os.environ['HOME'], '.cache', 'alchemy', 'tests',
             #                                           test_name + '.pickle')
             cached_trajectory_filename = None
-            f = partial(overlap_check, test_system.modified_rf_system, alchemical_system, test_system.positions,
+            f = partial(overlap_check, test_system.system, alchemical_system, test_system.positions,
                         cached_trajectory_filename=cached_trajectory_filename, name=test_name)
             f.description = "Testing reference/alchemical overlap for {}".format(test_name)
             yield f
diff --git a/openmmtools/tests/test_cache.py b/openmmtools/tests/test_cache.py
index 8704183..ef8a8ef 100644
--- a/openmmtools/tests/test_cache.py
+++ b/openmmtools/tests/test_cache.py
@@ -260,19 +260,9 @@ class TestContextCache(object):
         assert state1.is_context_compatible(context)
         assert isinstance(integrator, type(self.verlet_2fs)), type(integrator)
 
-        # When it has a choice, ContextCache picks the same context
-        # in consecutive calls with same thermodynamic state.
-        # First add another integrator so that ContextCache has 2 possible options.
-        assert type(self.langevin_2fs_310k) is not type(default_integrator)  # test precondition
-        cache.get_context(state1, copy.deepcopy(self.langevin_2fs_310k))
-        assert len(cache) == 3
-        context, integrator = cache.get_context(state1)
-        for i in range(5):  # 5 attempts to make the test fail
-            assert cache.get_context(state1)[0] is context
-
         # With an incompatible state, a new Context is created.
         cache.get_context(state2)
-        assert len(cache) == 4
+        assert len(cache) == 3
 
     def test_cache_capacity_ttl(self):
         """Check that the cache capacity and time_to_live work as expected."""
diff --git a/openmmtools/tests/test_forces.py b/openmmtools/tests/test_forces.py
new file mode 100644
index 0000000..d076da6
--- /dev/null
+++ b/openmmtools/tests/test_forces.py
@@ -0,0 +1,169 @@
+#!/usr/bin/env python
+
+# =============================================================================
+# MODULE DOCSTRING
+# =============================================================================
+
+"""
+Test Force classes in forces.py.
+
+"""
+
+# =============================================================================
+# GLOBAL IMPORTS
+# =============================================================================
+
+from functools import partial
+
+import numpy as np
+
+from openmmtools.forcefactories import *
+from openmmtools import testsystems
+
+
+# =============================================================================
+# CONSTANTS
+# =============================================================================
+
+MAX_FORCE_RELATIVE_ERROR = 1.0e-6  # maximum allowable relative force error
+GLOBAL_FORCE_UNIT = unit.kilojoules_per_mole / unit.nanometers  # controls printed units
+GLOBAL_FORCES_PLATFORM = None  # This is used in every calculation.
+
+
+# =============================================================================
+# TESTING UTILITIES
+# =============================================================================
+
+def create_context(system, integrator, platform=None):
+    """Create a Context.
+
+    If platform is None, GLOBAL_ALCHEMY_PLATFORM is used.
+
+    """
+    if platform is None:
+        platform = GLOBAL_FORCES_PLATFORM
+    if platform is not None:
+        context = openmm.Context(system, integrator, platform)
+    else:
+        context = openmm.Context(system, integrator)
+    return context
+
+
+# =============================================================================
+# UTILITY FUNCTIONS
+# =============================================================================
+
+def compute_forces(system, positions, platform=None, force_group=-1):
+    """Compute forces of the system in the given positions.
+
+    Parameters
+    ----------
+    platform : simtk.openmm.Platform or None, optional
+        If None, the global GLOBAL_ALCHEMY_PLATFORM will be used.
+    force_group : int flag or set of int, optional
+        Passed to the groups argument of Context.getState().
+
+    """
+    timestep = 1.0 * unit.femtoseconds
+    integrator = openmm.VerletIntegrator(timestep)
+    context = create_context(system, integrator, platform)
+    context.setPositions(positions)
+    state = context.getState(getForces=True, groups=force_group)
+    forces = state.getForces(asNumpy=True)
+    del context, integrator, state
+    return forces
+
+
+def compare_system_forces(reference_system, alchemical_system, positions, name="", platform=None):
+    """Check that the forces of reference and modified systems are close.
+
+    Parameters
+    ---------
+    reference_system : simtk.openmm.System
+        Reference System
+    alchemical_system : simtk.openmm.System
+        System to compare to reference
+    positions : simtk.unit.Quantity of shape [nparticles,3] with units of distance
+        The particle positions to use
+    name : str, optional, default=""
+        System name to use for debugging.
+    platform : simtk.openmm.Platform, optional, default=None
+        If specified, use this platform
+
+    """
+    # Compute forces
+    reference_force = compute_forces(reference_system, positions, platform=platform) / GLOBAL_FORCE_UNIT
+    alchemical_force = compute_forces(alchemical_system, positions, platform=platform) / GLOBAL_FORCE_UNIT
+
+    # Check that error is small.
+    def magnitude(vec):
+        return np.sqrt(np.mean(np.sum(vec**2, axis=1)))
+
+    relative_error = magnitude(alchemical_force - reference_force) / magnitude(reference_force)
+    if np.any(np.abs(relative_error) > MAX_FORCE_RELATIVE_ERROR):
+        err_msg = ("Maximum allowable relative force error exceeded (was {:.8f}; allowed {:.8f}).\n"
+                   "alchemical_force = {:.8f}, reference_force = {:.8f}, difference = {:.8f}")
+        raise Exception(err_msg.format(relative_error, MAX_FORCE_RELATIVE_ERROR, magnitude(alchemical_force),
+                                       magnitude(reference_force), magnitude(alchemical_force-reference_force)))
+
+
+def generate_new_positions(system, positions, platform=None, nsteps=50):
+    """Generate new positions by taking a few steps from the old positions.
+
+    Parameters
+    ----------
+    platform : simtk.openmm.Platform or None, optional
+        If None, the global GLOBAL_ALCHEMY_PLATFORM will be used.
+    nsteps : int, optional, default=50
+        Number of steps of dynamics to take.
+
+    Returns
+    -------
+    new_positions : simtk.unit.Quantity of shape [nparticles,3] with units compatible with distance
+        New positions
+
+    """
+    temperature = 300 * unit.kelvin
+    collision_rate = 90 / unit.picoseconds
+    timestep = 1.0 * unit.femtoseconds
+    integrator = openmm.LangevinIntegrator(temperature, collision_rate, timestep)
+    context = create_context(system, integrator, platform)
+    context.setPositions(positions)
+    integrator.step(nsteps)
+    new_positions = context.getState(getPositions=True).getPositions(asNumpy=True)
+    del context, integrator
+    return new_positions
+
+
+# =============================================================================
+# TEST LRU CACHE
+# =============================================================================
+
+def test_replace_reaction_field():
+    """Check that replacing reaction-field electrostatics with Custom*Force
+    yields minimal force differences with original system.
+
+    Note that we cannot test for energy consistency or energy overlap because
+    which atoms are within the cutoff will cause energy difference to vary wildly.
+
+    """
+    test_cases = [
+        testsystems.AlanineDipeptideExplicit(nonbondedMethod=openmm.app.CutoffPeriodic),
+        testsystems.HostGuestExplicit(nonbondedMethod=openmm.app.CutoffPeriodic)
+    ]
+    platform = openmm.Platform.getPlatformByName('Reference')
+    for test_system in test_cases:
+        test_name = test_system.__class__.__name__
+
+        # Replace reaction field.
+        modified_rf_system = replace_reaction_field(test_system.system, switch_width=None)
+
+        # Make sure positions are not at minimum.
+        positions = generate_new_positions(test_system.system, test_system.positions)
+
+        # Test forces.
+        f = partial(compare_system_forces, test_system.system, modified_rf_system, positions,
+                    name=test_name, platform=platform)
+        f.description = "Testing replace_reaction_field on system {}".format(test_name)
+        yield f
+