From cb770745659f0a2f112299cff00c2057ae66399d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iv=C3=A1n=20Pulido?= <2949729+ijpulidos@users.noreply.github.com> Date: Fri, 26 May 2023 15:56:56 -0400 Subject: [PATCH 1/6] Default to current temperature. Add docstring. Raise meaningful exception. --- perses/samplers/multistate.py | 42 ++++++++++++++++++++++++++++++++--- 1 file changed, 39 insertions(+), 3 deletions(-) diff --git a/perses/samplers/multistate.py b/perses/samplers/multistate.py index 179cfa954..fb14ff9c9 100644 --- a/perses/samplers/multistate.py +++ b/perses/samplers/multistate.py @@ -32,8 +32,40 @@ def __init__(self, *args, hybrid_factory=None, **kwargs): # TODO: Should this overload the create() method from parent instead of being setup()? def setup(self, n_states, temperature, storage_file, minimisation_steps=100, n_replicas=None, lambda_schedule=None, - lambda_protocol=None, endstates=True, t_max=300 * unit.kelvin): - + lambda_protocol=None, endstates=True, t_max=None): + """ + Set up the simulation with the specified parameters. + + Parameters: + ----------- + n_states : int + The number of alchemical states to simulate. + temperature : openmm.unit.Quantity + The temperature of the simulation in Kelvin. + storage_file : str + The path to the storage file to store the simulation results. + minimisation_steps : int, optional + The number of minimisation steps to perform before simulation. Default is 100. + n_replicas : int, optional + The number of replicas for replica exchange. If not specified, it will be set to `n_states`. + lambda_schedule : array-like, optional + The schedule of lambda values for the alchemical states. Default is a linear schedule from 0 to 1. + lambda_protocol : object, optional + The lambda protocol object that defines the alchemical transformation protocol. Default is None. + endstates : bool, optional + Whether to generate unsampled endstates. Default is True. + t_max : openmm.unit.Quantity, optional + The maximum temperature for REST scaling. Default is None. + + Raises: + ------- + ValueError + If the hybrid factory name is not supported. + + Returns: + -------- + None + """ from perses.dispersed import feptasks # Retrieve class name, hybrid system, and hybrid positions @@ -50,11 +82,15 @@ def setup(self, n_states, temperature, storage_file, minimisation_steps=100, lambda_zero_alchemical_state = RESTCapableRelativeAlchemicalState.from_system(hybrid_system) lambda_protocol = RESTCapableLambdaProtocol() if lambda_protocol is None else lambda_protocol + # Default to current temperature if t_max is not specified (no REST scaling) + if t_max is None: + t_max = temperature + # Set beta_0 and beta_m beta_0 = 1 / (kB * temperature) beta_m = 1 / (kB * t_max) else: - raise Exception(f"{factory_name} not supported") + raise ValueError(f"{factory_name} not supported") # Create reference compound thermodynamic state thermostate = ThermodynamicState(hybrid_system, temperature=temperature) From 6e84d0c59a43e486bd94eaa1c7ba5b93dc3df26f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iv=C3=A1n=20Pulido?= <2949729+ijpulidos@users.noreply.github.com> Date: Fri, 26 May 2023 15:57:20 -0400 Subject: [PATCH 2/6] Use openmmtools CODATA 2018 constants --- perses/annihilation/relative.py | 9 +-------- 1 file changed, 1 insertion(+), 8 deletions(-) diff --git a/perses/annihilation/relative.py b/perses/annihilation/relative.py index 3f69f07aa..10cd3f3a4 100644 --- a/perses/annihilation/relative.py +++ b/perses/annihilation/relative.py @@ -2833,14 +2833,7 @@ class RESTCapableHybridTopologyFactory(HybridTopologyFactory): _new_system_exceptions : dict of key: tuple of ints, value: list of floats key: new system indices of the atoms in the exception, value: chargeProd (units of the proton charge squared), sigma (nm), and epsilon (kJ/mol) for the exception """ - - # Constants copied from: https://github.com/openmm/openmm/blob/master/platforms/reference/include/SimTKOpenMMRealType.h#L89. These will be imported directly once we have addresssed https://github.com/choderalab/openmmtools/issues/522 - M_PI = 3.14159265358979323846 - E_CHARGE = (1.602176634e-19) - AVOGADRO = (6.02214076e23) - EPSILON0 = (1e-6*8.8541878128e-12/(E_CHARGE*E_CHARGE*AVOGADRO)) - ONE_4PI_EPS0 = (1/(4*M_PI*EPSILON0)) - #from openmmtools.constants import ONE_4PI_EPS0 + from openmmtools.constants import ONE_4PI_EPS0 _default_electrostatics_expression_list = [ From 99f7384304b590f4167d1de0238c6ca2af490929 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iv=C3=A1n=20Pulido?= <2949729+ijpulidos@users.noreply.github.com> Date: Fri, 26 May 2023 16:04:42 -0400 Subject: [PATCH 3/6] Handling rest max temperature YAML parameter in CLI --- examples/new-cli/template.yaml | 6 ++++ perses/app/setup_relative_calculation.py | 40 +++++++++++++----------- 2 files changed, 27 insertions(+), 19 deletions(-) diff --git a/examples/new-cli/template.yaml b/examples/new-cli/template.yaml index e4d92cc2a..aa38c19ed 100644 --- a/examples/new-cli/template.yaml +++ b/examples/new-cli/template.yaml @@ -36,6 +36,8 @@ pressure: 1 # atmsopheres temperature: 300 # kelvin timestep: 4 # femtoseconds ionic_strength: 0.15 # molar +# Max temperature for REST scaling (optional) +max_temperature: 600 # Kelvin # Atom mapping specification atom_expression: @@ -78,3 +80,7 @@ phases: # Use geometry-derived mapping use_given_geometries: true given_geometries_tolerance: 0.4 # angstroms + +# Realtime analysis frequency +offline-freq: 100 + diff --git a/perses/app/setup_relative_calculation.py b/perses/app/setup_relative_calculation.py index 8deada142..d82da626f 100644 --- a/perses/app/setup_relative_calculation.py +++ b/perses/app/setup_relative_calculation.py @@ -3,7 +3,7 @@ import pickle import os import sys -import simtk.unit as unit +from openmm import unit import logging import warnings from cloudpathlib import AnyPath @@ -469,26 +469,28 @@ def run_setup(setup_options, serialize_systems=True, build_samplers=True): else: measure_shadow_work = False _logger.info(f"\tno measure_shadow_work specified: defaulting to False.") - if isinstance(setup_options['pressure'], (float, int)): - pressure = setup_options['pressure'] * unit.atmosphere - else: - pressure = setup_options['pressure'] - if isinstance(setup_options['temperature'], (float, int)): - temperature = setup_options['temperature'] * unit.kelvin - else: - temperature = setup_options['temperature'] - if isinstance(setup_options['solvent_padding'], (float, int)): - solvent_padding_angstroms = setup_options['solvent_padding'] * unit.angstrom - else: - solvent_padding_angstroms = setup_options['solvent_padding'] - if isinstance(setup_options['ionic_strength'], (float, int)): - ionic_strength = setup_options['ionic_strength'] * unit.molar - else: - ionic_strength = setup_options['ionic_strength'] + # Read simulation options/parameters and assign units if needed + pressure = setup_options['pressure'] + temperature = setup_options['temperature'] + solvent_padding_angstroms = setup_options['solvent_padding'] + ionic_strength = setup_options['ionic_strength'] + max_temperature = setup_options.get('max_temperature') + if isinstance(pressure, (float, int)): + pressure *= unit.atmosphere + if isinstance(temperature, (float, int)): + temperature *= unit.kelvin + if isinstance(solvent_padding_angstroms, (float, int)): + solvent_padding_angstroms *= unit.angstrom + if isinstance(ionic_strength, (float, int)): + ionic_strength *= unit.molar + if isinstance(max_temperature, (float, int)): + max_temperature *= unit.kelvin + _logger.info(f"\tsetting pressure: {pressure}.") - _logger.info(f"\tsetting temperature: {temperature}.") + _logger.info(f"\tsetting temperature: {temperature}K.") _logger.info(f"\tsetting solvent padding: {solvent_padding_angstroms}A.") _logger.info(f"\tsetting ionic strength: {ionic_strength}M.") + _logger.info(f"\tsetting max temperature: {max_temperature}K.") setup_pickle_file = setup_options['save_setup_pickle_as'] if 'save_setup_pickle_as' in list(setup_options) else None _logger.info(f"\tsetup pickle file: {setup_pickle_file}") @@ -710,7 +712,7 @@ def run_setup(setup_options, serialize_systems=True, build_samplers=True): hybrid_factory=htf[phase], online_analysis_interval=setup_options['offline-freq'], ) hss[phase].setup(n_states=n_states, temperature=temperature, storage_file=reporter, - endstates=endstates) + endstates=endstates, t_max=max_temperature) # We need to specify contexts AFTER setup hss[phase].energy_context_cache = energy_context_cache hss[phase].sampler_context_cache = sampler_context_cache From 97c97c0d1cb48a80e4dba9b46bc39e7286f07876 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iv=C3=A1n=20Pulido?= <2949729+ijpulidos@users.noreply.github.com> Date: Fri, 26 May 2023 16:11:15 -0400 Subject: [PATCH 4/6] Adding HTF input parameter to yaml template --- examples/new-cli/template.yaml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/examples/new-cli/template.yaml b/examples/new-cli/template.yaml index aa38c19ed..8a0878327 100644 --- a/examples/new-cli/template.yaml +++ b/examples/new-cli/template.yaml @@ -84,3 +84,5 @@ given_geometries_tolerance: 0.4 # angstroms # Realtime analysis frequency offline-freq: 100 +# Advanced (optional) -- Specify HTF (HybridTopologyFactory or RESTCapableHybridTopologyFactory) +# hybrid_topology_factory: RESTCapableHybridTopologyFactory \ No newline at end of file From e09f5e540ac9d2d7b21fb43ea720d6da1920d255 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iv=C3=A1n=20Pulido?= <2949729+ijpulidos@users.noreply.github.com> Date: Fri, 26 May 2023 17:26:19 -0400 Subject: [PATCH 5/6] handling units for REST parameters. REST parameters in YAML template. --- examples/new-cli/template.yaml | 9 ++++++--- perses/app/setup_relative_calculation.py | 6 ++++-- 2 files changed, 10 insertions(+), 5 deletions(-) diff --git a/examples/new-cli/template.yaml b/examples/new-cli/template.yaml index 8a0878327..87bb3fd9d 100644 --- a/examples/new-cli/template.yaml +++ b/examples/new-cli/template.yaml @@ -36,8 +36,7 @@ pressure: 1 # atmsopheres temperature: 300 # kelvin timestep: 4 # femtoseconds ionic_strength: 0.15 # molar -# Max temperature for REST scaling (optional) -max_temperature: 600 # Kelvin + # Atom mapping specification atom_expression: @@ -85,4 +84,8 @@ given_geometries_tolerance: 0.4 # angstroms offline-freq: 100 # Advanced (optional) -- Specify HTF (HybridTopologyFactory or RESTCapableHybridTopologyFactory) -# hybrid_topology_factory: RESTCapableHybridTopologyFactory \ No newline at end of file +#hybrid_topology_factory: RESTCapableHybridTopologyFactory +# REST specific parameters (optional) +#max_temperature: 600 # Kelvin +#rest_radius: 0.3 # nanometers +#w_lifting: 0.3 # nanometers diff --git a/perses/app/setup_relative_calculation.py b/perses/app/setup_relative_calculation.py index d82da626f..9e2bdc8ba 100644 --- a/perses/app/setup_relative_calculation.py +++ b/perses/app/setup_relative_calculation.py @@ -1095,11 +1095,13 @@ def _generate_htf(phase: str, topology_proposal_dictionary: dict, setup_options: # Add/use specified REST HTF parameters if present rest_specific_options = dict() try: - rest_specific_options.update({'rest_radius': setup_options['rest_radius']}) + rest_radius = setup_options['rest_radius'] * unit.nanometer + rest_specific_options.update({'rest_radius': rest_radius}) except KeyError: _logger.info("'rest_radius' not specified. Using default value.") try: - rest_specific_options.update({'w_lifting': setup_options['w_lifting']}) + w_lifting = setup_options['w_lifting'] * unit.nanometer + rest_specific_options.update({'w_lifting': w_lifting}) except KeyError: _logger.info("'w_lifting' not specified. Using default value.") From 05c204beb0ddac928601a6aa7fa5cc87539ee813 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iv=C3=A1n=20Pulido?= <2949729+ijpulidos@users.noreply.github.com> Date: Mon, 7 Aug 2023 22:39:15 -0400 Subject: [PATCH 6/6] Using not water instead for rest region --- perses/annihilation/relative.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/perses/annihilation/relative.py b/perses/annihilation/relative.py index 10cd3f3a4..ded48f7ba 100644 --- a/perses/annihilation/relative.py +++ b/perses/annihilation/relative.py @@ -3213,8 +3213,10 @@ def _generate_rest_region(self): # Retrieve neighboring atoms within self._rest_radius nm of the query atoms traj = md.Trajectory(np.array(self._hybrid_positions), self._hybrid_topology) - solute_atoms = list(traj.topology.select("is_protein")) + solute_atoms = list(traj.topology.select("not water")) rest_atoms = list(md.compute_neighbors(traj, self._rest_radius.value_in_unit_system(unit.md_unit_system), query_indices, haystack_indices=solute_atoms)[0]) + # rest_atoms = list( + # md.compute_neighbors(traj, self._rest_radius.value_in_unit_system(unit.md_unit_system), query_indices)[0]) # Retrieve full residues for all atoms in rest region residues = [atom.residue.index for atom in traj.topology.atoms if atom.index in rest_atoms]