choderalab · andrrizzi · Sep 17, 2018 · Jul 21, 2018 · Jul 21, 2018 · Jul 21, 2018
diff --git a/docs/releasehistory.rst b/docs/releasehistory.rst
@@ -1,13 +1,29 @@
 Release History
 ===============
 
-0.15.1 - Something New This Way Comes [WIP]
-===========================================
+0.16.0 - GlobalParameterState class, SamplerState reads CVs
+============================================================
 
 New features
 ------------
-- Add ability for ``SamplerState`` to access new `OpenMM Custom CV Force Variables <http://docs.openmm.org/development/api-python/generated/simtk.openmm.openmm.CustomCVForce.html#simtk.openmm.openmm.CustomCVForce>`_
-- ``SamplerState.update_from_context`` now has keywords to support finer grain updating from the Context. This is only recommended for advanced users.
+- Add ability for ``SamplerState`` to access new `OpenMM Custom CV Force Variables
+  <http://docs.openmm.org/development/api-python/generated/simtk.openmm.openmm.CustomCVForce.html#simtk.openmm.openmm.CustomCVForce>`_ 
+  (`#362 <https://github.com/choderalab/openmmtools/pull/362>`_).
+- ``SamplerState.update_from_context`` now has keywords to support finer grain updating from the Context. This is only
+  recommended for advanced users (`#362 <https://github.com/choderalab/openmmtools/pull/362>`_).
+- Added the new class ``states.GlobalParameterState`` designed to simplify the implementation of composable states that
+  control global variables (`#363 <https://github.com/choderalab/openmmtools/pull/363>`_).
+- Allow restraint force classes to be controlled by a parameter other than ``lambda_restraints``. This will enable
+  multi-restraints simulations (`#363 <https://github.com/choderalab/openmmtools/pull/363>`_).
+
+Deprecated
+----------
+- Python2 is officially deprecated. Support will be dropped in future versions.
+- Deprecated the signature of ``IComposableState._on_setattr`` to fix a bug where the objects were temporarily left in
+  an inconsistent state when an exception was raised and caught.
+- Deprecated ``update_alchemical_charges`` in ``AlchemicalState`` in anticipation of the new implementation of the
+  exact PME that will be based on the ``NonbondedForce`` offsets rather than ``updateParametersInContext()``.
+
 
 0.15.0 - Restraint forces
 =========================

diff --git a/docs/states.rst b/docs/states.rst
@@ -8,6 +8,7 @@ The module :mod:`openmmtools.states` contains classes to maintain a consistent s
  - :class:`ThermodynamicState`: Represent and manipulate the thermodynamic state of OpenMM `System <http://docs.openmm.org/7.1.0/api-python/generated/simtk.openmm.openmm.System.html>`_ and `Context <http://docs.openmm.org/7.1.0/api-python/generated/simtk.openmm.openmm.Context.html>`_ objects.
  - :class:`SamplerState`: Represent and cache the state of the simulation that changes when the `System <http://docs.openmm.org/7.1.0/api-python/generated/simtk.openmm.openmm.System.html>`_ is integrated.
  - :class:`CompoundThermodynamicState`: Extend the :class:`ThermodynamicState` to handle parameters other than temperature and pressure through the implementations of the :class:`IComposableState` abstract class.
+ - :class:`GlobalParameterState`: An implementation of :class:`IComposableState` specialized to control OpenMM forces' global parameters.
 
 States
 ------
@@ -21,3 +22,4 @@ States
     SamplerState
     CompoundThermodynamicState
     IComposableState
+    GlobalParameterState
diff --git a/openmmtools/alchemy.py b/openmmtools/alchemy.py
@@ -117,7 +117,7 @@ class AlchemicalState(object):
     lambda_torsions : float, optional
         Scaling factor for alchemically-softened torsions (default is 1.0).
     update_alchemical_charges : bool, optional
