diff --git a/docs/releasehistory.rst b/docs/releasehistory.rst index e145ed28a..6c52e1e70 100644 --- a/docs/releasehistory.rst +++ b/docs/releasehistory.rst @@ -7,6 +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` (`#437 `_) +- Added support for platform properties in ContextCache (e.g. for mixed and double precision CUDA in multistate sampler) (`#437 `_) Bugfixes -------- diff --git a/openmmtools/cache.py b/openmmtools/cache.py index 3bd0ffe94..619e45530 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,10 @@ class ContextCache(object): """ - def __init__(self, platform=None, **kwargs): + def __init__(self, platform=None, platform_properties=None, **kwargs): + self._validate_platform_properties(platform, platform_properties) self._platform = platform + self._platform_properties = platform_properties self._lru = LRUCache(**kwargs) def __len__(self): @@ -324,8 +329,18 @@ 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): """The maximum number of Context cached. @@ -429,7 +444,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() @@ -442,18 +457,24 @@ 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: 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): + # this serialization format was introduced in openmmtools > 0.18.3 (pull request #437) 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']) # ------------------------------------------------------------------------- @@ -630,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/states.py b/openmmtools/states.py index aad42e21c..df351a888 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_xy(box_vectors): + """Return the xy-area of the box vectors. + + Parameters + ---------- + box_vectors : simtk.unit.Quantity + Vectors defining the box. + + Returns + ------- + + area_xy : simtk.unit.Quantity + The box area in units of length^2. + """ + return box_vectors[0][0] * box_vectors[1][1] + + # ============================================================================= # CUSTOM EXCEPTIONS # ============================================================================= @@ -281,11 +298,13 @@ class ThermodynamicsError(Exception): MULTIPLE_BAROSTATS, NO_BAROSTAT, UNSUPPORTED_BAROSTAT, + UNSUPPORTED_ANISOTROPIC_BAROSTAT, + SURFACE_TENSION_NOT_SUPPORTED, 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 +313,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.", + SURFACE_TENSION_NOT_SUPPORTED: + "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.", @@ -362,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 @@ -382,12 +413,18 @@ 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 ---------- system temperature pressure + surface_tension volume n_particles @@ -461,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): @@ -600,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): @@ -677,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 @@ -684,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) @@ -696,10 +741,12 @@ 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._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 @@ -750,6 +797,18 @@ 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""" + return self._surface_tension + + @surface_tension.setter + def surface_tension(self, gamma): + 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. @@ -766,15 +825,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 @@ -783,6 +844,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 ---------- @@ -827,17 +889,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_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_xy # 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): @@ -878,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()} @@ -904,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. @@ -1018,7 +1085,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 @@ -1042,6 +1109,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 ------- @@ -1053,6 +1123,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 -------- @@ -1093,9 +1165,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. @@ -1141,7 +1217,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) # ------------------------------------------------------------------------- @@ -1193,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 @@ -1227,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) @@ -1235,10 +1312,19 @@ 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. + # 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() + 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: thermostat = self._find_thermostat(system) @@ -1253,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) @@ -1268,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} @@ -1283,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. @@ -1331,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. @@ -1343,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) @@ -1386,7 +1477,9 @@ 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) + 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): @@ -1407,7 +1500,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()) @@ -1415,14 +1508,22 @@ def _set_context_barostat(self, context, update_pressure, update_temperature): 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) - context.setParameter(barostat.Pressure(), 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: self._set_barostat_temperature(barostat, self.temperature) # TODO remove try except when drop openmm7.0 support @@ -1461,9 +1562,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) + 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) @@ -1544,7 +1646,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 +1675,15 @@ 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 active] + if any(abs(pressure - active_pressures[0]) > 0 for pressure in active_pressures): + raise ThermodynamicsError(ThermodynamicsError.UNSUPPORTED_ANISOTROPIC_BAROSTAT) if get_index: return force_idx, barostat return barostat @@ -1595,16 +1706,29 @@ 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 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 = barostat.getDefaultPressure() - is_consistent = utils.is_quantity_close(barostat_temperature, self.temperature) - is_consistent = is_consistent and utils.is_quantity_close(barostat_pressure, - self.pressure) + barostat_pressure = self._get_barostat_pressure(barostat) + barostat_surface_tension = self._get_barostat_surface_tension(barostat) + + 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) + else: + is_consistent = is_consistent and (barostat_surface_tension == self._surface_tension) # both None return is_consistent def _set_system_pressure(self, system, pressure): @@ -1644,13 +1768,73 @@ 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): """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) + 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.""" + # 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: + raise ThermodynamicsError(ThermodynamicsError.INCOMPATIBLE_ENSEMBLE) + context.setParameter(barostat.SurfaceTension(), surface_tension) + # ------------------------------------------------------------------------- # Internal-usage: thermostat handling # ------------------------------------------------------------------------- @@ -1715,12 +1899,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_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_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): @@ -1985,6 +2171,11 @@ def volume(self): """The volume of the box (read-only)""" return _box_vectors_volume(self.box_vectors) + @property + def area_xy(self): + """The xy-area of the box (read-only)""" + return _box_vectors_area_xy(self.box_vectors) + @property def n_particles(self): """Number of particles (read-only).""" diff --git a/openmmtools/tests/test_cache.py b/openmmtools/tests/test_cache.py index bc74fa2f8..7fbe4ce3b 100644 --- a/openmmtools/tests/test_cache.py +++ b/openmmtools/tests/test_cache.py @@ -357,3 +357,69 @@ 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): + # Failure tests + # no platform specified + platform_properties = {"CpuThreads": "2"} + with nose.tools.assert_raises(ValueError) as cm: + ContextCache(platform=None, platform_properties=platform_properties) + # 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 + + # Functionality test + cache = ContextCache( + platform=cpu_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 diff --git a/openmmtools/tests/test_states.py b/openmmtools/tests/test_states.py index 008a94a7d..937fe2001 100644 --- a/openmmtools/tests/test_states.py +++ b/openmmtools/tests/test_states.py @@ -59,6 +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 = ThermodynamicState._STANDARD_SURFACE_TENSION + cls.modified_surface_tension = 100.0 * unit.bar * unit.nanometer alanine_explicit = testsystems.AlanineDipeptideExplicit() cls.alanine_positions = alanine_explicit.positions @@ -94,13 +96,37 @@ 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.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.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) + + # 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.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.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) @@ -166,10 +192,16 @@ 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.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_barostat_alanine, TE.UNSUPPORTED_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) @@ -205,6 +237,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.membrane_barostat_gamma_zero) + def test_method_set_system_temperature(self): """ThermodynamicState._set_system_temperature() method.""" system = copy.deepcopy(self.alanine_no_thermostat) @@ -222,29 +257,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.membrane_barostat_alanine_gamma_zero]: + 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.membrane_barostat_alanine_gamma_zero]: + 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 +309,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.membrane_barostat_alanine_gamma_zero] 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,12 +363,50 @@ 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 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 + 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, 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, 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.modified_surface_tension) + def test_property_volume(self): """Check that volume is computed correctly.""" # For volume-fluctuating systems volume is None. @@ -355,6 +437,9 @@ 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.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), (inconsistent_barostat_temperature, TE.INCONSISTENT_BAROSTAT)] @@ -422,7 +507,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) + ] for i, (system, err_code) in enumerate(test_cases): with nose.tools.assert_raises(TE) as cm: ThermodynamicState(system=system, temperature=self.std_temperature) @@ -584,6 +670,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.""" @@ -620,6 +727,17 @@ 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 + ) + 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.""" friction = 5.0/unit.picosecond @@ -681,9 +799,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.modified_surface_tension + state3.apply_to_context(gamma_context) + assert (gamma_context.getParameter(self.membrane_barostat_gamma_nonzero.SurfaceTension()) + == self.modified_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) @@ -692,6 +827,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 @@ -723,6 +861,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.modified_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) + assert np.isclose(reduced_potential, state.reduced_potential(context)) + def test_method_reduced_potential_at_states(self): """ThermodynamicState.reduced_potential_at_states() method. @@ -974,6 +1129,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_xy == sampler_state.area_xy # Energies are undefined for as subset of atoms. assert sliced_sampler_state.kinetic_energy is None