Skip to content

Fix bug when computing reduced potential in multiple states #385

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 3 commits into from
Nov 9, 2018
Merged

Fix bug when computing reduced potential in multiple states #385

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
8188264
Fix bug when computing reduced potential in multiple states.
andrrizzi Nov 9, 2018
37aaa34
Use logger.warning instead of print
andrrizzi Nov 9, 2018
86c4a18
Update release history
andrrizzi Nov 9, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 9 additions & 1 deletion docs/releasehistory.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,9 +16,17 @@ considerable speed improvement over the previous implementation (`#380 <https://
- It is now possible to have multiple composable states exposing the same attributes/getter/setter in a
``CompoundThermodynamicState`` (`#380 <https://github.com/choderalab/openmmtools/pull/380>`_).

Bug fixes
---------
- Fixed a bug involving the ``NoseHooverChainVelocityVerletIntegrator`` with ``System`` with constraints. The constraints
were not taken into account when calculating the number of degrees of freedom resulting in the temperature not converging
to the target value. (`#384 <https://github.com/choderalab/openmmtools/pull/384>`_)
- Fixed a bug affecting ``reduced_potential_at_states`` when computing the reduced potential of systems in different
``AlchemicalState``s when the same alchemical parameter appeared in force objects split in different force groups. (`#385 <https://github.com/choderalab/openmmtools/pull/385>`_)

Deprecated and API breaks
-------------------------
- Python 2, which was deprecated in 0.16.0, is now unsupported.
- Python 2 is not supported anymore.
- The ``update_alchemical_charges`` attribute of ``AlchemicalState`, which was deprecated in 0.16.0, has now been removed
since it doesn't make sense with the new parameter offset implementation.
- The methods ``AlchemicalState.get_alchemical_variable`` and ``AlchemicalState.set_alchemical_variable`` have been
Expand Down
12 changes: 6 additions & 6 deletions openmmtools/integrators.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -630,8 +630,9 @@ def __init__(self, system=None, temperature=298*unit.kelvin, collision_frequency

system: openmm.app.System instance
Required to extract the system's number of degrees of freedom.
If the system is not passed to the constructor,
the temperature will converge to the wrong value.
If the system is not passed to the constructor and has constraints
or a ``CMMotionRemover`` force, the temperature will converge to the
wrong value.

temperature: unit.Quantity compatible with kelvin, default=298*unit.kelvin
The target temperature for the thermostat.
Expand Down Expand Up @@ -676,7 +677,7 @@ def __init__(self, system=None, temperature=298*unit.kelvin, collision_frequency
if chain_length < 0:
raise Exception("Nosé-Hoover chain length must be at least 0")
if chain_length == 0:
print("WARNING: Nosé-Hoover chain length is 0; falling back to regular velocity verlet algorithm.")
logger.warning('Nosé-Hoover chain length is 0; falling back to regular velocity verlet algorithm.')
self.M = chain_length

# Define the "mass" of the thermostat particles (multiply by ndf for particle 0)
Expand All @@ -688,8 +689,8 @@ def __init__(self, system=None, temperature=298*unit.kelvin, collision_frequency
# Compute the number of degrees of freedom.
#
if system is None:
print("SEVERE WARNING: The system was not passed to the NoseHooverChainVelocityVerletIntegrator. "
"For systems with constraints, the simulation will run at the wrong temperature.")
logger.warning('The system was not passed to the NoseHooverChainVelocityVerletIntegrator. '
'For systems with constraints, the simulation will run at the wrong temperature.')
# Fall back to old scheme, which only works for unconstrained systems
self.addGlobalVariable("ndf", 0)
self.addPerDofVariable("ones", 1.0)
Expand Down Expand Up @@ -750,7 +751,6 @@ def __init__(self, system=None, temperature=298*unit.kelvin, collision_frequency
# Compute heat bath energies
self.computeEnergies()


def propagateNHC(self):
""" Propagate the Nosé-Hoover chain """
self.addComputeGlobal("scale", "1.0")
Expand Down
20 changes: 9 additions & 11 deletions openmmtools/states.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3435,27 +3435,25 @@ def _find_force_groups_to_update(self, context, current_context_state, memo):
`current_context_state`.
"""
# Cache information about system force groups.
# We create a dictionary "memo" mapping parameter_name -> list of force groups to update.
if len(memo) == 0:
parameters_found = set()
system = context.getSystem()
system_parameters = self._get_system_controlled_parameters(system, self._parameters_name_suffix)
for force, parameter_name, _ in system_parameters:
if parameter_name not in parameters_found:
parameters_found.add(parameter_name)
# Keep track of valid parameters only.
if self._parameters[parameter_name] is not None:
memo[parameter_name] = force.getForceGroup()
# Break the loop if we have found all the parameters.
if len(parameters_found) == len(self._parameters):
break
# Keep track of valid parameters only.
if self._parameters[parameter_name] is not None:
try:
memo[parameter_name].append(force.getForceGroup())
except KeyError:
memo[parameter_name] = [force.getForceGroup()]

# Find lambda parameters that will change.
force_groups_to_update = set()
for parameter_name, force_group in memo.items():
for parameter_name, force_groups in memo.items():
self_parameter_value = getattr(self, parameter_name)
current_parameter_value = getattr(current_context_state, parameter_name)
if self_parameter_value != current_parameter_value:
force_groups_to_update.add(force_group)
force_groups_to_update.update(force_groups)
return force_groups_to_update

# -------------------------------------------------------------------------
Expand Down
7 changes: 3 additions & 4 deletions openmmtools/tests/test_integrators.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -224,10 +224,9 @@ def test_nose_hoover_integrator():

# Coarse check for target temperature
mean_temperature = np.mean(temperatures)
print(mean_temperature)
assert abs(mean_temperature - temperature.value_in_unit(unit.kelvin)) < 10.0


assert abs(mean_temperature - temperature.value_in_unit(unit.kelvin)) < 10.0, mean_temperature


def test_pretty_formatting():
"""
Test pretty-printing and pretty-formatting of integrators.
Expand Down
35 changes: 26 additions & 9 deletions openmmtools/tests/test_states.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1358,11 +1358,13 @@ def setup_class(cls):
system = openmm.System()
system.addParticle(40.0*unit.amu)
system.addParticle(40.0*unit.amu)
# Add a force defining lambda_bonds and gamma global parameters.
custom_force = openmm.CustomBondForce('lambda_bonds^gamma*60000*(r-{})^2;'.format(r0_nanometers))
custom_force.addGlobalParameter('lambda_bonds', cls.parameters_default_values['lambda_bonds'])
custom_force.addGlobalParameter('gamma', cls.parameters_default_values['gamma'])
custom_force.addBond(0, 1, [])
system.addForce(custom_force)
# Add a force defining the lambda_bonds_mysuffix global parameters.
custom_force_suffix = openmm.CustomBondForce('lambda_bonds_mysuffix*20000*(r-{})^2;'.format(r0_nanometers))
custom_force_suffix.addGlobalParameter('lambda_bonds_mysuffix', cls.parameters_default_values['lambda_bonds_mysuffix'])
custom_force_suffix.addBond(0, 1, [])
Expand All @@ -1374,6 +1376,14 @@ def setup_class(cls):
pos2 = [0.0, 0.0, r0_nanometers]
cls.diatomic_molecule_ss = SamplerState(positions=np.array([pos1, pos2]) * unit.nanometers)

# Create a system with a duplicate force to test handling forces
# defining the same parameters in different force groups.
custom_force = copy.deepcopy(custom_force_suffix)
custom_force.setForceGroup(30)
system_force_groups = copy.deepcopy(system)
system_force_groups.addForce(custom_force)
cls.diatomic_molecule_force_groups_ts = ThermodynamicState(system_force_groups, temperature=300.0*unit.kelvin)

# Create few incompatible systems for testing. An incompatible state
# has a different set of defined global parameters.
cls.incompatible_systems = [system]
Expand All @@ -1383,6 +1393,9 @@ def setup_class(cls):
cls.incompatible_systems.append(copy.deepcopy(system))
cls.incompatible_systems[i+1].removeForce(i)

# System with the global parameters duplicated in two different force groups.
cls.incompatible_systems.append(copy.deepcopy(system_force_groups))

# System with both lambda_bonds_suffix and gamma_bond_suffix defined (instead of only the former).
cls.incompatible_systems.append(copy.deepcopy(system))
custom_force = copy.deepcopy(cls.incompatible_systems[-1].getForce(1))
Expand Down Expand Up @@ -1664,19 +1677,23 @@ def check_is_standard(states, is_standard):

def test_find_force_groups_to_update(self):
"""Test method GlobalParameterState._find_force_groups_to_update."""
system = self.diatomic_molecule_ts.system
system = self.diatomic_molecule_force_groups_ts.system
integrator = openmm.VerletIntegrator(2.0*unit.femtoseconds)
# Test cases are (force_group, force_group_suffix)
test_cases = [(0, 0), (1, 5), (9, 4)]
# Test cases are (force_groups, force_groups_suffix)
test_cases = [
([0], [0, 0]),
([1], [5, 5]),
([9], [4, 2])
]

for force_groups in test_cases:
for i, force_group in enumerate(force_groups):
for test_case in test_cases:
for i, force_group in enumerate(test_case[0] + test_case[1]):
system.getForce(i).setForceGroup(force_group)
states = self.read_system_state(system)
context = openmm.Context(system, copy.deepcopy(integrator))

# No force group should be updated if we don't change the global parameter.
for state, force_group in zip(states, force_groups):
for state, force_groups in zip(states, test_case):
assert state._find_force_groups_to_update(context, state, memo={}) == set()

# Change the lambdas one by one and check that the method
Expand All @@ -1690,7 +1707,7 @@ def test_find_force_groups_to_update(self):

# Change the current state.
setattr(current_state, parameter_name, parameter_value / 2)
assert state._find_force_groups_to_update(context, current_state, memo={}) == {force_group}
assert state._find_force_groups_to_update(context, current_state, memo={}) == set(force_groups)
setattr(current_state, parameter_name, parameter_value) # Reset current state.
del context

Expand Down Expand Up @@ -1876,8 +1893,8 @@ def test_reduced_potential_compound_state(self):
positions = copy.deepcopy(self.diatomic_molecule_ss.positions)
# Build a mixed collection of compatible and incompatible thermodynamic states.
thermodynamic_states = [
ThermodynamicState(self.incompatible_systems[0], temperature=300*unit.kelvin),
ThermodynamicState(self.incompatible_systems[-1], temperature=300*unit.kelvin)
copy.deepcopy(self.diatomic_molecule_ts),
copy.deepcopy(self.diatomic_molecule_force_groups_ts)
]

compound_states = []
Expand Down