-        If True, ``lambda_electrostatics`` changes in alchemical systems
+        If False, ``lambda_electrostatics`` changes in alchemical systems
         that use exact treatment of PME electrostatics will be considered
         incompatible. This means that a new ``Context`` will be required
         for each `lambda_electrostatics`` state.
@@ -211,6 +211,12 @@ class AlchemicalState(object):
     # -------------------------------------------------------------------------
 
     def __init__(self, **kwargs):
+        if 'update_alchemical_charges' in kwargs:
+            # TODO Drop support for the parameter and remove deprecation warning from 0.17 on
+            # TODO after implementing new way for exact PME based on the NonbondedForce offsets.
+            import warnings
+            warnings.warn('The update_alchemical_charges in AlchemicalState.__init__ has been '
+                          'deprecated, and future versions of openmmtools may not support it.')
         self._initialize(**kwargs)
 
     @classmethod
@@ -251,16 +257,19 @@ def from_system(cls, system):
                 alchemical_parameters[parameter_name] = parameter_value
 
         # Handle the update parameters flag.
-        update_alchemical_charges = bool(alchemical_parameters.pop(_UPDATE_ALCHEMICAL_CHARGES_PARAMETER,
-                                                                   cls._UPDATE_ALCHEMICAL_CHARGES_DEFAULT))
+        update_alchemical_charges = alchemical_parameters.pop(_UPDATE_ALCHEMICAL_CHARGES_PARAMETER, None)
 
         # Check that the system is alchemical.
         if len(alchemical_parameters) == 0:
             raise AlchemicalStateError('System has no lambda parameters.')
 
+        # Avoid passing update_alchemical_charges if not necessary to not raise deprecation warning.
+        if (update_alchemical_charges is not None and
+                    update_alchemical_charges != cls._UPDATE_ALCHEMICAL_CHARGES_DEFAULT):
+            alchemical_parameters['update_alchemical_charges'] = bool(update_alchemical_charges)
+
         # Create and return the AlchemicalState.
-        return AlchemicalState(update_alchemical_charges=update_alchemical_charges,
-                               **alchemical_parameters)
+        return AlchemicalState(**alchemical_parameters)
 
     # -------------------------------------------------------------------------
     # Lambda properties
@@ -392,7 +401,8 @@ def __setstate__(self, serialization):
         parameters = serialization['parameters']
         alchemical_variables = serialization['alchemical_variables']
         # New attribute in OpenMMTools 0.14.0.
-        update_alchemical_charges = serialization.get('update_alchemical_charges', True)
+        update_alchemical_charges = serialization.get('update_alchemical_charges',
+                                                      self._UPDATE_ALCHEMICAL_CHARGES_DEFAULT)
         alchemical_functions = dict()
 
         # Temporarily store alchemical functions.
@@ -554,7 +564,7 @@ def _standardize_system(self, system, set_lambda_electrostatics=False):
         # states with different settings must be incompatible.
         alchemical_state._apply_to_system(system, set_update_charges_flag=False)
 
-    def _on_setattr(self, standard_system, attribute_name):
+    def _on_setattr(self, standard_system, attribute_name, old_attribute_value):
         """Check if the standard system needs changes after a state attribute is set.
 
         Parameters
@@ -563,6 +573,8 @@ def _on_setattr(self, standard_system, attribute_name):
             The standard system before setting the attribute.
         attribute_name : str
             The name of the attribute that has just been set or retrieved.
+        old_attribute_value : float
+            The value of the attribute retrieved before being set.
 
         Returns
         -------

diff --git a/openmmtools/forcefactories.py b/openmmtools/forcefactories.py
@@ -50,7 +50,7 @@ def replace_reaction_field(reference_system, switch_width=1.0*unit.angstrom,
         Setting it to `False` speeds up the function execution but modifies
         the `reference_system` object.
     shifted : bool, optional, default=False
-        If `True`, a shited reaction-field will be used.
+        If `True`, a shifted reaction-field will be used.
 
     Returns
     -------

diff --git a/openmmtools/forces.py b/openmmtools/forces.py
@@ -254,14 +254,22 @@ class RadiallySymmetricRestraintForce(utils.RestorableOpenMMObject):
         The indices of the first group of atoms to restrain.
     restrained_atom_indices2 : iterable of int
         The indices of the second group of atoms to restrain.
+    controlling_parameter_name : str
+        The name of the global parameter controlling the energy function.
     *args, **kwargs
         Parameters to pass to the super constructor.
 
+    Attributes
+    ----------
+    controlling_parameter_name
+
     """
 
     def __init__(self, restraint_parameters, restrained_atom_indices1,
-                 restrained_atom_indices2, *args, **kwargs):
+                 restrained_atom_indices2, controlling_parameter_name,
+                 *args, **kwargs):
         super(RadiallySymmetricRestraintForce, self).__init__(*args, **kwargs)
+        self._controlling_parameter_name = controlling_parameter_name
 
         # Unzip bond parameters names and values from dict.
         assert len(restraint_parameters) == 1 or isinstance(restraint_parameters, collections.OrderedDict)
@@ -271,7 +279,7 @@ def __init__(self, restraint_parameters, restrained_atom_indices1,
         self._create_bond(parameter_values, restrained_atom_indices1, restrained_atom_indices2)
 
         # Add parameters.
-        self.addGlobalParameter('lambda_restraints', 1.0)
+        self.addGlobalParameter(self._controlling_parameter_name, 1.0)
         for parameter in parameter_names:
             self.addPerBondParameter(parameter)
 
@@ -318,6 +326,11 @@ def restraint_parameters(self):
                                 for parameter_idx, parameter_value in enumerate(parameter_values)]
         return collections.OrderedDict(restraint_parameters)
 
+    @property
+    def controlling_parameter_name(self):
+        """str: The name of the global parameter controlling the energy function (read-only)."""
+        return self._controlling_parameter_name
+
     def distance_at_energy(self, potential_energy):
         """Compute the distance at which the potential energy is ``potential_energy``.
 
@@ -650,16 +663,17 @@ class RadiallySymmetricCentroidRestraintForce(RadiallySymmetricRestraintForce,
 
     The restraint is applied between the centers of mass of two groups
     of atoms. The restraint strength is controlled by a global context
-    parameter called 'lambda_restraints'.
+    parameter whose name is passed on construction through the optional
+    argument ``controlling_parameter_name``.
 
     With OpenCL, only on 64bit platforms are supported.
 
     Parameters
     ----------
     energy_function : str
         The energy function to pass to ``CustomCentroidBondForce``. The
-        global parameter 'lambda_restraint' will be prepended to this
-        expression.
+        name of the controlling global parameter  will be prepended to
+        this expression.
     restraint_parameters : OrderedDict
         An ordered dictionary containing the bond parameters in the form
         parameter_name: parameter_value. The order is important to make
@@ -669,23 +683,28 @@ class RadiallySymmetricCentroidRestraintForce(RadiallySymmetricRestraintForce,
         The indices of the first group of atoms to restrain.
     restrained_atom_indices2 : iterable of int
         The indices of the second group of atoms to restrain.
+    controlling_parameter_name : str, optional
+        The name of the global parameter controlling the energy function.
+        The default value is 'lambda_restraints'.
 
     Attributes
     ----------
     restraint_parameters
     restrained_atom_indices1
     restrained_atom_indices2
+    controlling_parameter_name
 
     """
 
     def __init__(self, energy_function, restraint_parameters,
-                 restrained_atom_indices1, restrained_atom_indices2):
+                 restrained_atom_indices1, restrained_atom_indices2,
+                 controlling_parameter_name='lambda_restraints'):
         # Initialize CustomCentroidBondForce.
-        energy_function = 'lambda_restraints * ' + energy_function
+        energy_function = controlling_parameter_name + ' * ' + energy_function
         custom_centroid_bond_force_args = [2, energy_function]
         super(RadiallySymmetricCentroidRestraintForce, self).__init__(
             restraint_parameters, restrained_atom_indices1, restrained_atom_indices2,
-            *custom_centroid_bond_force_args)
+            controlling_parameter_name, *custom_centroid_bond_force_args)
 
     @property
     def restrained_atom_indices1(self):
@@ -726,12 +745,14 @@ class RadiallySymmetricBondRestraintForce(RadiallySymmetricRestraintForce,
     """
 
     def __init__(self, energy_function, restraint_parameters,
-                 restrained_atom_index1, restrained_atom_index2):
+                 restrained_atom_index1, restrained_atom_index2,
+                 controlling_parameter_name='lambda_restraints'):
         # Initialize CustomBondForce.
         energy_function = energy_function.replace('distance(g1,g2)', 'r')
-        energy_function = 'lambda_restraints * ' + energy_function
-        super(RadiallySymmetricBondRestraintForce, self).__init__(restraint_parameters,
-                 [restrained_atom_index1], [restrained_atom_index2], energy_function)
+        energy_function = controlling_parameter_name + ' * ' + energy_function
+        super(RadiallySymmetricBondRestraintForce, self).__init__(
+            restraint_parameters, [restrained_atom_index1], [restrained_atom_index2],
+            controlling_parameter_name, energy_function)
 
     # -------------------------------------------------------------------------
     # Public properties.
@@ -832,12 +853,11 @@ class HarmonicRestraintForce(HarmonicRestraintForceMixIn,
 
     The energy expression of the restraint is given by
 
-       ``E = lambda_restraints * (K/2)*r^2``
+       ``E = controlling_parameter * (K/2)*r^2``
 
     where `K` is the spring constant, `r` is the distance between the
-    two group centroids, and `lambda_restraints` is a scale factor that
-    can be used to control the strength of the restraint. You can control
-    ``lambda_restraints`` through :class:`RestraintState` class.
+    two group centroids, and `controlling_parameter` is a scale factor that
+    can be used to control the strength of the restraint.
 
     With OpenCL, only on 64bit platforms are supported.
 
@@ -850,13 +870,17 @@ class HarmonicRestraintForce(HarmonicRestraintForceMixIn,
         The indices of the first group of atoms to restrain.
     restrained_atom_indices2 : iterable of int
         The indices of the second group of atoms to restrain.
+    controlling_parameter_name : str, optional
+        The name of the global parameter controlling the energy function.
+        The default value is 'lambda_restraints'.
 
     Attributes
     ----------
     spring_constant
     restrained_atom_indices1
     restrained_atom_indices2
     restraint_parameters
+    controlling_parameter_name
 
     """
     # All the methods are provided by the mix-ins.
@@ -879,13 +903,17 @@ class HarmonicRestraintBondForce(HarmonicRestraintForceMixIn,
         The index of the first atom to restrain.
     restrained_atom_index2 : int
         The index of the second atom to restrain.
+    controlling_parameter_name : str, optional
+        The name of the global parameter controlling the energy function.
+        The default value is 'lambda_restraints'.
 
     Attributes
     ----------
     spring_constant
     restrained_atom_indices1
     restrained_atom_indices2
     restraint_parameters
+    controlling_parameter_name
 
     """
     # All the methods are provided by the mix-ins.
@@ -986,14 +1014,13 @@ class FlatBottomRestraintForce(FlatBottomRestraintForceMixIn,
 
     More precisely, the energy expression of the restraint is given by
 
-        ``E = lambda_restraints * step(r-r0) * (K/2)*(r-r0)^2``
+        ``E = controlling_parameter * step(r-r0) * (K/2)*(r-r0)^2``
 
     where ``K`` is the spring constant, ``r`` is the distance between the
     restrained atoms, ``r0`` is another parameter defining the distance
-    at which the restraint is imposed, and ``lambda_restraints``
+    at which the restraint is imposed, and ``controlling_parameter``
     is a scale factor that can be used to control the strength of the
-    restraint. You can control ``lambda_restraints`` through the class
-    :class:`RestraintState`.
+    restraint.
 
     With OpenCL, only on 64bit platforms are supported.
 
@@ -1009,6 +1036,9 @@ class FlatBottomRestraintForce(FlatBottomRestraintForceMixIn,
         The indices of the first group of atoms to restrain.
     restrained_atom_indices2 : iterable of int
         The indices of the second group of atoms to restrain.
+    controlling_parameter_name : str, optional
+        The name of the global parameter controlling the energy function.
+        The default value is 'lambda_restraints'.
 
     Attributes
     ----------
@@ -1017,6 +1047,7 @@ class FlatBottomRestraintForce(FlatBottomRestraintForceMixIn,
     restrained_atom_indices1
     restrained_atom_indices2
     restraint_parameters
+    controlling_parameter_name
 
     """
     # All the methods are provided by the mix-ins.
@@ -1042,6 +1073,9 @@ class FlatBottomRestraintBondForce(FlatBottomRestraintForceMixIn,
         The index of the first group of atoms to restrain.
     restrained_atom_index2 : int
         The index of the second group of atoms to restrain.
+    controlling_parameter_name : str, optional
+        The name of the global parameter controlling the energy function.
+        The default value is 'lambda_restraints'.
 
     Attributes
     ----------
@@ -1050,6 +1084,7 @@ class FlatBottomRestraintBondForce(FlatBottomRestraintForceMixIn,
     restrained_atom_indices1
     restrained_atom_indices2
     restraint_parameters
+    controlling_parameter_name
 
     """
     # All the methods are provided by the mix-ins.