From b11fc6a6fe8cfea632ef0fdad831c3d446e02c8b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20Kr=C3=A4mer?= Date: Tue, 24 Sep 2019 11:55:45 -0400 Subject: [PATCH 01/12] ThermodynamicState support for more barostats - supports MonteCarloAnistropicBarostat with same pressure along all axes - supports MonteCarloMembraneBarostat with zero surface tension --- openmmtools/states.py | 67 +++++++++++++--- openmmtools/tests/test_states.py | 126 +++++++++++++++++++++++-------- 2 files changed, 151 insertions(+), 42 deletions(-) diff --git a/openmmtools/states.py b/openmmtools/states.py index aad42e21c..0456aa55f 100644 --- a/openmmtools/states.py +++ b/openmmtools/states.py @@ -281,11 +281,13 @@ class ThermodynamicsError(Exception): MULTIPLE_BAROSTATS, NO_BAROSTAT, UNSUPPORTED_BAROSTAT, + UNSUPPORTED_ANISOTROPIC_BAROSTAT, + UNSUPPORTED_MEMBRANE_BAROSTAT, INCONSISTENT_BAROSTAT, BAROSTATED_NONPERIODIC, INCONSISTENT_INTEGRATOR, INCOMPATIBLE_SAMPLER_STATE, - INCOMPATIBLE_ENSEMBLE) = range(12) + INCOMPATIBLE_ENSEMBLE) = range(14) error_messages = { MULTIPLE_THERMOSTATS: "System has multiple thermostats.", @@ -294,6 +296,10 @@ class ThermodynamicsError(Exception): INCONSISTENT_THERMOSTAT: "System thermostat is inconsistent with thermodynamic state.", MULTIPLE_BAROSTATS: "System has multiple barostats.", UNSUPPORTED_BAROSTAT: "Found unsupported barostat {} in system.", + UNSUPPORTED_ANISOTROPIC_BAROSTAT: + "MonteCarloAnisotropicBarostat is only supported if the pressure along all scaled axes is the same.", + UNSUPPORTED_MEMBRANE_BAROSTAT: + "MonteCarloAnisotropicBarostat is only supported if the surface tension is zero.", NO_BAROSTAT: "System does not have a barostat specifying the pressure.", INCONSISTENT_BAROSTAT: "System barostat is inconsistent with thermodynamic state.", BAROSTATED_NONPERIODIC: "Non-periodic systems cannot have a barostat.", @@ -696,7 +702,7 @@ def barostat(self, new_barostat): # Update the internally stored standard system, and restore the old # pressure if something goes wrong (e.g. the system is not periodic). try: - self._pressure = new_barostat.getDefaultPressure() + self._pressure = self._get_barostat_pressure(new_barostat) self._unsafe_set_system(system, fix_state=False) except ThermodynamicsError: self._pressure = old_pressure @@ -1235,9 +1241,12 @@ def _initialize(self, system, temperature=None, pressure=None): # If pressure is None, we try to infer the pressure from the barostat. barostat = self._find_barostat(system) if pressure is None and barostat is not None: - self._pressure = barostat.getDefaultPressure() + self._pressure = self._get_barostat_pressure(barostat) else: self._pressure = pressure # Pressure here can also be None. + self._barostat_type = openmm.MonteCarloBarostat + if barostat is not None: + self._barostat_type = type(barostat) # If temperature is None, we infer the temperature from a thermostat. if temperature is None: @@ -1386,7 +1395,7 @@ def _standardize_system(self, system): # Here we push the barostat at the end. barostat = self._pop_barostat(system) if barostat is not None: - barostat.setDefaultPressure(self._STANDARD_PRESSURE) + self._set_barostat_pressure(barostat, self._STANDARD_PRESSURE) system.addForce(barostat) def _compute_standard_system_hash(self, standard_system): @@ -1422,7 +1431,7 @@ def _set_context_barostat(self, context, update_pressure, update_temperature): # Apply pressure and temperature to barostat. if update_pressure: self._set_barostat_pressure(barostat, self.pressure) - context.setParameter(barostat.Pressure(), self.pressure) + self._set_barostat_pressure_in_context(barostat, self.pressure, context) if update_temperature: self._set_barostat_temperature(barostat, self.temperature) # TODO remove try except when drop openmm7.0 support @@ -1544,7 +1553,7 @@ def set_temp(_integrator): # Internal-usage: barostat handling # ------------------------------------------------------------------------- - _SUPPORTED_BAROSTATS = {'MonteCarloBarostat'} + _SUPPORTED_BAROSTATS = {'MonteCarloBarostat', 'MonteCarloAnisotropicBarostat', 'MonteCarloMembraneBarostat'} @classmethod def _find_barostat(cls, system, get_index=False): @@ -1573,6 +1582,19 @@ def _find_barostat(cls, system, get_index=False): if barostat.__class__.__name__ not in cls._SUPPORTED_BAROSTATS: raise ThermodynamicsError(ThermodynamicsError.UNSUPPORTED_BAROSTAT, barostat.__class__.__name__) + elif isinstance(barostat, openmm.MonteCarloAnisotropicBarostat): + # support only if pressure in all scaled directions is equal + pressures = barostat.getDefaultPressure().value_in_unit(unit.bar) + scaled = [barostat.getScaleX(), barostat.getScaleY(), barostat.getScaleY()] + if sum(scaled) == 0: + raise ThermodynamicsError(ThermodynamicsError.UNSUPPORTED_ANISOTROPIC_BAROSTAT) + active_pressures = [pressure for pressure, active in zip(pressures, scaled)] + if any(abs(pressure - active_pressures[0]) > 0 for pressure in active_pressures): + raise ThermodynamicsError(ThermodynamicsError.UNSUPPORTED_ANISOTROPIC_BAROSTAT) + elif isinstance(barostat, openmm.MonteCarloMembraneBarostat): + surface_tension = barostat.getDefaultSurfaceTension() + if abs(surface_tension.value_in_unit(unit.bar*unit.nanometer)) > 0: + raise ThermodynamicsError(ThermodynamicsError.UNSUPPORTED_MEMBRANE_BAROSTAT) if get_index: return force_idx, barostat return barostat @@ -1596,13 +1618,14 @@ def _pop_barostat(cls, system): return None def _is_barostat_consistent(self, barostat): - """Check the barostat's temperature and pressure.""" + """Check the barostat's type, temperature and pressure.""" try: barostat_temperature = barostat.getDefaultTemperature() except AttributeError: # versions previous to OpenMM 7.1 barostat_temperature = barostat.getTemperature() - barostat_pressure = barostat.getDefaultPressure() - is_consistent = utils.is_quantity_close(barostat_temperature, self.temperature) + barostat_pressure = self._get_barostat_pressure(barostat) + is_consistent = (type(barostat) is self._barostat_type) + is_consistent = is_consistent and utils.is_quantity_close(barostat_temperature, self.temperature) is_consistent = is_consistent and utils.is_quantity_close(barostat_pressure, self.pressure) return is_consistent @@ -1644,7 +1667,31 @@ def _set_system_pressure(self, system, pressure): @staticmethod def _set_barostat_pressure(barostat, pressure): """Set barostat pressure.""" - barostat.setDefaultPressure(pressure) + if isinstance(pressure, unit.Quantity): + pressure = pressure.value_in_unit(unit.bar) + if isinstance(barostat, openmm.MonteCarloAnisotropicBarostat): + barostat.setDefaultPressure(openmm.Vec3(pressure, pressure, pressure)*unit.bar) + else: + barostat.setDefaultPressure(pressure*unit.bar) + + @staticmethod + def _set_barostat_pressure_in_context(barostat, pressure, context): + """Set barostat pressure.""" + if isinstance(barostat, openmm.MonteCarloAnisotropicBarostat): + p = pressure.value_in_unit(unit.bar) + context.setParameter(barostat.Pressure(), openmm.Vec3(p, p, p)*unit.bar) + else: + context.setParameter(barostat.Pressure(), pressure) + + @staticmethod + def _get_barostat_pressure(barostat): + """Set barostat pressure.""" + if isinstance(barostat, openmm.MonteCarloAnisotropicBarostat): + scaled = [barostat.getScaleX(), barostat.getScaleY(), barostat.getScaleZ()] + first_scaled_axis = scaled.index(True) + return barostat.getDefaultPressure()[first_scaled_axis] + else: + return barostat.getDefaultPressure() @staticmethod def _set_barostat_temperature(barostat, temperature): diff --git a/openmmtools/tests/test_states.py b/openmmtools/tests/test_states.py index 008a94a7d..18ba4d02c 100644 --- a/openmmtools/tests/test_states.py +++ b/openmmtools/tests/test_states.py @@ -94,13 +94,36 @@ def setup_class(cls): cls.multiple_barostat_alanine.addForce(barostat) # A system with an unsupported MonteCarloAnisotropicBarostat - cls.unsupported_barostat_alanine = copy.deepcopy(cls.alanine_explicit) + cls.unsupported_anisotropic_barostat_alanine = copy.deepcopy(cls.alanine_explicit) + pressure_in_bars = cls.std_pressure / unit.bar + anisotropic_pressure = openmm.Vec3(pressure_in_bars, pressure_in_bars, + 1.0+pressure_in_bars) + cls.unsupported_anisotropic_barostat = openmm.MonteCarloAnisotropicBarostat(anisotropic_pressure, cls.std_temperature) + cls.unsupported_anisotropic_barostat_alanine.addForce(cls.unsupported_anisotropic_barostat) + + # A system with an unsupported MonteCarloMembraneBarostat + cls.unsupported_membrane_barostat_alanine = copy.deepcopy(cls.alanine_explicit) + cls.unsupported_membrane_barostat = openmm.MonteCarloMembraneBarostat( + cls.std_pressure, 1.0 * unit.bar * unit.nanometer, cls.std_temperature, + openmm.MonteCarloMembraneBarostat.XYIsotropic, openmm.MonteCarloMembraneBarostat.ZFree + ) + cls.unsupported_membrane_barostat_alanine.addForce(cls.unsupported_membrane_barostat) + + # A system with a supported MonteCarloAnisotropicBarostat + cls.supported_anisotropic_barostat_alanine = copy.deepcopy(cls.alanine_explicit) pressure_in_bars = cls.std_pressure / unit.bar anisotropic_pressure = openmm.Vec3(pressure_in_bars, pressure_in_bars, pressure_in_bars) - cls.anisotropic_barostat = openmm.MonteCarloAnisotropicBarostat(anisotropic_pressure, - cls.std_temperature) - cls.unsupported_barostat_alanine.addForce(cls.anisotropic_barostat) + cls.supported_anisotropic_barostat = openmm.MonteCarloAnisotropicBarostat(anisotropic_pressure, cls.std_temperature) + cls.supported_anisotropic_barostat_alanine.addForce(cls.supported_anisotropic_barostat) + + # A system with a supported MonteCarloMembraneBarostat + cls.supported_membrane_barostat_alanine = copy.deepcopy(cls.alanine_explicit) + cls.supported_membrane_barostat = openmm.MonteCarloMembraneBarostat( + cls.std_pressure, 0.0 * unit.bar * unit.nanometer, cls.std_temperature, + openmm.MonteCarloMembraneBarostat.XYIsotropic, openmm.MonteCarloMembraneBarostat.ZFree + ) + cls.supported_membrane_barostat_alanine.addForce(cls.supported_membrane_barostat) # A system with an inconsistent pressure in the barostat. cls.inconsistent_pressure_alanine = copy.deepcopy(cls.alanine_explicit) @@ -166,10 +189,17 @@ def test_method_find_barostat(self): barostat = ThermodynamicState._find_barostat(self.barostated_alanine) assert isinstance(barostat, openmm.MonteCarloBarostat) + barostat = ThermodynamicState._find_barostat(self.supported_anisotropic_barostat_alanine) + assert isinstance(barostat, openmm.MonteCarloAnisotropicBarostat) + + barostat = ThermodynamicState._find_barostat(self.supported_membrane_barostat_alanine) + assert isinstance(barostat, openmm.MonteCarloMembraneBarostat) + # Raise exception if multiple or unsupported barostats found TE = ThermodynamicsError # shortcut test_cases = [(self.multiple_barostat_alanine, TE.MULTIPLE_BAROSTATS), - (self.unsupported_barostat_alanine, TE.UNSUPPORTED_BAROSTAT)] + (self.unsupported_anisotropic_barostat_alanine, TE.UNSUPPORTED_ANISOTROPIC_BAROSTAT), + (self.unsupported_membrane_barostat_alanine, TE.UNSUPPORTED_MEMBRANE_BAROSTAT)] for system, err_code in test_cases: with nose.tools.assert_raises(ThermodynamicsError) as cm: ThermodynamicState._find_barostat(system) @@ -205,6 +235,9 @@ def test_method_is_barostat_consistent(self): barostat = openmm.MonteCarloBarostat(pressure, temperature + 10*unit.kelvin) assert not state._is_barostat_consistent(barostat) + assert not state._is_barostat_consistent(self.supported_anisotropic_barostat) + assert not state._is_barostat_consistent(self.supported_membrane_barostat) + def test_method_set_system_temperature(self): """ThermodynamicState._set_system_temperature() method.""" system = copy.deepcopy(self.alanine_no_thermostat) @@ -222,29 +255,35 @@ def test_method_set_system_temperature(self): def test_property_temperature(self): """ThermodynamicState.temperature property.""" - state = ThermodynamicState(self.barostated_alanine, - self.std_temperature) - assert state.temperature == self.std_temperature - - temperature = self.std_temperature + 10.0*unit.kelvin - state.temperature = temperature - assert state.temperature == temperature - assert get_barostat_temperature(state.barostat) == temperature - - # Setting temperature to None raise error. - with nose.tools.assert_raises(ThermodynamicsError) as cm: - state.temperature = None - assert cm.exception.code == ThermodynamicsError.NONE_TEMPERATURE + for system in [self.barostated_alanine, + self.supported_anisotropic_barostat_alanine, + self.supported_membrane_barostat_alanine]: + state = ThermodynamicState(system, + self.std_temperature) + assert state.temperature == self.std_temperature + + temperature = self.std_temperature + 10.0*unit.kelvin + state.temperature = temperature + assert state.temperature == temperature + assert get_barostat_temperature(state.barostat) == temperature + + # Setting temperature to None raise error. + with nose.tools.assert_raises(ThermodynamicsError) as cm: + state.temperature = None + assert cm.exception.code == ThermodynamicsError.NONE_TEMPERATURE def test_method_set_system_pressure(self): """ThermodynamicState._set_system_pressure() method.""" - state = ThermodynamicState(self.alanine_explicit, self.std_temperature) - system = state.system - assert state._find_barostat(system) is None - state._set_system_pressure(system, self.std_pressure) - assert state._find_barostat(system).getDefaultPressure() == self.std_pressure - state._set_system_pressure(system, None) - assert state._find_barostat(system) is None + for system in [self.barostated_alanine, + self.supported_anisotropic_barostat_alanine, + self.supported_membrane_barostat_alanine]: + state = ThermodynamicState(self.alanine_explicit, self.std_temperature) + system = state.system + assert state._find_barostat(system) is None + state._set_system_pressure(system, self.std_pressure) + assert state._find_barostat(system).getDefaultPressure() == self.std_pressure + state._set_system_pressure(system, None) + assert state._find_barostat(system) is None def test_property_pressure_barostat(self): """ThermodynamicState.pressure and barostat properties.""" @@ -268,11 +307,14 @@ def test_property_pressure_barostat(self): assert state.barostat is None # Correctly reads and set system pressures - periodic_testcases = [self.alanine_explicit] + periodic_testcases = [self.alanine_explicit, + self.supported_anisotropic_barostat_alanine, + self.supported_membrane_barostat_alanine] for system in periodic_testcases: - state = ThermodynamicState(system, self.std_temperature) - assert state.pressure is None - assert state.barostat is None + if system is self.alanine_explicit: + state = ThermodynamicState(system, self.std_temperature) + assert state.pressure is None + assert state.barostat is None # Setting pressure adds a barostat state.pressure = self.std_pressure @@ -319,8 +361,23 @@ def test_property_pressure_barostat(self): # Assign incompatible barostat raise error with nose.tools.assert_raises(ThermodynamicsError) as cm: - state.barostat = self.anisotropic_barostat - assert cm.exception.code == ThermodynamicsError.UNSUPPORTED_BAROSTAT + state.barostat = self.unsupported_anisotropic_barostat + assert cm.exception.code == ThermodynamicsError.UNSUPPORTED_ANISOTROPIC_BAROSTAT + + # Assign non-standard barostat raise error + with nose.tools.assert_raises(ThermodynamicsError) as cm: + state.barostat = self.supported_anisotropic_barostat + assert cm.exception.code == ThermodynamicsError.INCONSISTENT_BAROSTAT + + # Assign incompatible barostat raise error + with nose.tools.assert_raises(ThermodynamicsError) as cm: + state.barostat = self.unsupported_membrane_barostat + assert cm.exception.code == ThermodynamicsError.UNSUPPORTED_MEMBRANE_BAROSTAT + + # Assign non-standard barostat raise error + with nose.tools.assert_raises(ThermodynamicsError) as cm: + state.barostat = self.supported_membrane_barostat + assert cm.exception.code == ThermodynamicsError.INCONSISTENT_BAROSTAT # After exception, state is left consistent assert state.pressure is None @@ -355,6 +412,10 @@ def test_property_system(self): test_cases = [(self.toluene_vacuum, TE.NO_BAROSTAT), (self.barostated_toluene, TE.BAROSTATED_NONPERIODIC), (self.multiple_barostat_alanine, TE.MULTIPLE_BAROSTATS), + (self.unsupported_anisotropic_barostat_alanine, TE.UNSUPPORTED_ANISOTROPIC_BAROSTAT), + (self.unsupported_membrane_barostat_alanine, TE.UNSUPPORTED_MEMBRANE_BAROSTAT), + (self.supported_anisotropic_barostat_alanine, TE.INCONSISTENT_BAROSTAT), + (self.supported_membrane_barostat_alanine, TE.INCONSISTENT_BAROSTAT), (self.inconsistent_pressure_alanine, TE.INCONSISTENT_BAROSTAT), (self.inconsistent_temperature_alanine, TE.INCONSISTENT_THERMOSTAT), (inconsistent_barostat_temperature, TE.INCONSISTENT_BAROSTAT)] @@ -422,7 +483,8 @@ def test_constructor_unsupported_barostat(self): TE = ThermodynamicsError # shortcut test_cases = [(self.barostated_toluene, TE.BAROSTATED_NONPERIODIC), (self.multiple_barostat_alanine, TE.MULTIPLE_BAROSTATS), - (self.unsupported_barostat_alanine, TE.UNSUPPORTED_BAROSTAT)] + (self.unsupported_anisotropic_barostat_alanine, TE.UNSUPPORTED_ANISOTROPIC_BAROSTAT), + (self.unsupported_membrane_barostat_alanine, TE.UNSUPPORTED_MEMBRANE_BAROSTAT)] for i, (system, err_code) in enumerate(test_cases): with nose.tools.assert_raises(TE) as cm: ThermodynamicState(system=system, temperature=self.std_temperature) From 37f4cf64fa0b9606269b17cb31ea1b3e0e69a449 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20Kr=C3=A4mer?= Date: Thu, 26 Sep 2019 16:00:55 -0400 Subject: [PATCH 02/12] added ThermodynamicState support for non-zero surface tension --- docs/releasehistory.rst | 5 ++ openmmtools/states.py | 107 ++++++++++++++++++++++++----- openmmtools/tests/test_states.py | 111 ++++++++++++++++++++++++------- 3 files changed, 181 insertions(+), 42 deletions(-) diff --git a/docs/releasehistory.rst b/docs/releasehistory.rst index f9d429a57..808ffa119 100644 --- a/docs/releasehistory.rst +++ b/docs/releasehistory.rst @@ -1,6 +1,11 @@ Release History *************** +Enhancements +------------ +- Added support for anisotropic and membrane barostats in `ThermodynamicState` + + 0.18.3 - Storage enhancements and bugfixes ========================================== diff --git a/openmmtools/states.py b/openmmtools/states.py index 0456aa55f..4b5b2d904 100644 --- a/openmmtools/states.py +++ b/openmmtools/states.py @@ -242,6 +242,23 @@ def _box_vectors_volume(box_vectors): return np.linalg.det(box_matrix) * a.unit**3 +def _box_vectors_area(box_vectors): + """Return the xy-area of the box vectors. + + Parameters + ---------- + box_vectors : simtk.unit.Quantity + Vectors defining the box. + + Returns + ------- + + area : simtk.unit.Quantity + The box area in units of length^2. + """ + return box_vectors[0][0] * box_vectors[1][1] + + # ============================================================================= # CUSTOM EXCEPTIONS # ============================================================================= @@ -282,7 +299,7 @@ class ThermodynamicsError(Exception): NO_BAROSTAT, UNSUPPORTED_BAROSTAT, UNSUPPORTED_ANISOTROPIC_BAROSTAT, - UNSUPPORTED_MEMBRANE_BAROSTAT, + SURFACE_TENSION_NOT_SUPPORTED, INCONSISTENT_BAROSTAT, BAROSTATED_NONPERIODIC, INCONSISTENT_INTEGRATOR, @@ -298,8 +315,8 @@ class ThermodynamicsError(Exception): UNSUPPORTED_BAROSTAT: "Found unsupported barostat {} in system.", UNSUPPORTED_ANISOTROPIC_BAROSTAT: "MonteCarloAnisotropicBarostat is only supported if the pressure along all scaled axes is the same.", - UNSUPPORTED_MEMBRANE_BAROSTAT: - "MonteCarloAnisotropicBarostat is only supported if the surface tension is zero.", + SURFACE_TENSION_NOT_SUPPORTED: + "Surface tension can only be set for states that have a system with a MonteCarloBarostat.", NO_BAROSTAT: "System does not have a barostat specifying the pressure.", INCONSISTENT_BAROSTAT: "System barostat is inconsistent with thermodynamic state.", BAROSTATED_NONPERIODIC: "Non-periodic systems cannot have a barostat.", @@ -368,11 +385,19 @@ class ThermodynamicState(object): state. It can be used to create new OpenMM Contexts, or to convert an existing Context to this particular thermodynamic state. - Only NVT and NPT ensembles are supported. The temperature must + NVT, NPT and NPgammaT ensembles are supported. The temperature must be specified in the constructor, either implicitly via a thermostat force in the system, or explicitly through the temperature parameter, which overrides an eventual thermostat indication. + To set a ThermodynamicState up in the NPgammaT ensemble, the system + passed to the constructor has to have a MonteCarloMembraneBarostat. + + To set a ThermodynamicState up with anisotropic pressure control, + the system passed to the constructor has to have a MonteCarloAnisotropicBarostat. + Currently the MonteCarloAnisotropicBarostat is only supported if + the pressure is equal for all axes that are under pressure control. + Parameters ---------- system : simtk.openmm.System @@ -394,6 +419,7 @@ class ThermodynamicState(object): system temperature pressure + surface_tension volume n_particles @@ -756,6 +782,26 @@ def is_periodic(self): """True if the system is in a periodic box (read-only).""" return self._standard_system.usesPeriodicBoundaryConditions() + @property + def surface_tension(self): + """Surface tension; None if the barostat is not a MonteCarloMembraneBarostat""" + barostat = self._find_barostat(self._standard_system) + if isinstance(barostat, openmm.MonteCarloMembraneBarostat): + return barostat.getDefaultSurfaceTension() + else: + return None + + @surface_tension.setter + def surface_tension(self, gamma): + # working around a bug in the unit conversion https://github.com/openmm/openmm/issues/2406 + if isinstance(gamma, unit.Quantity): + gamma = gamma.value_in_unit(unit.bar * unit.nanometer) + barostat = self._find_barostat(self._standard_system) + if isinstance(barostat, openmm.MonteCarloMembraneBarostat): + barostat.setDefaultSurfaceTension(gamma) + else: + raise ThermodynamicsError(ThermodynamicsError.SURFACE_TENSION_NOT_SUPPORTED) + def reduced_potential(self, context_state): """Reduced potential in this thermodynamic state. @@ -772,15 +818,17 @@ def reduced_potential(self, context_state): Notes ----- - The reduced potential is defined as in Ref. [1] + The reduced potential is defined as in Ref. [1], + with a additional term for the surface tension - u = \beta [U(x) + p V(x) + \mu N(x)] + u = \beta [U(x) + p V(x) + \mu N(x) - \gamma A] where the thermodynamic parameters are \beta = 1/(kB T) is the inverse temperature p is the pressure \mu is the chemical potential + \gamma is the surface tension and the configurational properties are @@ -789,6 +837,7 @@ def reduced_potential(self, context_state): V(x) is the instantaneous box volume N(x) the numbers of various particle species (e.g. protons of titratible groups) + A(x) is the xy-area of the box. References ---------- @@ -833,17 +882,18 @@ def reduced_potential(self, context_state): openmm_state = context_state.getState(getEnergy=True) potential_energy = openmm_state.getPotentialEnergy() volume = openmm_state.getPeriodicBoxVolume() + area = _box_vectors_area(openmm_state.getPeriodicBoxVectors()) else: n_particles = context_state.n_particles potential_energy = context_state.potential_energy volume = context_state.volume + area = context_state.area # Check compatibility. if n_particles != self.n_particles: raise ThermodynamicsError(ThermodynamicsError.INCOMPATIBLE_SAMPLER_STATE) - return self._compute_reduced_potential(potential_energy, self.temperature, - volume, self.pressure) + volume, self.pressure, area, self.surface_tension) @classmethod def reduced_potential_at_states(cls, context, thermodynamic_states): @@ -1147,7 +1197,7 @@ def apply_to_context(self, context): 310.0 """ - self._set_context_barostat(context, update_pressure=True, update_temperature=True) + self._set_context_barostat(context, update_pressure=True, update_temperature=True, update_surface_tension=True) self._set_context_thermostat(context) # ------------------------------------------------------------------------- @@ -1416,7 +1466,7 @@ def _update_standard_system(self, standard_system): # Internal-usage: context handling # ------------------------------------------------------------------------- - def _set_context_barostat(self, context, update_pressure, update_temperature): + def _set_context_barostat(self, context, update_pressure, update_temperature, update_surface_tension): """Set the barostat parameters in the Context.""" barostat = self._find_barostat(context.getSystem()) @@ -1432,6 +1482,9 @@ def _set_context_barostat(self, context, update_pressure, update_temperature): if update_pressure: self._set_barostat_pressure(barostat, self.pressure) self._set_barostat_pressure_in_context(barostat, self.pressure, context) + if self.surface_tension is not None and update_surface_tension: + self._set_barostat_surface_tension_in_context(barostat, self.surface_tension, context) + if update_temperature: self._set_barostat_temperature(barostat, self.temperature) # TODO remove try except when drop openmm7.0 support @@ -1470,9 +1523,10 @@ def _apply_to_context_in_state(self, context, thermodynamic_state): """ update_pressure = self.pressure != thermodynamic_state.pressure update_temperature = self.temperature != thermodynamic_state.temperature + update_surface_tension = self.surface_tension != thermodynamic_state.surface_tension if update_pressure or update_temperature: - self._set_context_barostat(context, update_pressure, update_temperature) + self._set_context_barostat(context, update_pressure, update_temperature, update_surface_tension) if update_temperature: self._set_context_thermostat(context) @@ -1588,13 +1642,13 @@ def _find_barostat(cls, system, get_index=False): scaled = [barostat.getScaleX(), barostat.getScaleY(), barostat.getScaleY()] if sum(scaled) == 0: raise ThermodynamicsError(ThermodynamicsError.UNSUPPORTED_ANISOTROPIC_BAROSTAT) - active_pressures = [pressure for pressure, active in zip(pressures, scaled)] + active_pressures = [pressure for pressure, active in zip(pressures, scaled) if active] if any(abs(pressure - active_pressures[0]) > 0 for pressure in active_pressures): raise ThermodynamicsError(ThermodynamicsError.UNSUPPORTED_ANISOTROPIC_BAROSTAT) - elif isinstance(barostat, openmm.MonteCarloMembraneBarostat): - surface_tension = barostat.getDefaultSurfaceTension() - if abs(surface_tension.value_in_unit(unit.bar*unit.nanometer)) > 0: - raise ThermodynamicsError(ThermodynamicsError.UNSUPPORTED_MEMBRANE_BAROSTAT) + #elif isinstance(barostat, openmm.MonteCarloMembraneBarostat): + # surface_tension = barostat.getDefaultSurfaceTension() + # if abs(surface_tension.value_in_unit(unit.bar*unit.nanometer)) > 0: + # raise ThermodynamicsError(ThermodynamicsError.UNSUPPORTED_MEMBRANE_BAROSTAT) if get_index: return force_idx, barostat return barostat @@ -1698,6 +1752,18 @@ def _set_barostat_temperature(barostat, temperature): """Set barostat temperature.""" barostat.setDefaultTemperature(temperature) + @staticmethod + def _set_barostat_surface_tension_in_context(barostat, surface_tension, context): + """Set barostat surface tension.""" + # work around a unit conversion issue in openmm + if isinstance(surface_tension, unit.Quantity): + surface_tension = surface_tension.value_in_unit(unit.nanometer*unit.bar) + try: + context.getParameter(barostat.SurfaceTension()) + except Exception as e: + raise ThermodynamicsError(ThermodynamicsError.INCOMPATIBLE_ENSEMBLE) + context.setParameter(barostat.SurfaceTension(), surface_tension) + # ------------------------------------------------------------------------- # Internal-usage: thermostat handling # ------------------------------------------------------------------------- @@ -1762,12 +1828,14 @@ def _set_system_temperature(cls, system, temperature): # ------------------------------------------------------------------------- @staticmethod - def _compute_reduced_potential(potential_energy, temperature, volume, pressure): + def _compute_reduced_potential(potential_energy, temperature, volume, pressure, area=None, surface_tension=None): """Convert potential energy into reduced potential.""" beta = 1.0 / (unit.BOLTZMANN_CONSTANT_kB * temperature) reduced_potential = potential_energy / unit.AVOGADRO_CONSTANT_NA if pressure is not None: reduced_potential += pressure * volume + if area is not None and surface_tension is not None: + reduced_potential -= surface_tension * area return beta * reduced_potential def _find_force_groups_to_update(self, context, thermodynamic_state, memo): @@ -2032,6 +2100,11 @@ def volume(self): """The volume of the box (read-only)""" return _box_vectors_volume(self.box_vectors) + @property + def area(self): + """The xy-area of the box (read-only)""" + return _box_vectors_area(self.box_vectors) + @property def n_particles(self): """Number of particles (read-only).""" diff --git a/openmmtools/tests/test_states.py b/openmmtools/tests/test_states.py index 18ba4d02c..f78a18d68 100644 --- a/openmmtools/tests/test_states.py +++ b/openmmtools/tests/test_states.py @@ -59,6 +59,7 @@ def setup_class(cls): """Create the test systems used in the test suite.""" cls.std_pressure = ThermodynamicState._STANDARD_PRESSURE cls.std_temperature = ThermodynamicState._STANDARD_TEMPERATURE + cls.std_surface_tension = 100.0 * unit.bar * unit.nanometer alanine_explicit = testsystems.AlanineDipeptideExplicit() cls.alanine_positions = alanine_explicit.positions @@ -102,12 +103,13 @@ def setup_class(cls): cls.unsupported_anisotropic_barostat_alanine.addForce(cls.unsupported_anisotropic_barostat) # A system with an unsupported MonteCarloMembraneBarostat - cls.unsupported_membrane_barostat_alanine = copy.deepcopy(cls.alanine_explicit) - cls.unsupported_membrane_barostat = openmm.MonteCarloMembraneBarostat( - cls.std_pressure, 1.0 * unit.bar * unit.nanometer, cls.std_temperature, + cls.membrane_barostat_alanine_gamma_nonzero = copy.deepcopy(cls.alanine_explicit) + # working around a bug in the unit conversion https://github.com/openmm/openmm/issues/2406 + cls.membrane_barostat_gamma_nonzero = openmm.MonteCarloMembraneBarostat( + cls.std_pressure, cls.std_surface_tension.value_in_unit(unit.bar*unit.nanometer), cls.std_temperature, openmm.MonteCarloMembraneBarostat.XYIsotropic, openmm.MonteCarloMembraneBarostat.ZFree ) - cls.unsupported_membrane_barostat_alanine.addForce(cls.unsupported_membrane_barostat) + cls.membrane_barostat_alanine_gamma_nonzero.addForce(cls.membrane_barostat_gamma_nonzero) # A system with a supported MonteCarloAnisotropicBarostat cls.supported_anisotropic_barostat_alanine = copy.deepcopy(cls.alanine_explicit) @@ -118,12 +120,12 @@ def setup_class(cls): cls.supported_anisotropic_barostat_alanine.addForce(cls.supported_anisotropic_barostat) # A system with a supported MonteCarloMembraneBarostat - cls.supported_membrane_barostat_alanine = copy.deepcopy(cls.alanine_explicit) - cls.supported_membrane_barostat = openmm.MonteCarloMembraneBarostat( - cls.std_pressure, 0.0 * unit.bar * unit.nanometer, cls.std_temperature, + cls.membrane_barostat_alanine_gamma_zero = copy.deepcopy(cls.alanine_explicit) + cls.membrane_barostat_gamma_zero = openmm.MonteCarloMembraneBarostat( + cls.std_pressure, 0.0, cls.std_temperature, openmm.MonteCarloMembraneBarostat.XYIsotropic, openmm.MonteCarloMembraneBarostat.ZFree ) - cls.supported_membrane_barostat_alanine.addForce(cls.supported_membrane_barostat) + cls.membrane_barostat_alanine_gamma_zero.addForce(cls.membrane_barostat_gamma_zero) # A system with an inconsistent pressure in the barostat. cls.inconsistent_pressure_alanine = copy.deepcopy(cls.alanine_explicit) @@ -192,14 +194,13 @@ def test_method_find_barostat(self): barostat = ThermodynamicState._find_barostat(self.supported_anisotropic_barostat_alanine) assert isinstance(barostat, openmm.MonteCarloAnisotropicBarostat) - barostat = ThermodynamicState._find_barostat(self.supported_membrane_barostat_alanine) + barostat = ThermodynamicState._find_barostat(self.membrane_barostat_alanine_gamma_zero) assert isinstance(barostat, openmm.MonteCarloMembraneBarostat) # Raise exception if multiple or unsupported barostats found TE = ThermodynamicsError # shortcut test_cases = [(self.multiple_barostat_alanine, TE.MULTIPLE_BAROSTATS), - (self.unsupported_anisotropic_barostat_alanine, TE.UNSUPPORTED_ANISOTROPIC_BAROSTAT), - (self.unsupported_membrane_barostat_alanine, TE.UNSUPPORTED_MEMBRANE_BAROSTAT)] + (self.unsupported_anisotropic_barostat_alanine, TE.UNSUPPORTED_ANISOTROPIC_BAROSTAT)] for system, err_code in test_cases: with nose.tools.assert_raises(ThermodynamicsError) as cm: ThermodynamicState._find_barostat(system) @@ -236,7 +237,7 @@ def test_method_is_barostat_consistent(self): assert not state._is_barostat_consistent(barostat) assert not state._is_barostat_consistent(self.supported_anisotropic_barostat) - assert not state._is_barostat_consistent(self.supported_membrane_barostat) + assert not state._is_barostat_consistent(self.membrane_barostat_gamma_zero) def test_method_set_system_temperature(self): """ThermodynamicState._set_system_temperature() method.""" @@ -257,7 +258,7 @@ def test_property_temperature(self): """ThermodynamicState.temperature property.""" for system in [self.barostated_alanine, self.supported_anisotropic_barostat_alanine, - self.supported_membrane_barostat_alanine]: + self.membrane_barostat_alanine_gamma_zero]: state = ThermodynamicState(system, self.std_temperature) assert state.temperature == self.std_temperature @@ -276,7 +277,7 @@ def test_method_set_system_pressure(self): """ThermodynamicState._set_system_pressure() method.""" for system in [self.barostated_alanine, self.supported_anisotropic_barostat_alanine, - self.supported_membrane_barostat_alanine]: + self.membrane_barostat_alanine_gamma_zero]: state = ThermodynamicState(self.alanine_explicit, self.std_temperature) system = state.system assert state._find_barostat(system) is None @@ -309,7 +310,7 @@ def test_property_pressure_barostat(self): # Correctly reads and set system pressures periodic_testcases = [self.alanine_explicit, self.supported_anisotropic_barostat_alanine, - self.supported_membrane_barostat_alanine] + self.membrane_barostat_alanine_gamma_zero] for system in periodic_testcases: if system is self.alanine_explicit: state = ThermodynamicState(system, self.std_temperature) @@ -369,19 +370,42 @@ def test_property_pressure_barostat(self): state.barostat = self.supported_anisotropic_barostat assert cm.exception.code == ThermodynamicsError.INCONSISTENT_BAROSTAT - # Assign incompatible barostat raise error - with nose.tools.assert_raises(ThermodynamicsError) as cm: - state.barostat = self.unsupported_membrane_barostat - assert cm.exception.code == ThermodynamicsError.UNSUPPORTED_MEMBRANE_BAROSTAT - # Assign non-standard barostat raise error with nose.tools.assert_raises(ThermodynamicsError) as cm: - state.barostat = self.supported_membrane_barostat + state.barostat = self.membrane_barostat_gamma_zero assert cm.exception.code == ThermodynamicsError.INCONSISTENT_BAROSTAT # After exception, state is left consistent assert state.pressure is None + def test_surface_tension(self): + """Test querying and setting surface tension""" + # test setting and getting surface tension for a system without barostat + state = ThermodynamicState(self.alanine_explicit, self.std_temperature) + assert state.surface_tension is None + with nose.tools.assert_raises(ThermodynamicsError) as cm: + state.surface_tension = self.std_surface_tension + assert cm.exception.code == ThermodynamicsError.SURFACE_TENSION_NOT_SUPPORTED + + # test setting and getting surface tension for a system with a non-membrane barostat + state = ThermodynamicState(self.barostated_alanine, self.std_temperature) + assert state.surface_tension is None + with nose.tools.assert_raises(ThermodynamicsError) as cm: + state.surface_tension = self.std_surface_tension + assert cm.exception.code == ThermodynamicsError.SURFACE_TENSION_NOT_SUPPORTED + + # test setting and getting surface tension + state = ThermodynamicState(self.membrane_barostat_alanine_gamma_zero, self.std_temperature) + assert utils.is_quantity_close(state.surface_tension, 0.0 * unit.bar * unit.nanometer) + state.surface_tension = self.std_surface_tension + assert utils.is_quantity_close(state.surface_tension, self.std_surface_tension) + state.surface_tension = 0.0 * unit.bar * unit.nanometer + assert utils.is_quantity_close(state.surface_tension, 0.0 * unit.bar * unit.nanometer) + + # test initial surface tension of nonzero-gamma barostat + state = ThermodynamicState(self.membrane_barostat_alanine_gamma_nonzero, self.std_temperature) + assert utils.is_quantity_close(state.surface_tension, self.std_surface_tension) + def test_property_volume(self): """Check that volume is computed correctly.""" # For volume-fluctuating systems volume is None. @@ -413,9 +437,8 @@ def test_property_system(self): (self.barostated_toluene, TE.BAROSTATED_NONPERIODIC), (self.multiple_barostat_alanine, TE.MULTIPLE_BAROSTATS), (self.unsupported_anisotropic_barostat_alanine, TE.UNSUPPORTED_ANISOTROPIC_BAROSTAT), - (self.unsupported_membrane_barostat_alanine, TE.UNSUPPORTED_MEMBRANE_BAROSTAT), (self.supported_anisotropic_barostat_alanine, TE.INCONSISTENT_BAROSTAT), - (self.supported_membrane_barostat_alanine, TE.INCONSISTENT_BAROSTAT), + (self.membrane_barostat_alanine_gamma_zero, TE.INCONSISTENT_BAROSTAT), (self.inconsistent_pressure_alanine, TE.INCONSISTENT_BAROSTAT), (self.inconsistent_temperature_alanine, TE.INCONSISTENT_THERMOSTAT), (inconsistent_barostat_temperature, TE.INCONSISTENT_BAROSTAT)] @@ -483,8 +506,8 @@ def test_constructor_unsupported_barostat(self): TE = ThermodynamicsError # shortcut test_cases = [(self.barostated_toluene, TE.BAROSTATED_NONPERIODIC), (self.multiple_barostat_alanine, TE.MULTIPLE_BAROSTATS), - (self.unsupported_anisotropic_barostat_alanine, TE.UNSUPPORTED_ANISOTROPIC_BAROSTAT), - (self.unsupported_membrane_barostat_alanine, TE.UNSUPPORTED_MEMBRANE_BAROSTAT)] + (self.unsupported_anisotropic_barostat_alanine, TE.UNSUPPORTED_ANISOTROPIC_BAROSTAT) + ] for i, (system, err_code) in enumerate(test_cases): with nose.tools.assert_raises(TE) as cm: ThermodynamicState(system=system, temperature=self.std_temperature) @@ -743,9 +766,26 @@ def test_method_apply_to_context(self): state2.apply_to_context(context) assert cm.exception.code == ThermodynamicsError.INCOMPATIBLE_ENSEMBLE + state3 = ThermodynamicState(self.membrane_barostat_alanine_gamma_zero, self.std_temperature) + with nose.tools.assert_raises(ThermodynamicsError) as cm: + state3.apply_to_context(context) + assert cm.exception.code == ThermodynamicsError.INCOMPATIBLE_ENSEMBLE + + # apply surface tension + gamma_context = openmm.Context(self.membrane_barostat_alanine_gamma_zero, openmm.VerletIntegrator(0.001)) + state3.apply_to_context(gamma_context) + assert gamma_context.getParameter(self.membrane_barostat_gamma_nonzero.SurfaceTension()) == 0.0 + state3.surface_tension = self.std_surface_tension + state3.apply_to_context(gamma_context) + assert (gamma_context.getParameter(self.membrane_barostat_gamma_nonzero.SurfaceTension()) + == self.std_surface_tension.value_in_unit(unit.nanometer*unit.bar)) + state3.surface_tension = 0.0 * unit.nanometer * unit.bar + + # Clean up contexts. del context, langevin_integrator del thermostated_context, verlet_integrator + del gamma_context verlet_integrator = openmm.VerletIntegrator(time_step) nvt_context = create_default_context(state2, verlet_integrator) @@ -754,6 +794,9 @@ def test_method_apply_to_context(self): assert cm.exception.code == ThermodynamicsError.INCOMPATIBLE_ENSEMBLE del nvt_context, verlet_integrator + + + def test_method_reduced_potential(self): """ThermodynamicState.reduced_potential() method.""" kj_mol = unit.kilojoule_per_mole @@ -785,6 +828,23 @@ def test_method_reduced_potential(self): state.reduced_potential(incompatible_sampler_state) assert cm.exception.code == ThermodynamicsError.INCOMPATIBLE_SAMPLER_STATE + # Compute constant surface tension reduced potential. + state = ThermodynamicState(self.membrane_barostat_alanine_gamma_nonzero, self.std_temperature) + integrator = openmm.VerletIntegrator(1.0 * unit.femtosecond) + context = create_default_context(state, integrator) + context.setPositions(self.alanine_positions) + sampler_state = SamplerState.from_context(context) + state.pressure = self.std_pressure + reduced_potential = state.reduced_potential(sampler_state) + pressure_volume_work = (self.std_pressure * sampler_state.volume * + unit.AVOGADRO_CONSTANT_NA) + surface_work = (self.std_surface_tension * sampler_state.area * + unit.AVOGADRO_CONSTANT_NA) + potential_energy = (reduced_potential / beta - pressure_volume_work + surface_work) / kj_mol + + assert np.isclose(sampler_state.potential_energy / kj_mol, potential_energy) + assert np.isclose(reduced_potential, state.reduced_potential(context)) + def test_method_reduced_potential_at_states(self): """ThermodynamicState.reduced_potential_at_states() method. @@ -1036,6 +1096,7 @@ def test_operator_getitem(self): # The other attributes are copied correctly. assert sliced_sampler_state.volume == sampler_state.volume + assert sliced_sampler_state.area == sampler_state.area # Energies are undefined for as subset of atoms. assert sliced_sampler_state.kinetic_energy is None From e75f21286082be53f66125c687c4824216b16381 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20Kr=C3=A4mer?= Date: Thu, 26 Sep 2019 16:04:52 -0400 Subject: [PATCH 03/12] removed old comment --- openmmtools/states.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/openmmtools/states.py b/openmmtools/states.py index 4b5b2d904..40ce870ce 100644 --- a/openmmtools/states.py +++ b/openmmtools/states.py @@ -1645,10 +1645,6 @@ def _find_barostat(cls, system, get_index=False): active_pressures = [pressure for pressure, active in zip(pressures, scaled) if active] if any(abs(pressure - active_pressures[0]) > 0 for pressure in active_pressures): raise ThermodynamicsError(ThermodynamicsError.UNSUPPORTED_ANISOTROPIC_BAROSTAT) - #elif isinstance(barostat, openmm.MonteCarloMembraneBarostat): - # surface_tension = barostat.getDefaultSurfaceTension() - # if abs(surface_tension.value_in_unit(unit.bar*unit.nanometer)) > 0: - # raise ThermodynamicsError(ThermodynamicsError.UNSUPPORTED_MEMBRANE_BAROSTAT) if get_index: return force_idx, barostat return barostat From a6f1f386aabd39ddea7082755893f7d439e5b046 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20Kr=C3=A4mer?= Date: Fri, 27 Sep 2019 10:43:15 -0400 Subject: [PATCH 04/12] enable platform properties in ThermodynamicState.create_context --- openmmtools/states.py | 13 +++++++++++-- openmmtools/tests/test_states.py | 21 +++++++++++++++++++++ 2 files changed, 32 insertions(+), 2 deletions(-) diff --git a/openmmtools/states.py b/openmmtools/states.py index 40ce870ce..3933b9d52 100644 --- a/openmmtools/states.py +++ b/openmmtools/states.py @@ -1074,7 +1074,7 @@ def is_context_compatible(self, context): is_compatible = self._standard_system_hash == context_system_hash return is_compatible - def create_context(self, integrator, platform=None): + def create_context(self, integrator, platform=None, platform_properties=None): """Create a context in this ThermodynamicState. The context contains a copy of the system. If the integrator @@ -1098,6 +1098,9 @@ def create_context(self, integrator, platform=None): platform : simtk.openmm.Platform, optional Platform to use. If None, OpenMM tries to select the fastest available platform. Default is None. + platform_properties : dict, optional + A dictionary of platform properties. Requires platform to be + specified. Returns ------- @@ -1109,6 +1112,8 @@ def create_context(self, integrator, platform=None): ThermodynamicsError If the integrator has a temperature different from this ThermodynamicState. + ValueError + If platform_properties is specified, but platform is None Examples -------- @@ -1149,9 +1154,13 @@ def create_context(self, integrator, platform=None): # Create context. if platform is None: + if platform_properties is not None: + raise ValueError("To set platform_properties, you need to also specify the platform.") return openmm.Context(system, integrator) - else: + elif platform_properties is None: return openmm.Context(system, integrator, platform) + else: + return openmm.Context(system, integrator, platform, platform_properties) def apply_to_context(self, context): """Apply this ThermodynamicState to the context. diff --git a/openmmtools/tests/test_states.py b/openmmtools/tests/test_states.py index f78a18d68..1c1dd91c1 100644 --- a/openmmtools/tests/test_states.py +++ b/openmmtools/tests/test_states.py @@ -669,6 +669,27 @@ def test_method_create_context(self): # Get rid of old context. This test can create a lot of them. del context, integrator + # test platform properties + state = ThermodynamicState(self.toluene_vacuum, self.std_temperature) + platform_properties = {"CpuThreads": "2"} + with nose.tools.assert_raises(ValueError) as cm: + state.create_context( + openmm.VerletIntegrator(0.001), + platform=None, + platform_properties=platform_properties + ) + assert str(cm.exception) == "To set platform_properties, you need to also specify the platform." + + platform = openmm.Platform.getPlatformByName("CPU") + context = state.create_context( + openmm.VerletIntegrator(0.001), + platform=platform, + platform_properties=platform_properties + ) + assert context.getPlatform().getPropertyValue(context, "CpuThreads") == "2" + del context + + def test_method_is_compatible(self): """ThermodynamicState context and state compatibility methods.""" From b6a60d6e6a21f681e7243686e5c6ba7a5fd33c06 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20Kr=C3=A4mer?= Date: Fri, 27 Sep 2019 11:01:17 -0400 Subject: [PATCH 05/12] enable platform properties in ContextCache --- docs/releasehistory.rst | 2 +- openmmtools/cache.py | 10 ++++++++-- openmmtools/tests/test_cache.py | 19 +++++++++++++++++++ 3 files changed, 28 insertions(+), 3 deletions(-) diff --git a/docs/releasehistory.rst b/docs/releasehistory.rst index 808ffa119..3aa6a6ba2 100644 --- a/docs/releasehistory.rst +++ b/docs/releasehistory.rst @@ -4,7 +4,7 @@ Release History Enhancements ------------ - Added support for anisotropic and membrane barostats in `ThermodynamicState` - +- Added support for platform properties in ContextCache (e.g. for mixed and double precision CUDA in multistate sampler) 0.18.3 - Storage enhancements and bugfixes ========================================== diff --git a/openmmtools/cache.py b/openmmtools/cache.py index 3bd0ffe94..12ab7449b 100644 --- a/openmmtools/cache.py +++ b/openmmtools/cache.py @@ -229,6 +229,9 @@ class ContextCache(object): platform : simtk.openmm.Platform, optional The OpenMM platform to use to create Contexts. If None, OpenMM tries to select the fastest one available (default is None). + platform_properties : dict, optional + A dictionary of platform properties for the OpenMM platform. + Only valid if the platform is not None (default is None). **kwargs Parameters to pass to the underlying LRUCache constructor such as capacity and time_to_live. @@ -303,8 +306,11 @@ class ContextCache(object): """ - def __init__(self, platform=None, **kwargs): + def __init__(self, platform=None, platform_properties=None, **kwargs): self._platform = platform + self._platform_properties = platform_properties + if platform_properties is not None and platform is None: + raise ValueError("To set platform_properties, you need to also specify the platform.") self._lru = LRUCache(**kwargs) def __len__(self): @@ -429,7 +435,7 @@ def get_context(self, thermodynamic_state, integrator=None): try: context = self._lru[context_id] except KeyError: - context = thermodynamic_state.create_context(integrator, self._platform) + context = thermodynamic_state.create_context(integrator, self._platform, self._platform_properties) self._lru[context_id] = context context_integrator = context.getIntegrator() diff --git a/openmmtools/tests/test_cache.py b/openmmtools/tests/test_cache.py index bc74fa2f8..eeeeca4e0 100644 --- a/openmmtools/tests/test_cache.py +++ b/openmmtools/tests/test_cache.py @@ -357,3 +357,22 @@ def test_platform_property(self): cache.get_context(self.compatible_states[0], integrator) with nose.tools.assert_raises(RuntimeError): cache.platform = platforms[0] + + def test_platform_properties(self): + platform_properties = {"CpuThreads": "2"} + with nose.tools.assert_raises(ValueError) as cm: + ContextCache(platform=None, platform_properties=platform_properties) + assert str(cm.exception) == "To set platform_properties, you need to also specify the platform." + + platform = openmm.Platform.getPlatformByName("CPU") + cache = ContextCache( + platform=platform, + platform_properties=platform_properties + ) + + thermodynamic_state = copy.deepcopy(self.water_300k) + integrator = integrators.LangevinIntegrator(temperature=300 * unit.kelvin, measure_heat=True, + measure_shadow_work=True) + context, _ = cache.get_context(thermodynamic_state, integrator) + assert context.getPlatform().getPropertyValue(context, "CpuThreads") == "2" + del context From 03158ef0af6b971bf30ea4d0c18a66c670721113 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20Kr=C3=A4mer?= Date: Fri, 27 Sep 2019 18:07:46 -0400 Subject: [PATCH 06/12] added platform properties to ContextCache serialization --- openmmtools/cache.py | 6 +++++- openmmtools/tests/test_cache.py | 13 ++++++++++++- 2 files changed, 17 insertions(+), 2 deletions(-) diff --git a/openmmtools/cache.py b/openmmtools/cache.py index 12ab7449b..62f1f86c7 100644 --- a/openmmtools/cache.py +++ b/openmmtools/cache.py @@ -453,13 +453,17 @@ def __getstate__(self): else: platform_serialization = None return dict(platform=platform_serialization, capacity=self.capacity, - time_to_live=self.time_to_live) + time_to_live=self.time_to_live, platform_properties=self._platform_properties) def __setstate__(self, serialization): if serialization['platform'] is None: self._platform = None else: self._platform = openmm.Platform.getPlatformByName(serialization['platform']) + if not 'platform_properties' in serialization: + self._platform_properties = None + else: + self._platform_properties = serialization["platform_properties"] self._lru = LRUCache(serialization['capacity'], serialization['time_to_live']) # ------------------------------------------------------------------------- diff --git a/openmmtools/tests/test_cache.py b/openmmtools/tests/test_cache.py index eeeeca4e0..1d62248eb 100644 --- a/openmmtools/tests/test_cache.py +++ b/openmmtools/tests/test_cache.py @@ -359,20 +359,31 @@ def test_platform_property(self): cache.platform = platforms[0] def test_platform_properties(self): + # Failure test platform_properties = {"CpuThreads": "2"} with nose.tools.assert_raises(ValueError) as cm: ContextCache(platform=None, platform_properties=platform_properties) assert str(cm.exception) == "To set platform_properties, you need to also specify the platform." + # test platform = openmm.Platform.getPlatformByName("CPU") cache = ContextCache( platform=platform, platform_properties=platform_properties ) - thermodynamic_state = copy.deepcopy(self.water_300k) integrator = integrators.LangevinIntegrator(temperature=300 * unit.kelvin, measure_heat=True, measure_shadow_work=True) context, _ = cache.get_context(thermodynamic_state, integrator) assert context.getPlatform().getPropertyValue(context, "CpuThreads") == "2" + + # test serialization + cache.__setstate__(cache.__getstate__()) + assert cache._platform_properties == {"CpuThreads": "2"} + + # test serialization no 2 + state = cache.__getstate__() + state["platform_properties"] = {"CpuThreads": "3"} + cache.__setstate__(state) + assert cache._platform_properties == {"CpuThreads": "3"} del context From ba91110bf032a425ba6621b781df09184582a744 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20Kr=C3=A4mer?= Date: Mon, 21 Oct 2019 17:15:18 -0400 Subject: [PATCH 07/12] Merge branch 'master' into other_barostats update release history refactor area->area_xy # Conflicts: # docs/releasehistory.rst --- docs/releasehistory.rst | 4 ++-- openmmtools/states.py | 18 +++++++++--------- openmmtools/tests/test_states.py | 6 +++--- 3 files changed, 14 insertions(+), 14 deletions(-) diff --git a/docs/releasehistory.rst b/docs/releasehistory.rst index 2b46a0baa..e24eaf8cd 100644 --- a/docs/releasehistory.rst +++ b/docs/releasehistory.rst @@ -7,8 +7,8 @@ Release History New features ------------ - Added support in ``AbsoluteAlchemicalFactory`` for handling multiple independent alchemical regions (`#438 `_). -- Added support for anisotropic and membrane barostats in `ThermodynamicState` -- Added support for platform properties in ContextCache (e.g. for mixed and double precision CUDA in multistate sampler) +- Added support for anisotropic and membrane barostats in `ThermodynamicState` (`#437 `_) +- Added support for platform properties in ContextCache (e.g. for mixed and double precision CUDA in multistate sampler) (`#437 `_) 0.18.3 - Storage enhancements and bugfixes ========================================== diff --git a/openmmtools/states.py b/openmmtools/states.py index 3933b9d52..55906fa2c 100644 --- a/openmmtools/states.py +++ b/openmmtools/states.py @@ -242,7 +242,7 @@ def _box_vectors_volume(box_vectors): return np.linalg.det(box_matrix) * a.unit**3 -def _box_vectors_area(box_vectors): +def _box_vectors_area_xy(box_vectors): """Return the xy-area of the box vectors. Parameters @@ -253,7 +253,7 @@ def _box_vectors_area(box_vectors): Returns ------- - area : simtk.unit.Quantity + area_xy : simtk.unit.Quantity The box area in units of length^2. """ return box_vectors[0][0] * box_vectors[1][1] @@ -882,12 +882,12 @@ def reduced_potential(self, context_state): openmm_state = context_state.getState(getEnergy=True) potential_energy = openmm_state.getPotentialEnergy() volume = openmm_state.getPeriodicBoxVolume() - area = _box_vectors_area(openmm_state.getPeriodicBoxVectors()) + area = _box_vectors_area_xy(openmm_state.getPeriodicBoxVectors()) else: n_particles = context_state.n_particles potential_energy = context_state.potential_energy volume = context_state.volume - area = context_state.area + area = context_state.area_xy # Check compatibility. if n_particles != self.n_particles: @@ -1833,14 +1833,14 @@ def _set_system_temperature(cls, system, temperature): # ------------------------------------------------------------------------- @staticmethod - def _compute_reduced_potential(potential_energy, temperature, volume, pressure, area=None, surface_tension=None): + def _compute_reduced_potential(potential_energy, temperature, volume, pressure, area_xy=None, surface_tension=None): """Convert potential energy into reduced potential.""" beta = 1.0 / (unit.BOLTZMANN_CONSTANT_kB * temperature) reduced_potential = potential_energy / unit.AVOGADRO_CONSTANT_NA if pressure is not None: reduced_potential += pressure * volume - if area is not None and surface_tension is not None: - reduced_potential -= surface_tension * area + if area_xy is not None and surface_tension is not None: + reduced_potential -= surface_tension * area_xy return beta * reduced_potential def _find_force_groups_to_update(self, context, thermodynamic_state, memo): @@ -2106,9 +2106,9 @@ def volume(self): return _box_vectors_volume(self.box_vectors) @property - def area(self): + def area_xy(self): """The xy-area of the box (read-only)""" - return _box_vectors_area(self.box_vectors) + return _box_vectors_area_xy(self.box_vectors) @property def n_particles(self): diff --git a/openmmtools/tests/test_states.py b/openmmtools/tests/test_states.py index 1c1dd91c1..083044714 100644 --- a/openmmtools/tests/test_states.py +++ b/openmmtools/tests/test_states.py @@ -859,8 +859,8 @@ def test_method_reduced_potential(self): reduced_potential = state.reduced_potential(sampler_state) pressure_volume_work = (self.std_pressure * sampler_state.volume * unit.AVOGADRO_CONSTANT_NA) - surface_work = (self.std_surface_tension * sampler_state.area * - unit.AVOGADRO_CONSTANT_NA) + surface_work = (self.std_surface_tension * sampler_state.area_xy * + unit.AVOGADRO_CONSTANT_NA) potential_energy = (reduced_potential / beta - pressure_volume_work + surface_work) / kj_mol assert np.isclose(sampler_state.potential_energy / kj_mol, potential_energy) @@ -1117,7 +1117,7 @@ def test_operator_getitem(self): # The other attributes are copied correctly. assert sliced_sampler_state.volume == sampler_state.volume - assert sliced_sampler_state.area == sampler_state.area + assert sliced_sampler_state.area_xy == sampler_state.area_xy # Energies are undefined for as subset of atoms. assert sliced_sampler_state.kinetic_energy is None From e4b6b06383281aadfe638907fd00d90ba4dfec71 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20Kr=C3=A4mer?= Date: Tue, 22 Oct 2019 16:17:44 -0400 Subject: [PATCH 08/12] ContextCache refactoring --- openmmtools/cache.py | 44 +++++++++++++++++++++++++++++-- openmmtools/tests/test_cache.py | 46 +++++++++++++++++++++++++++++---- 2 files changed, 83 insertions(+), 7 deletions(-) diff --git a/openmmtools/cache.py b/openmmtools/cache.py index 62f1f86c7..619e45530 100644 --- a/openmmtools/cache.py +++ b/openmmtools/cache.py @@ -307,10 +307,9 @@ class ContextCache(object): """ def __init__(self, platform=None, platform_properties=None, **kwargs): + self._validate_platform_properties(platform, platform_properties) self._platform = platform self._platform_properties = platform_properties - if platform_properties is not None and platform is None: - raise ValueError("To set platform_properties, you need to also specify the platform.") self._lru = LRUCache(**kwargs) def __len__(self): @@ -330,7 +329,17 @@ def platform(self): def platform(self, new_platform): if len(self._lru) > 0: raise RuntimeError('Cannot change platform of a non-empty ContextCache') + if new_platform is None: + self._platform_properties = None + self._validate_platform_properties(new_platform, self._platform_properties) + self._platform = new_platform + + def set_platform(self, new_platform, platform_properties=None): + if len(self._lru) > 0: + raise RuntimeError('Cannot change platform of a non-empty ContextCache') + self._validate_platform_properties(new_platform, platform_properties) self._platform = new_platform + self._platform_properties = platform_properties @property def capacity(self): @@ -448,6 +457,7 @@ def get_context(self, thermodynamic_state, integrator=None): return context, context_integrator def __getstate__(self): + # this serialization format was introduced in openmmtools > 0.18.3 (pull request #437) if self.platform is not None: platform_serialization = self.platform.getName() else: @@ -456,6 +466,7 @@ def __getstate__(self): time_to_live=self.time_to_live, platform_properties=self._platform_properties) def __setstate__(self, serialization): + # this serialization format was introduced in openmmtools > 0.18.3 (pull request #437) if serialization['platform'] is None: self._platform = None else: @@ -640,6 +651,35 @@ def _default_integrator_id(cls): return cls._cached_default_integrator_id _cached_default_integrator_id = None + @staticmethod + def _validate_platform_properties(platform=None, platform_properties=None): + """Check if platform properties are valid for the platform; else raise ValueError.""" + if platform_properties is None: + return True + if platform_properties is not None and platform is None: + raise ValueError("To set platform_properties, you need to also specify the platform.") + if not isinstance(platform_properties, dict): + raise ValueError("platform_properties must be a dictionary") + for key, value in platform_properties.items(): + if not isinstance(value, str): + raise ValueError( + "All platform properties must be strings. You supplied {}: {} of type {}".format( + key, value, type(value) + ) + ) + # create a context to check if all properties are + dummy_system = openmm.System() + dummy_system.addParticle(1) + dummy_integrator = openmm.VerletIntegrator(1.0*unit.femtoseconds) + try: + openmm.Context(dummy_system, dummy_integrator, platform, platform_properties) + return True + except Exception as e: + if "Illegal property name" in str(e): + raise ValueError("Invalid platform property for this platform. {}".format(e)) + else: + raise e + # ============================================================================= # DUMMY CONTEXT CACHE diff --git a/openmmtools/tests/test_cache.py b/openmmtools/tests/test_cache.py index 1d62248eb..7fbe4ce3b 100644 --- a/openmmtools/tests/test_cache.py +++ b/openmmtools/tests/test_cache.py @@ -359,16 +359,52 @@ def test_platform_property(self): cache.platform = platforms[0] def test_platform_properties(self): - # Failure test + # Failure tests + # no platform specified platform_properties = {"CpuThreads": "2"} with nose.tools.assert_raises(ValueError) as cm: ContextCache(platform=None, platform_properties=platform_properties) - assert str(cm.exception) == "To set platform_properties, you need to also specify the platform." + # non-string value in properties + cpu_platform = openmm.Platform.getPlatformByName("CPU") + ref_platform = openmm.Platform.getPlatformByName("Reference") + with nose.tools.assert_raises(ValueError) as cm: + ContextCache(platform=cpu_platform, platform_properties={"CpuThreads": 2}) + assert "All platform properties must be strings." in str(cm.exception) + # non-dict properties + with nose.tools.assert_raises(ValueError) as cm: + ContextCache(platform=cpu_platform, platform_properties="jambalaya") + assert str(cm.exception) == "platform_properties must be a dictionary" + # invalid property + with nose.tools.assert_raises(ValueError) as cm: + ContextCache(platform=cpu_platform, platform_properties={"jambalaya": "2"}) + assert "Invalid platform property for this platform." in str(cm.exception) + + # setter + cache = ContextCache( + platform=cpu_platform, + platform_properties=platform_properties + ) + with nose.tools.assert_raises(ValueError) as cm: + cache.platform = ref_platform + assert "Invalid platform property for this platform." in str(cm.exception) + # this should work + cache.set_platform(ref_platform) + assert cache.platform == ref_platform + # assert errors are checked in set_platform + with nose.tools.assert_raises(ValueError) as cm: + cache.set_platform(cpu_platform, platform_properties={"jambalaya": "2"}) + assert "Invalid platform property for this platform." in str(cm.exception) + # assert that resetting the platform resets the properties + cache = ContextCache( + platform=cpu_platform, + platform_properties=platform_properties + ) + cache.platform = None + assert cache._platform_properties is None - # test - platform = openmm.Platform.getPlatformByName("CPU") + # Functionality test cache = ContextCache( - platform=platform, + platform=cpu_platform, platform_properties=platform_properties ) thermodynamic_state = copy.deepcopy(self.water_300k) From fb94572263a1ae5897025ea99aea65c7a2e15be5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20Kr=C3=A4mer?= Date: Wed, 23 Oct 2019 10:30:24 -0400 Subject: [PATCH 09/12] ThermodynamicState refactoring for surface tension --- openmmtools/states.py | 130 +++++++++++++++++++++++-------- openmmtools/tests/test_states.py | 28 ++++--- 2 files changed, 114 insertions(+), 44 deletions(-) diff --git a/openmmtools/states.py b/openmmtools/states.py index 55906fa2c..5125ea8fb 100644 --- a/openmmtools/states.py +++ b/openmmtools/states.py @@ -316,7 +316,7 @@ class ThermodynamicsError(Exception): UNSUPPORTED_ANISOTROPIC_BAROSTAT: "MonteCarloAnisotropicBarostat is only supported if the pressure along all scaled axes is the same.", SURFACE_TENSION_NOT_SUPPORTED: - "Surface tension can only be set for states that have a system with a MonteCarloBarostat.", + "Surface tension can only be set for states that have a system with a MonteCarloMembraneBarostat.", NO_BAROSTAT: "System does not have a barostat specifying the pressure.", INCONSISTENT_BAROSTAT: "System barostat is inconsistent with thermodynamic state.", BAROSTATED_NONPERIODIC: "Non-periodic systems cannot have a barostat.", @@ -413,6 +413,11 @@ class ThermodynamicState(object): or just set to this pressure in case it already exists. If None, the pressure is inferred from the system barostat, and NVT ensemble is assumed if there is no barostat. + surface_tension : simtk.unit.Quantity, optional + The surface tension for the system at constant surface tension. + If this is specified, the system must have a MonteCarloMembraneBarostat. + If None, the surface_tension is inferred from the barostat and + NPT/NVT ensemble is assumed if there is no MonteCarloMembraneBarostat. Attributes ---------- @@ -493,8 +498,8 @@ class ThermodynamicState(object): # Public interface # ------------------------------------------------------------------------- - def __init__(self, system, temperature=None, pressure=None): - self._initialize(system, temperature, pressure) + def __init__(self, system, temperature=None, pressure=None, surface_tension=None): + self._initialize(system, temperature, pressure, surface_tension) @property def system(self): @@ -632,6 +637,7 @@ def get_system(self, remove_thermostat=False, remove_barostat=False): self._pop_barostat(system) else: # Set pressure of standard barostat. self._set_system_pressure(system, self.pressure) + self._set_system_surface_tension(system, self.surface_tension) # Set temperature of standard thermostat and barostat. if not (remove_barostat and remove_thermostat): @@ -709,6 +715,8 @@ def barostat(self): if barostat is not None: # NPT ensemble. self._set_barostat_pressure(barostat, self.pressure) self._set_barostat_temperature(barostat, self.temperature) + if self.surface_tension is not None: + self._set_barostat_surface_tension(barostat, self.surface_tension) return barostat @barostat.setter @@ -716,10 +724,15 @@ def barostat(self, new_barostat): # If None, just remove the barostat from the standard system. if new_barostat is None: self.pressure = None + self.surface_tension = None return - # Remember old pressure in case something goes wrong. + # Remember old pressure and surface tension in case something goes wrong. old_pressure = self.pressure + old_surface_tension = self.surface_tension + # make sure that the barostat type does not change + if self.barostat is not None and type(new_barostat) is not type(self.barostat): + raise ThermodynamicsError(ThermodynamicsError.INCONSISTENT_BAROSTAT) # Build the system with the new barostat. system = self.get_system(remove_barostat=True) @@ -729,9 +742,11 @@ def barostat(self, new_barostat): # pressure if something goes wrong (e.g. the system is not periodic). try: self._pressure = self._get_barostat_pressure(new_barostat) + self._surface_tension = self._get_barostat_surface_tension(new_barostat) self._unsafe_set_system(system, fix_state=False) except ThermodynamicsError: self._pressure = old_pressure + self._surface_tension = old_surface_tension raise @property @@ -784,23 +799,15 @@ def is_periodic(self): @property def surface_tension(self): - """Surface tension; None if the barostat is not a MonteCarloMembraneBarostat""" - barostat = self._find_barostat(self._standard_system) - if isinstance(barostat, openmm.MonteCarloMembraneBarostat): - return barostat.getDefaultSurfaceTension() - else: - return None + """Surface tension""" + return self._surface_tension @surface_tension.setter def surface_tension(self, gamma): - # working around a bug in the unit conversion https://github.com/openmm/openmm/issues/2406 - if isinstance(gamma, unit.Quantity): - gamma = gamma.value_in_unit(unit.bar * unit.nanometer) - barostat = self._find_barostat(self._standard_system) - if isinstance(barostat, openmm.MonteCarloMembraneBarostat): - barostat.setDefaultSurfaceTension(gamma) - else: + if (self._surface_tension is None) != (gamma is None): raise ThermodynamicsError(ThermodynamicsError.SURFACE_TENSION_NOT_SUPPORTED) + else: + self._surface_tension = gamma def reduced_potential(self, context_state): """Reduced potential in this thermodynamic state. @@ -934,7 +941,9 @@ def reduced_potential_at_states(cls, context, thermodynamic_states): # In NPT, we'll need also the volume. is_npt = thermodynamic_states[0].pressure is not None + is_npgammat = thermodynamic_states[0].surface_tension is not None volume = None + area_xy = None energy_by_force_group = {force.getForceGroup(): 0.0*unit.kilocalories_per_mole for force in context.getSystem().getForces()} @@ -960,11 +969,13 @@ def reduced_potential_at_states(cls, context, thermodynamic_states): # Compute volume if this is the first time we obtain a state. if is_npt and volume is None: volume = openmm_state.getPeriodicBoxVolume() + if is_npgammat and area_xy is None: + area_xy = _box_vectors_area_xy(openmm_state.getPeriodicBoxVectors()) # Compute the new total reduced potential. potential_energy = unit.sum(list(energy_by_force_group.values())) reduced_potential = cls._compute_reduced_potential(potential_energy, state.temperature, - volume, state.pressure) + volume, state.pressure, area_xy, state.surface_tension) reduced_potentials[state_idx] = reduced_potential # Update groups to compute for next states. @@ -1258,12 +1269,13 @@ def __getstate__(self, skip_system=False): serialized_system = openmm.XmlSerializer.serialize(self._standard_system) serialized_system = zlib.compress(serialized_system.encode(self._ENCODING)) return dict(standard_system=serialized_system, temperature=self.temperature, - pressure=self.pressure) + pressure=self.pressure, surface_tension=self._surface_tension) def __setstate__(self, serialization): """Set the state from a dictionary representation.""" self._temperature = serialization['temperature'] self._pressure = serialization['pressure'] + self._surface_tension = serialization['surface_tension'] serialized_system = serialization['standard_system'] # Decompress system, if need be @@ -1292,7 +1304,7 @@ def __setstate__(self, serialization): # Internal-usage: initialization # ------------------------------------------------------------------------- - def _initialize(self, system, temperature=None, pressure=None): + def _initialize(self, system, temperature=None, pressure=None, surface_tension=None): """Initialize the thermodynamic state.""" # Avoid modifying the original system when setting temperature and pressure. system = copy.deepcopy(system) @@ -1303,9 +1315,15 @@ def _initialize(self, system, temperature=None, pressure=None): self._pressure = self._get_barostat_pressure(barostat) else: self._pressure = pressure # Pressure here can also be None. - self._barostat_type = openmm.MonteCarloBarostat - if barostat is not None: - self._barostat_type = type(barostat) + + # If surface tension is None, we try to infer the pressure from the barostat. + barostat_type = type(barostat) + if surface_tension is None and barostat_type == openmm.MonteCarloMembraneBarostat: + self._surface_tension = barostat.getDefaultSurfaceTension() + elif surface_tension is not None and barostat_type != openmm.MonteCarloMembraneBarostat: + raise ThermodynamicsError(ThermodynamicsError.INCOMPATIBLE_ENSEMBLE) + else: + self._surface_tension = surface_tension # If temperature is None, we infer the temperature from a thermostat. if temperature is None: @@ -1321,6 +1339,8 @@ def _initialize(self, system, temperature=None, pressure=None): self._set_system_temperature(system, temperature) if pressure is not None: self._set_system_pressure(system, pressure) + if surface_tension is not None: + self._set_system_surface_tension(system, surface_tension) # We can use the unsafe set_system since the system has been copied. self._unsafe_set_system(system, fix_state=False) @@ -1336,6 +1356,7 @@ def _initialize(self, system, temperature=None, pressure=None): # caused by unit conversion. _STANDARD_PRESSURE = 1.0*unit.bar _STANDARD_TEMPERATURE = 273.0*unit.kelvin + _STANDARD_SURFACE_TENSION = 0.0*unit.nanometer*unit.bar _NONPERIODIC_NONBONDED_METHODS = {openmm.NonbondedForce.NoCutoff, openmm.NonbondedForce.CutoffNonPeriodic} @@ -1351,12 +1372,13 @@ def _unsafe_set_system(self, system, fix_state): # Configure temperature and pressure. if fix_state: # We just need to add/remove the barostat according to the ensemble. - # Temperature and pressure of thermostat and barostat will be set + # Temperature, pressure, surface tension of thermostat and barostat will be set # to their standard value afterwards. self._set_system_pressure(system, self.pressure) + self._set_system_surface_tension(system, self.surface_tension) else: # If the flag is deactivated, we check that temperature - # and pressure of the system are correct. + # pressure, and surface tension of the system are correct. self._check_system_consistency(system) # Update standard system. @@ -1455,6 +1477,8 @@ def _standardize_system(self, system): barostat = self._pop_barostat(system) if barostat is not None: self._set_barostat_pressure(barostat, self._STANDARD_PRESSURE) + if isinstance(barostat, openmm.MonteCarloMembraneBarostat): + self._set_barostat_surface_tension(barostat, self._STANDARD_SURFACE_TENSION) system.addForce(barostat) def _compute_standard_system_hash(self, standard_system): @@ -1483,15 +1507,20 @@ def _set_context_barostat(self, context, update_pressure, update_temperature, up if (barostat is None) != (self._pressure is None): raise ThermodynamicsError(ThermodynamicsError.INCOMPATIBLE_ENSEMBLE) + if (type(barostat) is openmm.MonteCarloMembraneBarostat) == (self._surface_tension is None): + raise ThermodynamicsError(ThermodynamicsError.INCOMPATIBLE_ENSEMBLE) + # No need to set the barostat if we are in NVT. if self._pressure is None: return - # Apply pressure and temperature to barostat. + # Apply pressure, surface tension, and temperature to barostat. if update_pressure: self._set_barostat_pressure(barostat, self.pressure) self._set_barostat_pressure_in_context(barostat, self.pressure, context) + if self.surface_tension is not None and update_surface_tension: + self._set_barostat_surface_tension(barostat, self.surface_tension) self._set_barostat_surface_tension_in_context(barostat, self.surface_tension, context) if update_temperature: @@ -1534,7 +1563,7 @@ def _apply_to_context_in_state(self, context, thermodynamic_state): update_temperature = self.temperature != thermodynamic_state.temperature update_surface_tension = self.surface_tension != thermodynamic_state.surface_tension - if update_pressure or update_temperature: + if update_pressure or update_temperature or update_surface_tension: self._set_context_barostat(context, update_pressure, update_temperature, update_surface_tension) if update_temperature: self._set_context_thermostat(context) @@ -1618,6 +1647,11 @@ def set_temp(_integrator): _SUPPORTED_BAROSTATS = {'MonteCarloBarostat', 'MonteCarloAnisotropicBarostat', 'MonteCarloMembraneBarostat'} + @property + def _barostat_type(self): + barostat = self._find_barostat(self._standard_system) + return type(barostat) + @classmethod def _find_barostat(cls, system, get_index=False): """Return the first barostat found in the system. @@ -1677,16 +1711,20 @@ def _pop_barostat(cls, system): return None def _is_barostat_consistent(self, barostat): - """Check the barostat's type, temperature and pressure.""" + """Check the barostat's, temperature, pressure, and surface_tension.""" try: barostat_temperature = barostat.getDefaultTemperature() except AttributeError: # versions previous to OpenMM 7.1 barostat_temperature = barostat.getTemperature() barostat_pressure = self._get_barostat_pressure(barostat) - is_consistent = (type(barostat) is self._barostat_type) - is_consistent = is_consistent and utils.is_quantity_close(barostat_temperature, self.temperature) - is_consistent = is_consistent and utils.is_quantity_close(barostat_pressure, - self.pressure) + barostat_surface_tension = self._get_barostat_surface_tension(barostat) + + is_consistent = utils.is_quantity_close(barostat_temperature, self.temperature) + is_consistent = is_consistent and utils.is_quantity_close(barostat_pressure, self.pressure) + if barostat is not None and self._surface_tension is not None: + is_consistent = is_consistent and utils.is_quantity_close(barostat_surface_tension, self._surface_tension) + else: + is_consistent = is_consistent and (barostat_surface_tension == self._surface_tension) # both None return is_consistent def _set_system_pressure(self, system, pressure): @@ -1757,6 +1795,32 @@ def _set_barostat_temperature(barostat, temperature): """Set barostat temperature.""" barostat.setDefaultTemperature(temperature) + def _set_system_surface_tension(self, system, gamma): + """Set system surface tension""" + if gamma is not None and not system.usesPeriodicBoundaryConditions(): + raise ThermodynamicsError(ThermodynamicsError.BAROSTATED_NONPERIODIC) + barostat = self._find_barostat(system) + if (gamma is None) == isinstance(barostat, openmm.MonteCarloMembraneBarostat): + raise ThermodynamicsError(ThermodynamicsError.INCOMPATIBLE_ENSEMBLE) + self._set_barostat_surface_tension(barostat, gamma) + + + def _set_barostat_surface_tension(self, barostat, gamma): + # working around a bug in the unit conversion https://github.com/openmm/openmm/issues/2406 + if isinstance(gamma, unit.Quantity): + gamma = gamma.value_in_unit(unit.bar * unit.nanometer) + # barostat = self._find_barostat(system) + if isinstance(barostat, openmm.MonteCarloMembraneBarostat): + barostat.setDefaultSurfaceTension(gamma) + elif gamma is not None: + raise ThermodynamicsError(ThermodynamicsError.SURFACE_TENSION_NOT_SUPPORTED) + + def _get_barostat_surface_tension(self, barostat): + if isinstance(barostat, openmm.MonteCarloMembraneBarostat): + return barostat.getDefaultSurfaceTension() + else: + return None + @staticmethod def _set_barostat_surface_tension_in_context(barostat, surface_tension, context): """Set barostat surface tension.""" @@ -1765,7 +1829,7 @@ def _set_barostat_surface_tension_in_context(barostat, surface_tension, context) surface_tension = surface_tension.value_in_unit(unit.nanometer*unit.bar) try: context.getParameter(barostat.SurfaceTension()) - except Exception as e: + except Exception: raise ThermodynamicsError(ThermodynamicsError.INCOMPATIBLE_ENSEMBLE) context.setParameter(barostat.SurfaceTension(), surface_tension) diff --git a/openmmtools/tests/test_states.py b/openmmtools/tests/test_states.py index 083044714..d5c0753dd 100644 --- a/openmmtools/tests/test_states.py +++ b/openmmtools/tests/test_states.py @@ -236,7 +236,7 @@ def test_method_is_barostat_consistent(self): barostat = openmm.MonteCarloBarostat(pressure, temperature + 10*unit.kelvin) assert not state._is_barostat_consistent(barostat) - assert not state._is_barostat_consistent(self.supported_anisotropic_barostat) + #assert not state._is_barostat_consistent(self.supported_anisotropic_barostat) assert not state._is_barostat_consistent(self.membrane_barostat_gamma_zero) def test_method_set_system_temperature(self): @@ -365,15 +365,15 @@ def test_property_pressure_barostat(self): state.barostat = self.unsupported_anisotropic_barostat assert cm.exception.code == ThermodynamicsError.UNSUPPORTED_ANISOTROPIC_BAROSTAT - # Assign non-standard barostat raise error - with nose.tools.assert_raises(ThermodynamicsError) as cm: - state.barostat = self.supported_anisotropic_barostat - assert cm.exception.code == ThermodynamicsError.INCONSISTENT_BAROSTAT - - # Assign non-standard barostat raise error - with nose.tools.assert_raises(ThermodynamicsError) as cm: - state.barostat = self.membrane_barostat_gamma_zero - assert cm.exception.code == ThermodynamicsError.INCONSISTENT_BAROSTAT + # Assign barostat with different type raise error + if state.barostat is not None and type(state.barostat) != type(self.supported_anisotropic_barostat): + with nose.tools.assert_raises(ThermodynamicsError) as cm: + state.barostat = self.supported_anisotropic_barostat + assert cm.exception.code == ThermodynamicsError.INCONSISTENT_BAROSTA + if state.barostat is not None and type(state.barostat) != type(self.membrane_barostat_gamma_zero): + with nose.tools.assert_raises(ThermodynamicsError) as cm: + state.barostat = self.membrane_barostat_gamma_zero + assert cm.exception.code == ThermodynamicsError.INCONSISTENT_BAROSTAT # After exception, state is left consistent assert state.pressure is None @@ -437,7 +437,7 @@ def test_property_system(self): (self.barostated_toluene, TE.BAROSTATED_NONPERIODIC), (self.multiple_barostat_alanine, TE.MULTIPLE_BAROSTATS), (self.unsupported_anisotropic_barostat_alanine, TE.UNSUPPORTED_ANISOTROPIC_BAROSTAT), - (self.supported_anisotropic_barostat_alanine, TE.INCONSISTENT_BAROSTAT), + #(self.supported_anisotropic_barostat_alanine, TE.INCONSISTENT_BAROSTAT), (self.membrane_barostat_alanine_gamma_zero, TE.INCONSISTENT_BAROSTAT), (self.inconsistent_pressure_alanine, TE.INCONSISTENT_BAROSTAT), (self.inconsistent_temperature_alanine, TE.INCONSISTENT_THERMOSTAT), @@ -726,6 +726,12 @@ def check_compatibility(state1, state2, is_compatible): barostated_alanine2.pressure = barostated_alanine.pressure + 0.2*unit.bars check_compatibility(barostated_alanine, barostated_alanine2, True) + check_compatibility( + ThermodynamicState(self.membrane_barostat_alanine_gamma_zero), + ThermodynamicState(self.membrane_barostat_alanine_gamma_nonzero), + True + ) + def test_method_apply_to_context(self): """ThermodynamicState.apply_to_context() method.""" friction = 5.0/unit.picosecond From 101915a7c56dd15bac45998a6af1d6023386eba8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20Kr=C3=A4mer?= Date: Wed, 23 Oct 2019 16:30:00 -0400 Subject: [PATCH 10/12] ThermodynamicState._is_barostat_type_consistent --- openmmtools/states.py | 19 ++++++++++++++----- openmmtools/tests/test_states.py | 7 ++++++- 2 files changed, 20 insertions(+), 6 deletions(-) diff --git a/openmmtools/states.py b/openmmtools/states.py index 5125ea8fb..a67dcdabe 100644 --- a/openmmtools/states.py +++ b/openmmtools/states.py @@ -1421,9 +1421,6 @@ def _check_system_consistency(self, system): # This line raises MULTIPLE_BAROSTATS and UNSUPPORTED_BAROSTAT. barostat = self._find_barostat(system) if barostat is not None: - if not self._is_barostat_consistent(barostat): - raise TE(TE.INCONSISTENT_BAROSTAT) - # Check that barostat is not added to non-periodic system. We # cannot use System.usesPeriodicBoundaryConditions() because # in OpenMM < 7.1 that returns True when a barostat is added. @@ -1433,6 +1430,10 @@ def _check_system_consistency(self, system): nonbonded_method = force.getNonbondedMethod() if nonbonded_method in self._NONPERIODIC_NONBONDED_METHODS: raise TE(TE.BAROSTATED_NONPERIODIC) + + if not self._is_barostat_consistent(barostat): + raise TE(TE.INCONSISTENT_BAROSTAT) + elif self.pressure is not None: raise TE(TE.NO_BAROSTAT) @@ -1710,8 +1711,15 @@ def _pop_barostat(cls, system): return barostat return None + def _is_barostat_type_consistent(self, barostat): + # during initialization (standard system not set), any barostat type is OK + if not hasattr(self, "_standard_system"): + return True + system_barostat = self._find_barostat(self._standard_system) + return type(barostat) == type(system_barostat) + def _is_barostat_consistent(self, barostat): - """Check the barostat's, temperature, pressure, and surface_tension.""" + """Check the barostat's temperature, pressure, and surface_tension.""" try: barostat_temperature = barostat.getDefaultTemperature() except AttributeError: # versions previous to OpenMM 7.1 @@ -1719,7 +1727,8 @@ def _is_barostat_consistent(self, barostat): barostat_pressure = self._get_barostat_pressure(barostat) barostat_surface_tension = self._get_barostat_surface_tension(barostat) - is_consistent = utils.is_quantity_close(barostat_temperature, self.temperature) + is_consistent = self._is_barostat_type_consistent(barostat) + is_consistent = is_consistent and utils.is_quantity_close(barostat_temperature, self.temperature) is_consistent = is_consistent and utils.is_quantity_close(barostat_pressure, self.pressure) if barostat is not None and self._surface_tension is not None: is_consistent = is_consistent and utils.is_quantity_close(barostat_surface_tension, self._surface_tension) diff --git a/openmmtools/tests/test_states.py b/openmmtools/tests/test_states.py index d5c0753dd..b0a5deb92 100644 --- a/openmmtools/tests/test_states.py +++ b/openmmtools/tests/test_states.py @@ -236,7 +236,7 @@ def test_method_is_barostat_consistent(self): barostat = openmm.MonteCarloBarostat(pressure, temperature + 10*unit.kelvin) assert not state._is_barostat_consistent(barostat) - #assert not state._is_barostat_consistent(self.supported_anisotropic_barostat) + assert not state._is_barostat_consistent(self.supported_anisotropic_barostat) assert not state._is_barostat_consistent(self.membrane_barostat_gamma_zero) def test_method_set_system_temperature(self): @@ -731,6 +731,11 @@ def check_compatibility(state1, state2, is_compatible): ThermodynamicState(self.membrane_barostat_alanine_gamma_nonzero), True ) + check_compatibility( + ThermodynamicState(self.barostated_alanine), + ThermodynamicState(self.membrane_barostat_alanine_gamma_nonzero), + False + ) def test_method_apply_to_context(self): """ThermodynamicState.apply_to_context() method.""" From d1c18125e764b148dc3e2925b53dfb08956a5a37 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20Kr=C3=A4mer?= Date: Sat, 26 Oct 2019 12:02:26 -0400 Subject: [PATCH 11/12] final cleanup of other barostats in ThermodynamicState --- openmmtools/states.py | 9 +-------- openmmtools/tests/test_states.py | 2 +- 2 files changed, 2 insertions(+), 9 deletions(-) diff --git a/openmmtools/states.py b/openmmtools/states.py index a67dcdabe..df351a888 100644 --- a/openmmtools/states.py +++ b/openmmtools/states.py @@ -1316,7 +1316,7 @@ def _initialize(self, system, temperature=None, pressure=None, surface_tension=N else: self._pressure = pressure # Pressure here can also be None. - # If surface tension is None, we try to infer the pressure from the barostat. + # If surface tension is None, we try to infer the surface tension from the barostat. barostat_type = type(barostat) if surface_tension is None and barostat_type == openmm.MonteCarloMembraneBarostat: self._surface_tension = barostat.getDefaultSurfaceTension() @@ -1648,11 +1648,6 @@ def set_temp(_integrator): _SUPPORTED_BAROSTATS = {'MonteCarloBarostat', 'MonteCarloAnisotropicBarostat', 'MonteCarloMembraneBarostat'} - @property - def _barostat_type(self): - barostat = self._find_barostat(self._standard_system) - return type(barostat) - @classmethod def _find_barostat(cls, system, get_index=False): """Return the first barostat found in the system. @@ -1813,12 +1808,10 @@ def _set_system_surface_tension(self, system, gamma): raise ThermodynamicsError(ThermodynamicsError.INCOMPATIBLE_ENSEMBLE) self._set_barostat_surface_tension(barostat, gamma) - def _set_barostat_surface_tension(self, barostat, gamma): # working around a bug in the unit conversion https://github.com/openmm/openmm/issues/2406 if isinstance(gamma, unit.Quantity): gamma = gamma.value_in_unit(unit.bar * unit.nanometer) - # barostat = self._find_barostat(system) if isinstance(barostat, openmm.MonteCarloMembraneBarostat): barostat.setDefaultSurfaceTension(gamma) elif gamma is not None: diff --git a/openmmtools/tests/test_states.py b/openmmtools/tests/test_states.py index b0a5deb92..f230aa93c 100644 --- a/openmmtools/tests/test_states.py +++ b/openmmtools/tests/test_states.py @@ -437,7 +437,7 @@ def test_property_system(self): (self.barostated_toluene, TE.BAROSTATED_NONPERIODIC), (self.multiple_barostat_alanine, TE.MULTIPLE_BAROSTATS), (self.unsupported_anisotropic_barostat_alanine, TE.UNSUPPORTED_ANISOTROPIC_BAROSTAT), - #(self.supported_anisotropic_barostat_alanine, TE.INCONSISTENT_BAROSTAT), + (self.supported_anisotropic_barostat_alanine, TE.INCONSISTENT_BAROSTAT), (self.membrane_barostat_alanine_gamma_zero, TE.INCONSISTENT_BAROSTAT), (self.inconsistent_pressure_alanine, TE.INCONSISTENT_BAROSTAT), (self.inconsistent_temperature_alanine, TE.INCONSISTENT_THERMOSTAT), From 46308b03aa42a816bf229215a39cac141a140b31 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20Kr=C3=A4mer?= Date: Tue, 29 Oct 2019 07:01:32 -0400 Subject: [PATCH 12/12] attempted fix of surface tension test --- openmmtools/tests/test_states.py | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/openmmtools/tests/test_states.py b/openmmtools/tests/test_states.py index f230aa93c..937fe2001 100644 --- a/openmmtools/tests/test_states.py +++ b/openmmtools/tests/test_states.py @@ -59,7 +59,8 @@ def setup_class(cls): """Create the test systems used in the test suite.""" cls.std_pressure = ThermodynamicState._STANDARD_PRESSURE cls.std_temperature = ThermodynamicState._STANDARD_TEMPERATURE - cls.std_surface_tension = 100.0 * unit.bar * unit.nanometer + cls.std_surface_tension = ThermodynamicState._STANDARD_SURFACE_TENSION + cls.modified_surface_tension = 100.0 * unit.bar * unit.nanometer alanine_explicit = testsystems.AlanineDipeptideExplicit() cls.alanine_positions = alanine_explicit.positions @@ -106,7 +107,7 @@ def setup_class(cls): cls.membrane_barostat_alanine_gamma_nonzero = copy.deepcopy(cls.alanine_explicit) # working around a bug in the unit conversion https://github.com/openmm/openmm/issues/2406 cls.membrane_barostat_gamma_nonzero = openmm.MonteCarloMembraneBarostat( - cls.std_pressure, cls.std_surface_tension.value_in_unit(unit.bar*unit.nanometer), cls.std_temperature, + cls.std_pressure, cls.modified_surface_tension.value_in_unit(unit.bar*unit.nanometer), cls.std_temperature, openmm.MonteCarloMembraneBarostat.XYIsotropic, openmm.MonteCarloMembraneBarostat.ZFree ) cls.membrane_barostat_alanine_gamma_nonzero.addForce(cls.membrane_barostat_gamma_nonzero) @@ -396,15 +397,15 @@ def test_surface_tension(self): # test setting and getting surface tension state = ThermodynamicState(self.membrane_barostat_alanine_gamma_zero, self.std_temperature) - assert utils.is_quantity_close(state.surface_tension, 0.0 * unit.bar * unit.nanometer) - state.surface_tension = self.std_surface_tension - assert utils.is_quantity_close(state.surface_tension, self.std_surface_tension) + assert utils.is_quantity_close(state.surface_tension, 0.0 * unit.bar * unit.nanometer, rtol=0.0, atol=1e-10) + state.surface_tension = self.modified_surface_tension + assert utils.is_quantity_close(state.surface_tension, self.modified_surface_tension) state.surface_tension = 0.0 * unit.bar * unit.nanometer - assert utils.is_quantity_close(state.surface_tension, 0.0 * unit.bar * unit.nanometer) + assert utils.is_quantity_close(state.surface_tension, 0.0 * unit.bar * unit.nanometer, rtol=0.0, atol=1e-10) # test initial surface tension of nonzero-gamma barostat state = ThermodynamicState(self.membrane_barostat_alanine_gamma_nonzero, self.std_temperature) - assert utils.is_quantity_close(state.surface_tension, self.std_surface_tension) + assert utils.is_quantity_close(state.surface_tension, self.modified_surface_tension) def test_property_volume(self): """Check that volume is computed correctly.""" @@ -807,10 +808,10 @@ def test_method_apply_to_context(self): gamma_context = openmm.Context(self.membrane_barostat_alanine_gamma_zero, openmm.VerletIntegrator(0.001)) state3.apply_to_context(gamma_context) assert gamma_context.getParameter(self.membrane_barostat_gamma_nonzero.SurfaceTension()) == 0.0 - state3.surface_tension = self.std_surface_tension + state3.surface_tension = self.modified_surface_tension state3.apply_to_context(gamma_context) assert (gamma_context.getParameter(self.membrane_barostat_gamma_nonzero.SurfaceTension()) - == self.std_surface_tension.value_in_unit(unit.nanometer*unit.bar)) + == self.modified_surface_tension.value_in_unit(unit.nanometer*unit.bar)) state3.surface_tension = 0.0 * unit.nanometer * unit.bar @@ -870,7 +871,7 @@ def test_method_reduced_potential(self): reduced_potential = state.reduced_potential(sampler_state) pressure_volume_work = (self.std_pressure * sampler_state.volume * unit.AVOGADRO_CONSTANT_NA) - surface_work = (self.std_surface_tension * sampler_state.area_xy * + surface_work = (self.modified_surface_tension * sampler_state.area_xy * unit.AVOGADRO_CONSTANT_NA) potential_energy = (reduced_potential / beta - pressure_volume_work + surface_work) / kj_mol