This package can be installed using conda (or mamba, a faster version of conda):
mamba install -c conda-forge femto
+
If you are running with MPI on an HPC cluster, you may need to instruct conda to use your local installation
depending on your setup
-mamba install -c conda-forge femto "openmpi=4.1.5=*external*"
+
Getting Started#
To get started, see the usage guide.
diff --git a/latest/search/search_index.json b/latest/search/search_index.json
index dd17d9f..60c18d4 100644
--- a/latest/search/search_index.json
+++ b/latest/search/search_index.json
@@ -1 +1 @@
-{"config":{"lang":["en"],"separator":"[\\s\\-]+","pipeline":["stopWordFilter"]},"docs":[{"location":"","title":"Overview","text":"ABsolute SOLVantion Free Energy Calculations Absolute solvation free energy calculations using OpenMM
The absolv framework aims to offer a simple API for computing the change in free energy when transferring a solute from one solvent to another, or to vacuum in the case of solvation free energy calculations.
It offers two routes to this end: standard equilibrium calculations and non-equilibrium switching type calculations, where the latter will be the main focus of this framework.
Warning
This code is currently experimental and under active development. If you are using this it, please be aware that it is not guaranteed to provide correct results, the documentation and testing may be incomplete, and the API may change without notice.
"},{"location":"#installation","title":"Installation","text":"This package can be installed using conda (or mamba, a faster version of conda):
mamba install -c conda-forge femto\n
If you are running with MPI on an HPC cluster, you may need to instruct conda to use your local installation depending on your setup
mamba install -c conda-forge femto \"openmpi=4.1.5=*external*\"\n
"},{"location":"#getting-started","title":"Getting Started","text":"To get started, see the usage guide.
"},{"location":"development/","title":"Development","text":"To create a development environment, you must have mamba installed.
A development conda environment can be created and activated with:
make env\nconda activate absolv\n
To format the codebase:
make format\n
To run the unit tests:
make test\n
To serve the documentation locally:
mkdocs serve\n
"},{"location":"reproducibility/","title":"Reproducibility","text":"While every effort has been made to ensure the 'correctness and reproducibility' of any results computed using this framework, achieving consistent free energies between different frameworks and simulation engines has been notoriously tricky 1.
In an attempt to ensure that this framework remains at least self-consistent between versions, and as consistent as possible with other packages, a suite of regression tests are provided in the main GitHub repository.
These include tests to ensure that:
- the energies of systems that been alchemically modified are consistent with the original system
- systems that contain custom non-bonded force are alchemically transformed in the exact same was as systems that contain the built-in OpenMM LJ non-bonded force
and most importantly, that computing the free energies using both the 'equilibrium' and 'non-equilibrium' methods supported in this framework are in agreement amongst themselves, with values computed using Yank, and with the GROMACS values reported by Loeffler et al 1.
"},{"location":"reproducibility/#regression-results","title":"Regression Results","text":"The results of running the free energy comparisons using the latest version of the framework are shown below:
Here the error bars show the standard deviation computed across three replicas.
-
Hannes H Loeffler, Stefano Bosisio, Guilherme Duarte Ramos Matos, Donghyuk Suh, Benoit Roux, David L Mobley, and Julien Michel. Reproducibility of free energy calculations across different molecular simulation software packages. Journal of chemical theory and computation, 14(11):5567\u20135582, 2018.\u00a0\u21a9\u21a9
"},{"location":"reference/","title":"Index","text":""},{"location":"reference/#absolv","title":"absolv","text":"Absolute solvation free energy calculations using OpenMM
Modules:
-
config \u2013 Configure free energy calculations.
-
fep \u2013 Prepare OpenMM systems for FEP calculations.
-
neq \u2013 Run non-equilibrium forward and reverse sampling.
-
runner \u2013 Run calculations defined by a config.
-
setup \u2013 Setup systems ready for calculations.
-
tests \u2013 -
utils \u2013 Common utils
"},{"location":"reference/SUMMARY/","title":"SUMMARY","text":" - absolv
- config
- fep
- neq
- runner
- setup
- utils
- openmm
- topology
"},{"location":"reference/config/","title":"Configure free energy calculations.
Classes:
-
System \u2013 Define the two solvents that solutes will be transferred between (a -> b),
-
MinimizationProtocol \u2013 Configure how a system should be energy minimized.
-
SimulationProtocol \u2013 Configure how a system should be evolved by molecular simulation.
-
HREMDProtocol \u2013 Configure how a system should be evolved by Hamiltonian replica exchange.
-
EquilibriumProtocol \u2013 Configure how an equilibrium (e.g. TI, MBAR) alchemical free energy
-
SwitchingProtocol \u2013 Configure non-reversibly driving a system between to alchemical states.
-
NonEquilibriumProtocol \u2013 Configure a non-equilibrium alchemical free energy calculation [1, 2].
-
Config \u2013 Configure a transfer free energy calculation.
-
Result \u2013 The result of a free energy calculation.
"},{"location":"reference/config/#absolv.config.System","title":"System pydantic-model","text":" Bases: BaseModel
Define the two solvents that solutes will be transferred between (a -> b), as well as the solutes themselves.
Fields:
-
solutes (dict[str, PositiveInt]) -
solvent_a (dict[str, PositiveInt] | None) -
solvent_b (dict[str, PositiveInt] | None) -
n_solute_molecules (int) -
n_solvent_molecules_a (int) -
n_solvent_molecules_b (int) -
components_a (list[tuple[str, int]]) -
components_b (list[tuple[str, int]])
Validators:
-
_validate_solvent_a \u2192 solvent_a -
_validate_solvent_b \u2192 solvent_b
"},{"location":"reference/config/#absolv.config.System.solutes","title":"solutes pydantic-field","text":"solutes: dict[str, PositiveInt]\n
A dictionary containing the SMILES patterns of each solute in the system as well as how many instances of each there should be.
"},{"location":"reference/config/#absolv.config.System.solvent_a","title":"solvent_a pydantic-field","text":"solvent_a: dict[str, PositiveInt] | None\n
A dictionary containing the SMILES patterns of each component in the first solvent as well as how many instances of each there should be.A value of None should be used to indicate vacuum.
"},{"location":"reference/config/#absolv.config.System.solvent_b","title":"solvent_b pydantic-field","text":"solvent_b: dict[str, PositiveInt] | None\n
A dictionary containing the SMILES patterns of each component in the second solvent as well as how many instances of each there should be. A value of None should be used to indicate vacuum.
"},{"location":"reference/config/#absolv.config.System.n_solute_molecules","title":"n_solute_molecules pydantic-field","text":"n_solute_molecules: int\n
Returns the total number of solute molecules that will be present.
"},{"location":"reference/config/#absolv.config.System.n_solvent_molecules_a","title":"n_solvent_molecules_a pydantic-field","text":"n_solvent_molecules_a: int\n
Returns the total number of solvent molecules that will be present in the first solution.
"},{"location":"reference/config/#absolv.config.System.n_solvent_molecules_b","title":"n_solvent_molecules_b pydantic-field","text":"n_solvent_molecules_b: int\n
Returns the total number of solvent molecules that will be present in the second solution.
"},{"location":"reference/config/#absolv.config.System.components_a","title":"components_a pydantic-field","text":"components_a: list[tuple[str, int]]\n
Returns the identities and counts of the molecules present in the first system.
Returns:
-
list[tuple[str, int]] \u2013 The SMILES representation and count of each component.
"},{"location":"reference/config/#absolv.config.System.components_b","title":"components_b pydantic-field","text":"components_b: list[tuple[str, int]]\n
Returns the identities and counts of the molecules present in the second system.
Returns:
-
list[tuple[str, int]] \u2013 The SMILES representation and count of each component.
"},{"location":"reference/config/#absolv.config.MinimizationProtocol","title":"MinimizationProtocol pydantic-model","text":" Bases: BaseModel
Configure how a system should be energy minimized.
Fields:
-
tolerance (OpenMMQuantity[kilojoule_per_mole / nanometers]) -
max_iterations (NonNegativeInt)
"},{"location":"reference/config/#absolv.config.MinimizationProtocol.tolerance","title":"tolerance pydantic-field","text":"tolerance: OpenMMQuantity[\n kilojoule_per_mole / nanometers\n] = (10.0 * kilojoule_per_mole / nanometers)\n
How precisely the energy minimum must be located [kj / mol / nm]
"},{"location":"reference/config/#absolv.config.MinimizationProtocol.max_iterations","title":"max_iterations pydantic-field","text":"max_iterations: NonNegativeInt = 0\n
The maximum number of iterations to perform. If this is 0, minimization is continued until the results converge.
"},{"location":"reference/config/#absolv.config.SimulationProtocol","title":"SimulationProtocol pydantic-model","text":" Bases: BaseModel
Configure how a system should be evolved by molecular simulation.
Fields:
-
integrator (LangevinIntegrator) -
n_steps (PositiveInt)
"},{"location":"reference/config/#absolv.config.SimulationProtocol.integrator","title":"integrator pydantic-field","text":"integrator: LangevinIntegrator = LangevinIntegrator(\n timestep=2.0 * femtosecond, friction=1.0 / picosecond\n)\n
The integrator to use for the simulation.
"},{"location":"reference/config/#absolv.config.SimulationProtocol.n_steps","title":"n_steps pydantic-field","text":"n_steps: PositiveInt\n
The number of steps to evolve the system by.
"},{"location":"reference/config/#absolv.config.HREMDProtocol","title":"HREMDProtocol pydantic-model","text":" Bases: BaseModel
Configure how a system should be evolved by Hamiltonian replica exchange.
Fields:
-
integrator (LangevinIntegrator) -
n_warmup_steps (int) -
n_steps_per_cycle (int) -
n_cycles (int)
"},{"location":"reference/config/#absolv.config.HREMDProtocol.integrator","title":"integrator pydantic-field","text":"integrator: LangevinIntegrator = LangevinIntegrator(\n timestep=2.0 * femtosecond, friction=1.0 / picosecond\n)\n
The integrator to use for the simulation.
"},{"location":"reference/config/#absolv.config.HREMDProtocol.n_warmup_steps","title":"n_warmup_steps pydantic-field","text":"n_warmup_steps: int = 50000\n
The number of steps to run each replica for before starting hremd trials. All energies gathered during this period will be discarded.
"},{"location":"reference/config/#absolv.config.HREMDProtocol.n_steps_per_cycle","title":"n_steps_per_cycle pydantic-field","text":"n_steps_per_cycle: int = 1000\n
The number of steps to propagate the system by before attempting an exchange.
"},{"location":"reference/config/#absolv.config.HREMDProtocol.n_cycles","title":"n_cycles pydantic-field","text":"n_cycles: int = 2500\n
The number of cycles of 'propagate the system' -> 'exchange replicas' to run.
"},{"location":"reference/config/#absolv.config.EquilibriumProtocol","title":"EquilibriumProtocol pydantic-model","text":" Bases: BaseModel
Configure how an equilibrium (e.g. TI, MBAR) alchemical free energy calculation.
Fields:
-
type (Literal['equilibrium']) -
minimization_protocol (MinimizationProtocol) -
equilibration_protocol (SimulationProtocol) -
production_protocol (HREMDProtocol) -
lambda_sterics (list[confloat(ge=0.0, le=1.0)]) -
lambda_electrostatics (list[confloat(ge=0.0, le=1.0)]) -
n_states (int)
"},{"location":"reference/config/#absolv.config.EquilibriumProtocol.minimization_protocol","title":"minimization_protocol pydantic-field","text":"minimization_protocol: MinimizationProtocol = (\n MinimizationProtocol()\n)\n
Whether to minimize the energy of the system prior to any simulations.
"},{"location":"reference/config/#absolv.config.EquilibriumProtocol.equilibration_protocol","title":"equilibration_protocol pydantic-field","text":"equilibration_protocol: SimulationProtocol = (\n SimulationProtocol(n_steps=100000)\n)\n
The (optional) protocol that describes the equilibration simulation to run prior to the production one.
"},{"location":"reference/config/#absolv.config.EquilibriumProtocol.production_protocol","title":"production_protocol pydantic-field","text":"production_protocol: HREMDProtocol = HREMDProtocol(\n n_steps_per_cycle=6250, n_cycles=160\n)\n
The protocol that describes the production to run.
"},{"location":"reference/config/#absolv.config.EquilibriumProtocol.lambda_sterics","title":"lambda_sterics pydantic-field","text":"lambda_sterics: list[confloat(ge=0.0, le=1.0)]\n
The alchemical pathway to transform the vdW interactions along. A value of 1.0 represents a fully interacting system while a value of 0.0 represents a system with the solute-solute and solute-solvent vdW interactions disabled.
"},{"location":"reference/config/#absolv.config.EquilibriumProtocol.lambda_electrostatics","title":"lambda_electrostatics pydantic-field","text":"lambda_electrostatics: list[confloat(ge=0.0, le=1.0)]\n
The alchemical pathway to transform the electrostatic interactions along. A value of 1.0 represents a fully interacting system while a value of 0.0 represents a system with the solute-solute and solute-solvent electrostatic interactions disabled.
"},{"location":"reference/config/#absolv.config.EquilibriumProtocol.n_states","title":"n_states pydantic-field","text":"n_states: int\n
Returns the number of lambda states that will be sampled at.
"},{"location":"reference/config/#absolv.config.SwitchingProtocol","title":"SwitchingProtocol pydantic-model","text":" Bases: BaseModel
Configure non-reversibly driving a system between to alchemical states.
Fields:
-
n_electrostatic_steps (NonNegativeInt) -
n_steps_per_electrostatic_step (NonNegativeInt) -
n_steric_steps (NonNegativeInt) -
n_steps_per_steric_step (NonNegativeInt)
"},{"location":"reference/config/#absolv.config.SwitchingProtocol.n_electrostatic_steps","title":"n_electrostatic_steps pydantic-field","text":"n_electrostatic_steps: NonNegativeInt\n
The number of steps to annihilate the electrostatics interactions over. The total time needed to annihilate the electrostatics interactions will be n_electrostatic_steps * n_steps_per_electrostatic_step * timestep
"},{"location":"reference/config/#absolv.config.SwitchingProtocol.n_steps_per_electrostatic_step","title":"n_steps_per_electrostatic_step pydantic-field","text":"n_steps_per_electrostatic_step: NonNegativeInt\n
The number of timesteps to evolve the system by each time the electrostatics interactions are modified. A value of 1 will give a 'smooth' transition between the each discrete lambda value whereas a value greater than 1 will yield a stepwise transition as shown in Figure 3 of doi:10.1063/1.4712028.
"},{"location":"reference/config/#absolv.config.SwitchingProtocol.n_steric_steps","title":"n_steric_steps pydantic-field","text":"n_steric_steps: NonNegativeInt\n
The number of steps to decouple the sterics interactions over once the electrostatics interactions have been annihilated. The total time needed to annihilate the sterics interactions will be n_steric_steps * n_steps_per_steric_step * timestep
"},{"location":"reference/config/#absolv.config.SwitchingProtocol.n_steps_per_steric_step","title":"n_steps_per_steric_step pydantic-field","text":"n_steps_per_steric_step: NonNegativeInt\n
The number of timesteps to evolve the system by each time the sterics interactions are modified. A value of 1 will give a 'smooth' transition between the each discrete lambda value whereas a value greater than 1 will yield a stepwise transition as shown in Figure 3 of doi:10.1063/1.4712028.
"},{"location":"reference/config/#absolv.config.NonEquilibriumProtocol","title":"NonEquilibriumProtocol pydantic-model","text":" Bases: BaseModel
Configure a non-equilibrium alchemical free energy calculation [1, 2].
It is expected that first the electrostatics interactions will be annihilated followed by a decoupling of the sterics interactions.
References [1] Ballard, Andrew J., and Christopher Jarzynski. \"Replica exchange with nonequilibrium switches: Enhancing equilibrium sampling by increasing replica overlap.\" The Journal of chemical physics 136.19 (2012): 194101.
[2] Gapsys, Vytautas, et al. \"Large scale relative protein ligand binding affinities using non-equilibrium alchemy.\" Chemical Science 11.4 (2020): 1140-1152.
Fields:
-
type (Literal['non-equilibrium']) -
minimization_protocol (MinimizationProtocol | None) -
equilibration_protocol (SimulationProtocol | None) -
production_protocol (SimulationProtocol) -
production_report_interval (PositiveInt) -
switching_protocol (SwitchingProtocol)
"},{"location":"reference/config/#absolv.config.NonEquilibriumProtocol.minimization_protocol","title":"minimization_protocol pydantic-field","text":"minimization_protocol: MinimizationProtocol | None = (\n MinimizationProtocol()\n)\n
The (optional) protocol to follow when minimizing the system in both the end states prior to running the equilibrium simulations.
"},{"location":"reference/config/#absolv.config.NonEquilibriumProtocol.equilibration_protocol","title":"equilibration_protocol pydantic-field","text":"equilibration_protocol: SimulationProtocol | None = (\n SimulationProtocol(n_steps=100000)\n)\n
The (optional) protocol to follow when equilibrating the system in both the end states prior to running the production equilibrium simulations.
"},{"location":"reference/config/#absolv.config.NonEquilibriumProtocol.production_protocol","title":"production_protocol pydantic-field","text":"production_protocol: SimulationProtocol = (\n SimulationProtocol(n_steps=6250 * 160)\n)\n
The protocol to follow when running the production equilibrium simulation in both the end states. The snapshots generated at the end of each iteration will be used for each non-equilibrium switch.
"},{"location":"reference/config/#absolv.config.NonEquilibriumProtocol.production_report_interval","title":"production_report_interval pydantic-field","text":"production_report_interval: PositiveInt = 6250\n
The interval at which to write out the simulation state during the production simulation.
"},{"location":"reference/config/#absolv.config.NonEquilibriumProtocol.switching_protocol","title":"switching_protocol pydantic-field","text":"switching_protocol: SwitchingProtocol\n
The protocol that describes how each snapshot generated during the production simulation should be driven from state 0 -> 1 and likewise 1 -> 0 in order to compute the non-equilibrium work distributions.
"},{"location":"reference/config/#absolv.config.Config","title":"Config pydantic-model","text":" Bases: BaseModel
Configure a transfer free energy calculation.
Fields:
-
temperature (OpenMMQuantity[kelvin]) -
pressure (OpenMMQuantity[atmosphere] | None) -
alchemical_protocol_a (AlchemicalProtocol) -
alchemical_protocol_b (AlchemicalProtocol)
"},{"location":"reference/config/#absolv.config.Config.temperature","title":"temperature pydantic-field","text":"temperature: OpenMMQuantity[kelvin]\n
The temperature to calculate at [K].
"},{"location":"reference/config/#absolv.config.Config.pressure","title":"pressure pydantic-field","text":"pressure: OpenMMQuantity[atmosphere] | None\n
The pressure to calculate at [atm].
"},{"location":"reference/config/#absolv.config.Config.alchemical_protocol_a","title":"alchemical_protocol_a pydantic-field","text":"alchemical_protocol_a: AlchemicalProtocol\n
The protocol that describes the alchemical pathway to transform the solute along in the first solvent.
"},{"location":"reference/config/#absolv.config.Config.alchemical_protocol_b","title":"alchemical_protocol_b pydantic-field","text":"alchemical_protocol_b: AlchemicalProtocol\n
The protocol that describes the alchemical pathway to transform the solute along in the second solvent.
"},{"location":"reference/config/#absolv.config.Result","title":"Result pydantic-model","text":" Bases: BaseModel
The result of a free energy calculation.
Fields:
-
dg_solvent_a (OpenMMQuantity[KCAL_MOL]) -
dg_std_solvent_a (OpenMMQuantity[KCAL_MOL]) -
dg_solvent_b (OpenMMQuantity[KCAL_MOL]) -
dg_std_solvent_b (OpenMMQuantity[KCAL_MOL]) -
dg (Quantity) -
dg_std (Quantity)
"},{"location":"reference/config/#absolv.config.Result.dg_solvent_a","title":"dg_solvent_a pydantic-field","text":"dg_solvent_a: OpenMMQuantity[KCAL_MOL]\n
The change in free energy of alchemically transforming the solute from an interacting to a non-interacting state in the first solvent.
"},{"location":"reference/config/#absolv.config.Result.dg_std_solvent_a","title":"dg_std_solvent_a pydantic-field","text":"dg_std_solvent_a: OpenMMQuantity[KCAL_MOL]\n
The standard error in dg_solvent_a.
"},{"location":"reference/config/#absolv.config.Result.dg_solvent_b","title":"dg_solvent_b pydantic-field","text":"dg_solvent_b: OpenMMQuantity[KCAL_MOL]\n
The change in free energy of alchemically transforming the solute from an interacting to a non-interacting state in the second solvent.
"},{"location":"reference/config/#absolv.config.Result.dg_std_solvent_b","title":"dg_std_solvent_b pydantic-field","text":"dg_std_solvent_b: OpenMMQuantity[KCAL_MOL]\n
The standard error in dg_solvent_b.
"},{"location":"reference/config/#absolv.config.Result.dg","title":"dg pydantic-field","text":"dg: Quantity\n
The change in free energy of transferring the solute from solvent-a to solvent-b in units of kT.
"},{"location":"reference/config/#absolv.config.Result.dg_std","title":"dg_std pydantic-field","text":"dg_std: Quantity\n
The standard error in ddg.
"},{"location":"reference/fep/","title":"Prepare OpenMM systems for FEP calculations.
Functions:
-
apply_fep \u2013 Generate a system whereby a number of the molecules can be alchemically
-
set_fep_lambdas \u2013 Set the values of the alchemical lambdas on an OpenMM context.
"},{"location":"reference/fep/#absolv.fep.apply_fep","title":"apply_fep","text":"apply_fep(\n system: System,\n alchemical_indices: list[set[int]],\n persistent_indices: list[set[int]],\n custom_alchemical_potential: str | None = None,\n) -> System\n
Generate a system whereby a number of the molecules can be alchemically transformed from a base chemical system.
Notes -
Currently only OpenMM systems that have:
-
vdW + electrostatics in a single built-in non-bonded force
- electrostatics in a built-in non-bonded force, vdW in a custom non-bonded force and vdW 1-4 interactions in a custom bond force
- all of the above sans any electrostatics
are supported. * By default a soft-core version of the LJ potential with a-b-c of 1-1-6 and alpha=0.5 that can be scaled by a global lambda_sterics parameter will be used for alchemical-chemical vdW interactions embedded in an OpenMM NonbondedForce while the energy expression set on a CustomNonbondedForce will be be modified to have the form \"lambda_sterics*({original_expression})\".
Parameters:
-
system (System) \u2013 The chemical system to generate the alchemical system from
-
alchemical_indices (list[set[int]]) \u2013 The atom indices corresponding to each molecule that should be alchemically transformable. The atom indices must correspond to all atoms in each molecule as alchemically transforming part of a molecule is not supported.
-
persistent_indices (list[set[int]]) \u2013 The atom indices corresponding to each molecule that should not be alchemically transformable.
-
custom_alchemical_potential (str | None, default: None ) \u2013 A custom expression to use for the potential energy function that describes the chemical-alchemical intermolecular interactions. See the Notes for information about the default value.
The expression must include \"lambda_sterics\".
Source code in absolv/fep.py def apply_fep(\n system: openmm.System,\n alchemical_indices: list[set[int]],\n persistent_indices: list[set[int]],\n custom_alchemical_potential: str | None = None,\n) -> openmm.System:\n \"\"\"Generate a system whereby a number of the molecules can be alchemically\n transformed from a base chemical system.\n\n Notes:\n * Currently only OpenMM systems that have:\n\n - vdW + electrostatics in a single built-in non-bonded force\n - electrostatics in a built-in non-bonded force, vdW in a custom\n **non-bonded** force and vdW 1-4 interactions in a custom **bond** force\n - all of the above sans any electrostatics\n\n are supported.\n * By default a soft-core version of the LJ potential with a-b-c of 1-1-6\n and alpha=0.5 that can be scaled by a global `lambda_sterics` parameter\n will be used for alchemical-chemical vdW interactions embedded in an\n OpenMM ``NonbondedForce`` while the energy expression set on a\n ``CustomNonbondedForce`` will be be modified to have the form\n ``\"lambda_sterics*({original_expression})\"``.\n\n Args:\n system: The chemical system to generate the alchemical system from\n alchemical_indices: The atom indices corresponding to each molecule that\n should be alchemically transformable. The atom indices **must**\n correspond to **all** atoms in each molecule as alchemically\n transforming part of a molecule is not supported.\n persistent_indices: The atom indices corresponding to each molecule that\n should **not** be alchemically transformable.\n custom_alchemical_potential: A custom expression to use for the potential\n energy function that describes the chemical-alchemical intermolecular\n interactions. See the Notes for information about the default value.\n\n The expression **must** include ``\"lambda_sterics\"``.\n \"\"\"\n\n system = copy.deepcopy(system)\n\n # Make sure we track v-sites attached to any solutes that may be alchemically\n # turned off. We do this as a post-process step as the OpenFF toolkit does not\n # currently expose a clean way to access this information.\n atom_indices = alchemical_indices + persistent_indices\n atom_indices = _find_v_sites(system, atom_indices)\n\n alchemical_indices = atom_indices[: len(alchemical_indices)]\n persistent_indices = atom_indices[len(alchemical_indices) :]\n\n (\n nonbonded_force,\n custom_nonbonded_force,\n custom_bond_force,\n ) = _find_nonbonded_forces(system)\n\n if nonbonded_force is not None:\n _add_electrostatics_lambda(nonbonded_force, alchemical_indices)\n\n if custom_nonbonded_force is not None:\n for alchemical_force in _add_custom_vdw_lambda(\n custom_nonbonded_force,\n alchemical_indices,\n persistent_indices,\n custom_alchemical_potential,\n ):\n system.addForce(alchemical_force)\n\n elif nonbonded_force is not None:\n for alchemical_force in _add_lj_vdw_lambda(\n nonbonded_force,\n alchemical_indices,\n persistent_indices,\n custom_alchemical_potential,\n ):\n system.addForce(alchemical_force)\n\n else:\n raise NotImplementedError\n\n return system\n
"},{"location":"reference/fep/#absolv.fep.set_fep_lambdas","title":"set_fep_lambdas","text":"set_fep_lambdas(\n context: Context,\n lambda_sterics: float | None = None,\n lambda_electrostatics: float | None = None,\n)\n
Set the values of the alchemical lambdas on an OpenMM context.
Parameters:
-
context (Context) \u2013 The context to update.
-
lambda_sterics (float | None, default: None ) \u2013 The (optional) value of the steric lambda.
-
lambda_electrostatics (float | None, default: None ) \u2013 The (optional) value of the electrostatics lambda.
Source code in absolv/fep.py def set_fep_lambdas(\n context: openmm.Context,\n lambda_sterics: float | None = None,\n lambda_electrostatics: float | None = None,\n):\n \"\"\"Set the values of the alchemical lambdas on an OpenMM context.\n\n Args:\n context: The context to update.\n lambda_sterics: The (optional) value of the steric lambda.\n lambda_electrostatics: The (optional) value of the electrostatics lambda.\n \"\"\"\n\n if lambda_sterics is not None:\n assert (\n 0.0 <= lambda_sterics <= 1.0\n ), f\"`{LAMBDA_STERICS}` must be between 0 and 1\"\n context.setParameter(LAMBDA_STERICS, lambda_sterics)\n\n if lambda_electrostatics is not None:\n assert (\n 0.0 <= lambda_electrostatics <= 1.0\n ), f\"`{LAMBDA_ELECTROSTATICS}` must be between 0 and 1\"\n\n context.setParameter(LAMBDA_ELECTROSTATICS, lambda_electrostatics)\n
"},{"location":"reference/neq/","title":"Run non-equilibrium forward and reverse sampling.
Functions:
-
run_neq \u2013 Run a non-equilibrium simulation with OpenMM, whereby a system is non-reversibly
"},{"location":"reference/neq/#absolv.neq.run_neq","title":"run_neq","text":"run_neq(\n simulation: Simulation,\n coords_0: State,\n coords_1: State,\n protocol: SwitchingProtocol,\n) -> tuple[float, float]\n
Run a non-equilibrium simulation with OpenMM, whereby a system is non-reversibly pulled along an alchemical pathway as described by Ballard and Jarzynski [1] (Figure 3) and Gapsys et al [2].
Both the forward and reverse directions will be simulated.
References [1] Ballard, Andrew J., and Christopher Jarzynski. \"Replica exchange with nonequilibrium switches: Enhancing equilibrium sampling by increasing replica overlap.\" The Journal of chemical physics 136.19 (2012): 194101.
[2] Gapsys, Vytautas, et al. \"Large scale relative protein ligand binding affinities using non-equilibrium alchemy.\" Chemical Science 11.4 (2020): 1140-1152.
Returns:
-
tuple[float, float] \u2013 The forward and reverse work values.
Source code in absolv/neq.py def run_neq(\n simulation: openmm.app.Simulation,\n coords_0: openmm.State,\n coords_1: openmm.State,\n protocol: absolv.config.SwitchingProtocol,\n) -> tuple[float, float]:\n \"\"\"Run a non-equilibrium simulation with OpenMM, whereby a system is non-reversibly\n pulled along an alchemical pathway as described by Ballard and Jarzynski [1]\n (Figure 3) and Gapsys et al [2].\n\n Both the forward and reverse directions will be simulated.\n\n References:\n [1] Ballard, Andrew J., and Christopher Jarzynski. \"Replica exchange with\n nonequilibrium switches: Enhancing equilibrium sampling by increasing replica\n overlap.\" The Journal of chemical physics 136.19 (2012): 194101.\n\n [2] Gapsys, Vytautas, et al. \"Large scale relative protein ligand binding\n affinities using non-equilibrium alchemy.\" Chemical Science 11.4 (2020):\n 1140-1152.\n\n Returns:\n The forward and reverse work values.\n \"\"\"\n\n forward_potentials = _simulate(\n simulation, coords_0, protocol, reverse_direction=False\n )\n reverse_potentials = _simulate(\n simulation, coords_1, protocol, reverse_direction=True\n )\n\n forward_work = (forward_potentials[:, 1] - forward_potentials[:, 0]).sum()\n reverse_work = (reverse_potentials[:, 1] - reverse_potentials[:, 0]).sum()\n\n return forward_work, reverse_work\n
"},{"location":"reference/runner/","title":"Run calculations defined by a config.
Classes:
-
PreparedSystem \u2013 A container for the prepared inputs for a particular solvent phase.
Functions:
-
setup \u2013 Prepare each system to be simulated, namely the ligand in each solvent.
-
run_eq \u2013 Perform a simulation at each lambda window and for each solvent.
-
run_neq \u2013 Performs the simulations required to estimate the free energy using a
"},{"location":"reference/runner/#absolv.runner.PreparedSystem","title":"PreparedSystem","text":" Bases: NamedTuple
A container for the prepared inputs for a particular solvent phase.
Attributes:
-
system (System) \u2013 The alchemically modified OpenMM system.
-
topology (Topology) \u2013 The OpenFF topology with any box vectors set.
-
coords (Quantity) \u2013 The coordinates of the system.
"},{"location":"reference/runner/#absolv.runner.PreparedSystem.system","title":"system instance-attribute","text":"system: System\n
The alchemically modified OpenMM system.
"},{"location":"reference/runner/#absolv.runner.PreparedSystem.topology","title":"topology instance-attribute","text":"topology: Topology\n
The OpenFF topology with any box vectors set.
"},{"location":"reference/runner/#absolv.runner.PreparedSystem.coords","title":"coords instance-attribute","text":"coords: Quantity\n
The coordinates of the system.
"},{"location":"reference/runner/#absolv.runner.setup","title":"setup","text":"setup(\n system: System,\n config: Config,\n force_field: ForceField | SystemGenerator,\n custom_alchemical_potential: str | None = None,\n) -> tuple[PreparedSystem, PreparedSystem]\n
Prepare each system to be simulated, namely the ligand in each solvent.
Parameters:
-
system (System) \u2013 The system to prepare.
-
config (Config) \u2013 The config defining the calculation to perform.
-
force_field (ForceField | SystemGenerator) \u2013 The force field, or a callable that transforms an OpenFF topology into an OpenMM system, without any alchemical modifications to run the calculations using.
If a callable is specified, it should take arguments of an OpenFF topology, a unit wrapped numpy array of atom coords, and a string literal with a value of either \"solvent-a\" or \"solvent-b\".
-
custom_alchemical_potential (str | None, default: None ) \u2013 A custom expression to use for the potential energy function that describes the chemical-alchemical intermolecular interactions.
See the absolv.fep.apply_fep function for more details.
Returns:
-
tuple[PreparedSystem, PreparedSystem] \u2013 The two prepared systems, corresponding to solvent-a and solvent-b respectively.
Source code in absolv/runner.py def setup(\n system: absolv.config.System,\n config: absolv.config.Config,\n force_field: openff.toolkit.ForceField | absolv.utils.openmm.SystemGenerator,\n custom_alchemical_potential: str | None = None,\n) -> tuple[PreparedSystem, PreparedSystem]:\n \"\"\"Prepare each system to be simulated, namely the ligand in each solvent.\n\n Args:\n system: The system to prepare.\n config: The config defining the calculation to perform.\n force_field: The force field, or a callable that transforms an OpenFF\n topology into an OpenMM system, **without** any alchemical modifications\n to run the calculations using.\n\n If a callable is specified, it should take arguments of an OpenFF\n topology, a unit wrapped numpy array of atom coords, and a string\n literal with a value of either ``\"solvent-a\"`` or ``\"solvent-b\"``.\n custom_alchemical_potential: A custom expression to use for the potential\n energy function that describes the chemical-alchemical intermolecular\n interactions.\n\n See the ``absolv.fep.apply_fep`` function for more details.\n\n Returns:\n The two prepared systems, corresponding to solvent-a and solvent-b respectively.\n \"\"\"\n\n solvated_a = _setup_solvent(\n \"solvent-a\",\n system.components_a,\n force_field,\n system.n_solute_molecules,\n system.n_solvent_molecules_a,\n custom_alchemical_potential,\n )\n solvated_b = _setup_solvent(\n \"solvent-b\",\n system.components_b,\n force_field,\n system.n_solute_molecules,\n system.n_solvent_molecules_b,\n custom_alchemical_potential,\n )\n\n if system.solvent_a is not None and config.pressure is not None:\n absolv.utils.openmm.add_barostat(\n solvated_a.system, config.temperature, config.pressure\n )\n if system.solvent_b is not None and config.pressure is not None:\n absolv.utils.openmm.add_barostat(\n solvated_b.system, config.temperature, config.pressure\n )\n\n return solvated_a, solvated_b\n
"},{"location":"reference/runner/#absolv.runner.run_eq","title":"run_eq","text":"run_eq(\n config: Config,\n prepared_system_a: PreparedSystem,\n prepared_system_b: PreparedSystem,\n platform: OpenMMPlatform = femto.md.constants.OpenMMPlatform.CUDA,\n output_dir: Path | None = None,\n) -> Result\n
Perform a simulation at each lambda window and for each solvent.
Parameters:
-
config (Config) \u2013 The config defining the calculation to perform.
-
prepared_system_a (PreparedSystem) \u2013 The prepared system a. See setup for more details.
-
prepared_system_b (PreparedSystem) \u2013 The prepared system b. See setup for more details.
-
platform (OpenMMPlatform, default: CUDA ) \u2013 The OpenMM platform to run using.
-
output_dir (Path | None, default: None ) \u2013 The (optional) directory to save HREMD samples to.
Source code in absolv/runner.py def run_eq(\n config: absolv.config.Config,\n prepared_system_a: PreparedSystem,\n prepared_system_b: PreparedSystem,\n platform: femto.md.constants.OpenMMPlatform = femto.md.constants.OpenMMPlatform.CUDA,\n output_dir: pathlib.Path | None = None,\n) -> absolv.config.Result:\n \"\"\"Perform a simulation at each lambda window and for each solvent.\n\n Args:\n config: The config defining the calculation to perform.\n prepared_system_a: The prepared system a. See ``setup`` for more details.\n prepared_system_b: The prepared system b. See ``setup`` for more details.\n platform: The OpenMM platform to run using.\n output_dir: The (optional) directory to save HREMD samples to.\n \"\"\"\n\n results_a, overlap_a = _run_phase_hremd(\n config.alchemical_protocol_a,\n config.temperature,\n prepared_system_a,\n platform,\n None if output_dir is None else output_dir / \"solvant-a\",\n )\n\n dg_a, dg_a_std = results_a[\"ddG_kcal_mol\"], results_a[\"ddG_error_kcal_mol\"]\n # overlap_a = overlap_a[\"overlap_0\"]\n\n results_b, overlap_b = _run_phase_hremd(\n config.alchemical_protocol_b,\n config.temperature,\n prepared_system_b,\n platform,\n None if output_dir is None else output_dir / \"solvant-b\",\n )\n\n dg_b, dg_b_std = results_b[\"ddG_kcal_mol\"], results_b[\"ddG_error_kcal_mol\"]\n # overlap_b = overlap_b[\"overlap_0\"]\n\n return absolv.config.Result(\n dg_solvent_a=dg_a * openmm.unit.kilocalorie_per_mole,\n dg_std_solvent_a=dg_a_std * openmm.unit.kilocalorie_per_mole,\n dg_solvent_b=dg_b * openmm.unit.kilocalorie_per_mole,\n dg_std_solvent_b=dg_b_std * openmm.unit.kilocalorie_per_mole,\n )\n
"},{"location":"reference/runner/#absolv.runner.run_neq","title":"run_neq","text":"run_neq(\n config: Config,\n prepared_system_a: PreparedSystem,\n prepared_system_b: PreparedSystem,\n platform: OpenMMPlatform = femto.md.constants.OpenMMPlatform.CUDA,\n) -> Result\n
Performs the simulations required to estimate the free energy using a non-equilibrium method.
These include equilibrium simulations at the two end states (i.e. fully interacting and fully de-coupled solute) for each solvent followed by non-equilibrium switching simulations between each end state to compute the forward and reverse work values.
Parameters:
-
config (Config) \u2013 The config defining the calculation to perform.
-
prepared_system_a (PreparedSystem) \u2013 The prepared system a. See setup for more details.
-
prepared_system_b (PreparedSystem) \u2013 The prepared system b. See setup for more details.
-
platform (OpenMMPlatform, default: CUDA ) \u2013 The OpenMM platform to run using.
Source code in absolv/runner.py def run_neq(\n config: absolv.config.Config,\n prepared_system_a: PreparedSystem,\n prepared_system_b: PreparedSystem,\n platform: femto.md.constants.OpenMMPlatform = femto.md.constants.OpenMMPlatform.CUDA,\n) -> absolv.config.Result:\n \"\"\"Performs the simulations required to estimate the free energy using a\n non-equilibrium method.\n\n These include **equilibrium** simulations at the two end states (i.e. fully\n interacting and fully de-coupled solute) for each solvent followed by\n non-equilibrium switching simulations between each end state to compute the\n forward and reverse work values.\n\n Args:\n config: The config defining the calculation to perform.\n prepared_system_a: The prepared system a. See ``setup`` for more details.\n prepared_system_b: The prepared system b. See ``setup`` for more details.\n platform: The OpenMM platform to run using.\n \"\"\"\n\n with tempfile.TemporaryDirectory() as tmp_dir:\n tmp_dir = pathlib.Path(tmp_dir)\n\n solvent_a_dir = tmp_dir / \"solvent-a\"\n solvent_b_dir = tmp_dir / \"solvent-b\"\n\n _run_phase_end_states(\n config.alchemical_protocol_a,\n config.temperature,\n prepared_system_a,\n solvent_a_dir,\n platform,\n )\n dg_a, dg_std_a = _run_switching(\n config.alchemical_protocol_a,\n config.temperature,\n prepared_system_a,\n solvent_a_dir,\n platform,\n )\n\n _run_phase_end_states(\n config.alchemical_protocol_b,\n config.temperature,\n prepared_system_b,\n solvent_b_dir,\n platform,\n )\n dg_b, dg_std_b = _run_switching(\n config.alchemical_protocol_b,\n config.temperature,\n prepared_system_b,\n solvent_b_dir,\n platform,\n )\n\n return absolv.config.Result(\n dg_solvent_a=dg_a.in_units_of(openmm.unit.kilocalories_per_mole),\n dg_std_solvent_a=dg_std_a.in_units_of(openmm.unit.kilocalories_per_mole),\n dg_solvent_b=dg_b.in_units_of(openmm.unit.kilocalories_per_mole),\n dg_std_solvent_b=dg_std_b.in_units_of(openmm.unit.kilocalories_per_mole),\n )\n
"},{"location":"reference/setup/","title":"Setup systems ready for calculations.
Functions:
-
setup_system \u2013 Generate a set of molecule coordinate by using the PACKMOL package.
"},{"location":"reference/setup/#absolv.setup.setup_system","title":"setup_system","text":"setup_system(\n components: list[tuple[str, int]],\n box_target_density: Quantity = 0.95 * _G_PER_ML,\n box_scale_factor: float = 1.1,\n box_padding: Quantity = 2.0 * openmm.unit.angstrom,\n tolerance: Quantity = 2.0 * openmm.unit.angstrom,\n) -> tuple[Topology, Quantity]\n
Generate a set of molecule coordinate by using the PACKMOL package.
Parameters:
-
components (list[tuple[str, int]]) \u2013 A list of the form components[i] = (smiles_i, count_i) where smiles_i is the SMILES representation of component i and count_i is the number of corresponding instances of that component to create.
-
box_target_density (Quantity, default: 0.95 * _G_PER_ML ) \u2013 Target mass density when approximating the box size for the final system with units compatible with g / mL.
-
box_scale_factor (float, default: 1.1 ) \u2013 The amount to scale the approximate box size by.
-
box_padding (Quantity, default: 2.0 * angstrom ) \u2013 The amount of extra padding to add to the box size to avoid PBC issues in units compatible with angstroms.
-
tolerance (Quantity, default: 2.0 * angstrom ) \u2013 The minimum spacing between molecules during packing in units compatible with angstroms.
Returns:
-
tuple[Topology, Quantity] \u2013 A topology containing the molecules the coordinates were generated for and a unit [A] wrapped numpy array of coordinates with shape=(n_atoms, 3).
Source code in absolv/setup.py def setup_system(\n components: list[tuple[str, int]],\n box_target_density: openmm.unit.Quantity = 0.95 * _G_PER_ML,\n box_scale_factor: float = 1.1,\n box_padding: openmm.unit.Quantity = 2.0 * openmm.unit.angstrom,\n tolerance: openmm.unit.Quantity = 2.0 * openmm.unit.angstrom,\n) -> tuple[openff.toolkit.Topology, openmm.unit.Quantity]:\n \"\"\"Generate a set of molecule coordinate by using the PACKMOL package.\n\n Args:\n components: A list of the form ``components[i] = (smiles_i, count_i)`` where\n ``smiles_i`` is the SMILES representation of component `i` and\n ``count_i`` is the number of corresponding instances of that component\n to create.\n box_target_density: Target mass density when approximating the box size for the\n final system with units compatible with g / mL.\n box_scale_factor: The amount to scale the approximate box size by.\n box_padding: The amount of extra padding to add to the box size to avoid PBC\n issues in units compatible with angstroms.\n tolerance: The minimum spacing between molecules during packing in units\n compatible with angstroms.\n\n Raises:\n * PACKMOLRuntimeError\n\n Returns:\n A topology containing the molecules the coordinates were generated for and\n a unit [A] wrapped numpy array of coordinates with shape=(n_atoms, 3).\n \"\"\"\n\n packmol_path = shutil.which(\"packmol\")\n\n if packmol_path is None:\n raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), \"packmol\")\n\n box_size = (\n _approximate_box_size_by_density(components, box_target_density)\n * box_scale_factor\n )\n\n molecules = {}\n\n for smiles, _ in components:\n if smiles in molecules:\n continue\n\n molecule = openff.toolkit.Molecule.from_smiles(smiles)\n molecule.generate_conformers(n_conformers=1)\n molecule.name = f\"component-{len(molecules)}.xyz\"\n molecules[smiles] = molecule\n\n with openff.utilities.temporary_cd():\n for molecule in molecules.values():\n molecule.to_file(molecule.name, \"xyz\")\n\n input_file_contents = _generate_input_file(\n [(molecules[smiles].name, count) for smiles, count in components],\n box_size,\n tolerance,\n )\n\n with open(\"input.txt\", \"w\") as file:\n file.write(input_file_contents)\n\n with open(\"input.txt\") as file:\n subprocess.run(packmol_path, stdin=file, check=True, capture_output=True)\n\n with open(\"output.xyz\") as file:\n output_lines = file.read().splitlines(False)\n\n coordinates = (\n numpy.array(\n [\n [float(coordinate) for coordinate in coordinate_line.split()[1:]]\n for coordinate_line in output_lines[2:]\n if len(coordinate_line) > 0\n ]\n )\n * openmm.unit.angstrom\n )\n\n topology = openff.toolkit.Topology.from_molecules(\n [molecules[smiles] for smiles, count in components for _ in range(count)]\n )\n topology.box_vectors = numpy.eye(3) * (box_size + box_padding)\n\n return topology, coordinates\n
"},{"location":"reference/utils/","title":"Index","text":""},{"location":"reference/utils/#absolv.utils","title":"utils","text":"Common utils
Modules:
-
openmm \u2013 Utilities to manipulate OpenMM objects.
-
topology \u2013 Utilities for manipulating OpenFF topology objects.
"},{"location":"reference/utils/openmm/","title":"Utilities to manipulate OpenMM objects.
Functions:
-
add_barostat \u2013 Add a barostat to a system in-place.
-
create_simulation \u2013 Creates an OpenMM simulation object.
-
create_system_generator \u2013 Creates a 'system generator' that can be used when setting up an alchemical
-
extract_frame \u2013 Extracts a frame from a trajectory as an OpenMM state object.
"},{"location":"reference/utils/openmm/#absolv.utils.openmm.add_barostat","title":"add_barostat","text":"add_barostat(\n system: System,\n temperature: Quantity,\n pressure: Quantity,\n frequency: int = 25,\n)\n
Add a barostat to a system in-place.
Parameters:
-
system (System) \u2013 The system to add the barostat to.
-
temperature (Quantity) \u2013 The temperature to simulate at.
-
pressure (Quantity) \u2013 The pressure to simulate at.
-
frequency (int, default: 25 ) \u2013 The frequency at which to apply the barostat.
Source code in absolv/utils/openmm.py def add_barostat(\n system: openmm.System,\n temperature: openmm.unit.Quantity,\n pressure: openmm.unit.Quantity,\n frequency: int = 25,\n):\n \"\"\"Add a barostat to a system in-place.\n\n Args:\n system: The system to add the barostat to.\n temperature: The temperature to simulate at.\n pressure: The pressure to simulate at.\n frequency: The frequency at which to apply the barostat.\n \"\"\"\n\n barostats = [\n force\n for force in system.getForces()\n if isinstance(force, openmm.MonteCarloBarostat)\n ]\n assert len(barostats) == 0, \"the system should not already contain a barostat\"\n\n system.addForce(openmm.MonteCarloBarostat(pressure, temperature, frequency))\n
"},{"location":"reference/utils/openmm/#absolv.utils.openmm.create_simulation","title":"create_simulation","text":"create_simulation(\n system: System,\n topology: Topology,\n coords: Quantity,\n integrator: Integrator,\n platform: OpenMMPlatform,\n) -> Simulation\n
Creates an OpenMM simulation object.
Parameters:
-
system (System) \u2013 The system to simulate
-
topology (Topology) \u2013 The topology being simulated.
-
coords (Quantity) \u2013 The initial coordinates. Box vectors (if any) will be taken from the topology.
-
integrator (Integrator) \u2013 The integrator to evolve the system with.
-
platform (OpenMMPlatform) \u2013 The accelerator to run using.
Returns:
-
Simulation \u2013 The created simulation.
Source code in absolv/utils/openmm.py def create_simulation(\n system: openmm.System,\n topology: openff.toolkit.Topology,\n coords: openmm.unit.Quantity,\n integrator: openmm.Integrator,\n platform: femto.md.constants.OpenMMPlatform,\n) -> openmm.app.Simulation:\n \"\"\"Creates an OpenMM simulation object.\n\n Args:\n system: The system to simulate\n topology: The topology being simulated.\n coords: The initial coordinates. Box vectors (if any) will be taken from the\n topology.\n integrator: The integrator to evolve the system with.\n platform: The accelerator to run using.\n\n Returns:\n The created simulation.\n \"\"\"\n platform_properties = (\n {\"Precision\": \"mixed\"} if platform.upper() in [\"CUDA\", \"OPENCL\"] else {}\n )\n platform = openmm.Platform.getPlatformByName(platform)\n\n if topology.box_vectors is not None:\n system.setDefaultPeriodicBoxVectors(*topology.box_vectors.to_openmm())\n\n simulation = openmm.app.Simulation(\n topology.to_openmm(), system, integrator, platform, platform_properties\n )\n\n if topology.box_vectors is not None:\n simulation.context.setPeriodicBoxVectors(*topology.box_vectors.to_openmm())\n\n simulation.context.setPositions(coords)\n simulation.context.setVelocitiesToTemperature(integrator.getTemperature())\n\n return simulation\n
"},{"location":"reference/utils/openmm/#absolv.utils.openmm.create_system_generator","title":"create_system_generator","text":"create_system_generator(\n force_field: ForceField,\n solvent_a_nonbonded_method: int,\n solvent_b_nonbonded_method: int,\n nonbonded_cutoff: Quantity = 1.0\n * openmm.unit.nanometer,\n constraints: int | None = None,\n rigid_water: bool | None = None,\n remove_cmm_motion: bool = True,\n hydrogen_mass: Quantity | None = None,\n switch_distance: Quantity | None = None,\n) -> SystemGenerator\n
Creates a 'system generator' that can be used when setting up an alchemical free energy calculation from an OpenMM force field.
Parameters:
-
force_field (ForceField) \u2013 The OpenMM force field to parameterize the topology using.
-
solvent_a_nonbonded_method (int) \u2013 The non-bonded method to use in solvent a.
-
solvent_b_nonbonded_method (int) \u2013 The non-bonded method to use in solvent b.
-
nonbonded_cutoff (Quantity, default: 1.0 * nanometer ) \u2013 The non-bonded cutoff to use.
-
constraints (int | None, default: None ) \u2013 The type of constraints to apply to the system.
-
rigid_water (bool | None, default: None ) \u2013 Whether to force rigid water.
-
remove_cmm_motion (bool, default: True ) \u2013 Whether to remove any CMM motion.
-
hydrogen_mass (Quantity | None, default: None ) \u2013 The mass to use for hydrogens.
-
switch_distance (Quantity | None, default: None ) \u2013 The switch distance to use.
Returns:
-
SystemGenerator \u2013 A callable that will create an OpenMM system from an OpenFF topology and the name of the solvent (i.e. \"solvent-a\" or \"solvent-b\") the system will be used for.
Source code in absolv/utils/openmm.py def create_system_generator(\n force_field: openmm.app.ForceField,\n solvent_a_nonbonded_method: int,\n solvent_b_nonbonded_method: int,\n nonbonded_cutoff: openmm.unit.Quantity = 1.0 * openmm.unit.nanometer,\n constraints: int | None = None,\n rigid_water: bool | None = None,\n remove_cmm_motion: bool = True,\n hydrogen_mass: openmm.unit.Quantity | None = None,\n switch_distance: openmm.unit.Quantity | None = None,\n) -> SystemGenerator:\n \"\"\"Creates a 'system generator' that can be used when setting up an alchemical\n free energy calculation from an OpenMM force field.\n\n Args:\n force_field: The OpenMM force field to parameterize the topology using.\n solvent_a_nonbonded_method: The non-bonded method to use in solvent a.\n solvent_b_nonbonded_method: The non-bonded method to use in solvent b.\n nonbonded_cutoff: The non-bonded cutoff to use.\n constraints: The type of constraints to apply to the system.\n rigid_water: Whether to force rigid water.\n remove_cmm_motion: Whether to remove any CMM motion.\n hydrogen_mass: The mass to use for hydrogens.\n switch_distance: The switch distance to use.\n\n Returns:\n A callable that will create an OpenMM system from an OpenFF topology and the\n name of the solvent (i.e. ``\"solvent-a\"`` or ``\"solvent-b\"``) the system will\n be used for.\n \"\"\"\n\n def system_generator(\n topology: openff.toolkit.Topology,\n coordinates: openmm.unit.Quantity,\n solvent_idx: typing.Literal[\"solvent-a\", \"solvent-b\"],\n ) -> openmm.System:\n openmm_topology = topology.to_openmm()\n\n if topology.box_vectors is not None:\n openmm_topology.setPeriodicBoxVectors(topology.box_vectors.to_openmm())\n\n # We need to fix the special case of water in order for OMM to correctly apply\n # a constraint between H atoms.\n for chain in openmm_topology.chains():\n for residue in chain.residues():\n if len(residue) != 3:\n continue\n\n symbols = sorted(atom.element.symbol for atom in residue.atoms())\n\n if symbols == [\"H\", \"H\", \"O\"]:\n residue.name = \"HOH\"\n\n from openmm.app import Modeller\n\n modeller = Modeller(\n openmm_topology,\n [\n openmm.Vec3(coordinate[0], coordinate[1], coordinate[2])\n for coordinate in coordinates.value_in_unit(openmm.unit.nanometers)\n ]\n * openmm.unit.nanometers,\n )\n modeller.addExtraParticles(force_field)\n\n system = force_field.createSystem(\n modeller.getTopology(),\n nonbondedMethod=(\n solvent_a_nonbonded_method\n if solvent_idx == \"solvent-a\"\n else solvent_b_nonbonded_method\n ),\n nonbondedCutoff=nonbonded_cutoff,\n constraints=constraints,\n rigidWater=rigid_water,\n removeCMMotion=remove_cmm_motion,\n hydrogenMass=hydrogen_mass,\n switchDistance=switch_distance,\n )\n\n return system\n\n return system_generator\n
"},{"location":"reference/utils/openmm/#absolv.utils.openmm.extract_frame","title":"extract_frame","text":"extract_frame(trajectory: Trajectory, idx: int) -> State\n
Extracts a frame from a trajectory as an OpenMM state object.
Parameters:
-
trajectory (Trajectory) \u2013 The trajectory to extract the frame from.
-
idx (int) \u2013 The index of the frame to extract.
Returns:
-
State \u2013 The extracted frame.
Source code in absolv/utils/openmm.py def extract_frame(trajectory: mdtraj.Trajectory, idx: int) -> openmm.State:\n \"\"\"Extracts a frame from a trajectory as an OpenMM state object.\n\n Args:\n trajectory: The trajectory to extract the frame from.\n idx: The index of the frame to extract.\n\n Returns:\n The extracted frame.\n \"\"\"\n\n system = openmm.System()\n\n for _ in range(trajectory.n_atoms):\n system.addParticle(1.0 * openmm.unit.dalton)\n\n context = openmm.Context(system, openmm.VerletIntegrator(0.0001))\n context.setPositions(trajectory.openmm_positions(idx))\n\n if trajectory.unitcell_vectors is not None:\n context.setPeriodicBoxVectors(*trajectory.openmm_boxes(idx))\n\n return context.getState(getPositions=True)\n
"},{"location":"reference/utils/topology/","title":"Utilities for manipulating OpenFF topology objects.
Functions:
-
topology_to_components \u2013 A helper method for condensing a topology down to a list of components
-
topology_to_atom_indices \u2013 A helper method for extracting the sets of atom indices associated with each
"},{"location":"reference/utils/topology/#absolv.utils.topology.topology_to_components","title":"topology_to_components","text":"topology_to_components(\n topology: Topology,\n) -> list[tuple[str, int]]\n
A helper method for condensing a topology down to a list of components and their counts.
Notes If the topology is not contiguous then the returned list may contain multiple tuples with the same smiles but different counts.
Parameters:
-
topology (Topology) \u2013 The topology to condense.
Returns:
-
list[tuple[str, int]] \u2013 A list of the form components[i] = (smiles_i, count_i) where smiles_i is the SMILES representation of component i and count_i is the number of corresponding instances of that component to create.
Source code in absolv/utils/topology.py def topology_to_components(topology: openff.toolkit.Topology) -> list[tuple[str, int]]:\n \"\"\"A helper method for condensing a topology down to a list of components\n and their counts.\n\n Notes:\n If the topology is not contiguous then the returned list may contain multiple\n tuples with the same smiles but different counts.\n\n Args:\n topology: The topology to condense.\n\n Returns:\n A list of the form ``components[i] = (smiles_i, count_i)`` where\n ``smiles_i`` is the SMILES representation of component `i` and\n ``count_i`` is the number of corresponding instances of that component\n to create.\n \"\"\"\n components = []\n\n current_smiles = None\n current_count = 0\n\n for molecule in topology.molecules:\n smiles = molecule.to_smiles()\n\n if smiles == current_smiles:\n current_count += 1\n continue\n\n if current_count > 0:\n components.append((current_smiles, current_count))\n\n current_smiles = smiles\n current_count = 1\n\n if current_count > 0:\n components.append((current_smiles, current_count))\n\n return components\n
"},{"location":"reference/utils/topology/#absolv.utils.topology.topology_to_atom_indices","title":"topology_to_atom_indices","text":"topology_to_atom_indices(\n topology: Topology,\n) -> list[set[int]]\n
A helper method for extracting the sets of atom indices associated with each molecule in a topology.
Parameters:
-
topology (Topology) \u2013 The topology to extract the atom indices from.
Returns:
-
list[set[int]] \u2013 The set of atoms indices associated with each molecule in the topology.
Source code in absolv/utils/topology.py def topology_to_atom_indices(topology: openff.toolkit.Topology) -> list[set[int]]:\n \"\"\"A helper method for extracting the sets of atom indices associated with each\n molecule in a topology.\n\n Args:\n topology: The topology to extract the atom indices from.\n\n Returns:\n The set of atoms indices associated with each molecule in the topology.\n \"\"\"\n\n atom_indices: list[set[int]] = []\n current_atom_idx = 0\n\n for molecule in topology.molecules:\n atom_indices.append({i + current_atom_idx for i in range(molecule.n_atoms)})\n current_atom_idx += molecule.n_atoms\n\n return atom_indices\n
"},{"location":"user-guide/overview/","title":"Overview","text":"To begin with, we define a configuration that encodes the entirety of the absolute free energy calculation.
This includes specifying the solutes and the two solvents that they will be transferred between:
import absolv.config\n\nsystem=absolv.config.System(\n solutes={\"CCO\": 1}, solvent_a=None, solvent_b={\"O\": 895}\n)\n
Each key should be a SMILES representation of a molecule, and the value should be the number of copies of that molecule to include in the system. There may be multiple in the case of, e.g., ion pairs like Na+ and Cl-. None may be used to specify vacuum, e.g., in the above case that the solute will be transferred from vacuum into bulk water.
The temperature and pressure that the calculation will be performed at must also be specified:
import openmm.unit\n\ntemperature=298.15 * openmm.unit.kelvin,\npressure=1.0 * openmm.unit.atmosphere\n
Finally, the alchemical pathway to transform the solute along in each solvent must be specified. This can either be a more traditional 'equilibrium' pathway, or a 'non-equilibrium' pathway:
EquilibriumNon-equilibrium import absolv.config\n\nalchemical_protocol_a=absolv.config.EquilibriumProtocol(\n lambda_sterics=[1.0, 1.0, 1.0, 1.0, 1.0],\n lambda_electrostatics=[1.0, 0.75, 0.5, 0.25, 0.0]\n)\nalchemical_protocol_b=absolv.config.EquilibriumProtocol(\n lambda_sterics=[\n 1.00, 1.00, 1.00, 1.00, 1.00, 0.95, 0.90, 0.80, 0.70, 0.60, 0.50, 0.40,\n 0.35, 0.30, 0.25, 0.20, 0.15, 0.10, 0.05, 0.00,\n ],\n lambda_electrostatics=[\n 1.00, 0.75, 0.50, 0.25, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00,\n 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00,\n ]\n)\n
Here the default lambda schedule from FreeSolv is used. A lambda of 1.0 indicates a fully interacting system, and a lambda of 0.0 indicates a fully decoupled system.
import femto.md.config\nimport openmm.unit\n\nimport absolv.config\n\nintegrator = femto.md.config.LangevinIntegrator(\n timestep=2.0 * openmm.unit.femtoseconds\n)\n\nalchemical_protocol_a=absolv.config.NonEquilibriumProtocol(\n # the protocol to use for the production run at each end state\n production_protocol=absolv.config.SimulationProtocol(\n integrator=integrator, n_steps=6250 * 160\n ),\n # the interval with which to store frames that NEQ switching\n # simulations will be launched from\n production_report_interval=6250,\n # define how the NEQ switching will be performed\n switching_protocol=absolv.config.SwitchingProtocol(\n n_electrostatic_steps=60,\n n_steps_per_electrostatic_step=100,\n # intra-molecular vdW interactions are not decoupled by default, so we\n # don't need to do any vdW decoupling in vacuum when there's only one solvent\n n_steric_steps=0,\n n_steps_per_steric_step=0\n )\n)\nalchemical_protocol_b=absolv.config.NonEquilibriumProtocol(\n production_protocol=absolv.config.SimulationProtocol(\n integrator=integrator, n_steps=6250 * 160,\n ),\n switching_protocol=absolv.config.SwitchingProtocol(\n # Annihilate the electrostatic interactions over the first 12 ps\n n_electrostatic_steps=60,\n n_steps_per_electrostatic_step=100,\n # followed by decoupling the vdW interactions over the next 38 ps\n n_steric_steps=190,\n n_steps_per_steric_step=100,\n )\n)\n
These individual components are then combined into a single configuration object:
import absolv.config\n\nconfig = absolv.config.Config(\n temperature=temperature,\n pressure=pressure,\n alchemical_protocol_a=alchemical_protocol_a,\n alchemical_protocol_b=alchemical_protocol_b,\n)\n
which can be used to trivially set up
import openff.toolkit\nforce_field = openff.toolkit.ForceField(\"openff-2.1.0.offxml\")\n\nimport absolv.runner\nprepared_system_a, prepared_system_b = absolv.runner.setup(system, config, force_field)\n
and run the calculation:
EquilibriumNon-equilibrium result = absolv.runner.run_eq(\n config, prepared_system_a, prepared_system_b, \"CUDA\"\n)\nprint(result)\n
result = absolv.runner.run_neq(\n config, prepared_system_a, prepared_system_b, \"CUDA\"\n)\nprint(result)\n
where the result will be a Result object.
"},{"location":"user-guide/theory/","title":"Theory","text":"The change in free energy of transferring a solute from one solvent to another \\(\\Delta G_{A->B}\\), can be readily computed by constructing a thermodynamic cycle in which the solute in each solvent is alchemically transformed into a non-interacting state as shown below:
Namely,
\\[ \\Delta G_{A->B} = \\Delta G_1 - \\Delta G_3 \\] where here \\(\\Delta G_1\\) and \\(\\Delta G_3\\) are the change in free energy of alchemically transforming a solute (or perhaps multiple solutes in the case of charged molecules with counter ions) so that it no longer interacts with the surrounding corresponding 'solvent'.
Computing the solvation (/ hydration) free energy is a special case of this cycle when the first 'solvent' is vacuum.
This framework currently offers two routes to computing such a change in free energy as the solute is being alchemically transformed, a more commonly used 'equilibrium route' and a 'non-equilibrium' route.
Note For more information about how the different inter- and intramolecular interactions are alchemically modified see the transformations page.
"},{"location":"user-guide/theory/#equilibrium-calculations","title":"Equilibrium Calculations","text":"Within this framework we refer to free energy calculations that involve dividing the alchemical pathway into discrete windows at each value of the coupling parameter \\(\\lambda\\), collecting equilibrium samples for each such discrete state, and using these samples to compute the free energy using an approach such as thermodynamic integration (TI) 1, BAR 2, and MBAR 3 as 'equilibrium' free energy calculations.
At present absolv does not offer functionality for computing the derivatives with respect to lambda required for TI, and only supports MBAR and technically BAR although this estimator is not recommended.
See the overview for more information on running equilibrium free energy calculations using absolv.
"},{"location":"user-guide/theory/#non-equilibrium-calculations","title":"Non-equilibrium Calculations","text":"Within this framework we refer to free energy calculations that proceed by generating equilibrium configurations at both end states (with the solute fully interacting with the solute and with the solute-solvent interactions fully decoupled / annihilated), and then driving each configuration non-reverisbly along the alchemical pathway by scaling the coupling factor \\(\\lambda\\) as a function of time 4.
From a practical perspective it is computationally more efficient and convenient to proceed along the alchemical pathway in a stepwise, rather than continuous fashion. More specifically, the protocol proceeds by making a perturbation to \\(\\lambda\\), followed by several relaxation steps, and repeating these two sub-stebs until the alchemical transformation is complete.
in this way the required to transform the solute from the interacting to the non-interacting state and likewise from the non-interacting to the interacting state can be computed according to
\\[W = \\sum_{i=0}^{n-1} \\left[u_{\\lambda_{i+1}}\\left(x_i\\right) - u_{\\lambda_{i}}\\left(x_i\\right)\\right]\\] where here \\(u_{\\lambda_i}\\left(x_i\\right)\\) is the reduced potential evaluated at configuration \\(i\\) and \\(\\lambda_i\\).
The free energy is then estimated by solving
\\[\\sum^N_{i=1}\\dfrac{1}{1+\\exp(\\beta(W^f_i-\\Delta F))} = \\sum^N_{i=1}\\dfrac{1}{1+\\exp(\\beta(W^r_i+\\Delta F))}\\] self consistently where \\(W^f_i\\) corresponds to work computed along the forward pathway going from the interacting to the non-interacting state and \\(W^r_i\\) to work computed along the reverse pathway going from the non-interacting to the interacting state. \\(N\\) refers to the total number of equilibrium snapshots that were generated at each end state.
-
TP Straatsma and JA McCammon. Computational alchemy. Annual Review of Physical Chemistry, 43(1):407\u2013435, 1992.\u00a0\u21a9
-
Charles H Bennett. Efficient estimation of free energy differences from monte carlo data. Journal of Computational Physics, 22(2):245\u2013268, 1976.\u00a0\u21a9
-
Michael R Shirts and John D Chodera. Statistically optimal analysis of samples from multiple equilibrium states. The Journal of chemical physics, 129(12):124105, 2008.\u00a0\u21a9
-
Andrew J Ballard and Christopher Jarzynski. Replica exchange with nonequilibrium switches: enhancing equilibrium sampling by increasing replica overlap. The Journal of chemical physics, 136(19):194101, 2012.\u00a0\u21a9
"},{"location":"user-guide/transformations/","title":"Transformations","text":"The absolv framework supports alchemically transforming the electrostatic and vdW interactions within an OpenMM System object via the absolv.fep.apply_fep function.
Due to the huge flexibility of the OpenMM system, allowing completely custom electrostatic and vdW forces to be implemented, it would be impossible to support all possible combinations of non-bonded forces. As such, the framework currently only systems that contain
- vdW + electrostatic interactions in a single built-in
NonbondedForce object - electrostatic interactions in a built-in
NonbondedForce object, vdW interactions in a CustomNonbondedForce object and vdW 1-4 interactions in a CustomBondForce - the above combinations sans any electrostatics interactions
are supported.
"},{"location":"user-guide/transformations/#electrostatics","title":"Electrostatics","text":"The electrostatic interactions in the system are alchemically transformed by linearly scaling all partial charges on particles in solute molecules by \\(\\lambda_{elec}\\), corresponding to a \"lambda_electrostatics\" context variable that will be added by the factory, such that
\\[q^{sol}_i \\left(\\lambda_{elec}\\right) = \\lambda_{elec} \\times q^{sol}_i\\] All intramolecular interactions will be switched off during the alchemical transformation. This is referred to as 'annihilating' the electrostatic interactions in other packages and some literature.
Because the charges are scaled directly, the energy contributions of the alchemically scaled electrostatic interactions will be
\\[U^E = \\lambda_{elec} U^E_{sol-solv} + \\lambda_{elec}^2 U^E_{sol-sol} + U^E_{solv-solv}\\] where \\(U^E_{sol-sol}\\), \\(U^E_{sol-solv}\\) and \\(U^E_{solv-solv}\\) are the un-scaled electrostatic contributions to the energies of the solute-solute, solute-solvent and solvent-solvent interactions respectively.
"},{"location":"user-guide/transformations/#vdw","title":"vdW","text":"Currently vdW interactions can only be transformed if they are stored in a standard NonbondedForce OR if they are split between a CustomBondForce (1-2, 1-3, and 1-4 interactions) and a CustomNonbondedForce (the remaining 1-n and intermolecular interactions).
The interactions will be transformed according to \\(\\lambda_{vdW}\\) which corresponds to a \"lambda_sterics\" context variable that will be added by the factory.
Only intermolecular vdW interactions will be alchemically scaled, while all intramolecular interactions will be left un-scaled. This is is referred to as 'decoupling' the vdW interactions in other packages and some literature.
"},{"location":"user-guide/transformations/#lennard-jones","title":"Lennard--Jones","text":"If the vdW interactions are stored in a standard NonbondedForce, then the alchemical factory will split them so that
-
the NonbondedForce force retains all interactions between solvent particles
-
all intermolecular alchemical (both solute-solvent and solute-solute) interactions are moved into a new CustomNonbondedForce
-
all solute intramolecular interactions are moved into a new CustomNonbondedForce
Note The intramolecular solute-solute interactions won't use any periodic boundary corrections such that the the decoupled state of the solute corresponds to the proper vacuum state without periodicity effects.
The CustomNonbondedForce will copy over all settings (including cut-off, long-range correction, etc) from the original NonbondedForce, but will replace the normal LJ energy expression with a soft-core version. By default, this takes the form:
\\[U^{vdW} \\left( \\lambda_{vdW} \\right) = \\lambda_{vdW} \\times 4 \\varepsilon \\left[ \\left( \\dfrac{\\sigma}{r_{eff}}\\right)^{12} - \\left( \\dfrac{\\sigma}{r_{eff}}\\right)^{6} \\right]\\] where
\\[r_{eff} = \\sigma \\left( \\dfrac{1}{2} \\left(1 - \\lambda_{vdW}\\right) + \\left( \\dfrac{r}{\\sigma} \\right) ^ 6 \\right) ^ \\frac{1}{6}\\]"},{"location":"user-guide/transformations/#custom-vdw-forms","title":"Custom vdW Forms","text":"If the vdW interactions are split across a CustomNonbondedForce and a CustomBondForce then the alchemical factory will further split them so that
-
the original CustomNonbondedForce force will retain all interactions between solvent particles
-
all solute intramolecular interactions are moved into a new CustomNonbondedForce
Note The intramolecular solute-solute interactions won't use any periodic boundary corrections such that the the decoupled state of the solute corresponds to the proper vacuum state without periodicity effects.
-
all intermolecular alchemical (both solute-solvent and solute-solute) interactions are moved into another new CustomNonbondedForce with a modified energy expression such that
\\[U^{vdW}_{sol-solv} \\left( \\lambda_{vdW} \\right) = \\lambda_{vdW} \\times U^{vdW}_{sol-solv}\\]
"}]}
\ No newline at end of file
+{"config":{"lang":["en"],"separator":"[\\s\\-]+","pipeline":["stopWordFilter"]},"docs":[{"location":"","title":"Overview","text":"ABsolute SOLVantion Free Energy Calculations Absolute solvation free energy calculations using OpenMM
The absolv framework aims to offer a simple API for computing the change in free energy when transferring a solute from one solvent to another, or to vacuum in the case of solvation free energy calculations.
It offers two routes to this end: standard equilibrium calculations and non-equilibrium switching type calculations, where the latter will be the main focus of this framework.
Warning
This code is currently experimental and under active development. If you are using this it, please be aware that it is not guaranteed to provide correct results, the documentation and testing may be incomplete, and the API may change without notice.
"},{"location":"#installation","title":"Installation","text":"This package can be installed using conda (or mamba, a faster version of conda):
mamba install -c conda-forge absolv\n
If you are running with MPI on an HPC cluster, you may need to instruct conda to use your local installation depending on your setup
mamba install -c conda-forge absolv \"openmpi=4.1.5=*external*\"\n
"},{"location":"#getting-started","title":"Getting Started","text":"To get started, see the usage guide.
"},{"location":"development/","title":"Development","text":"To create a development environment, you must have mamba installed.
A development conda environment can be created and activated with:
make env\nconda activate absolv\n
To format the codebase:
make format\n
To run the unit tests:
make test\n
To serve the documentation locally:
mkdocs serve\n
"},{"location":"reproducibility/","title":"Reproducibility","text":"While every effort has been made to ensure the 'correctness and reproducibility' of any results computed using this framework, achieving consistent free energies between different frameworks and simulation engines has been notoriously tricky 1.
In an attempt to ensure that this framework remains at least self-consistent between versions, and as consistent as possible with other packages, a suite of regression tests are provided in the main GitHub repository.
These include tests to ensure that:
- the energies of systems that been alchemically modified are consistent with the original system
- systems that contain custom non-bonded force are alchemically transformed in the exact same was as systems that contain the built-in OpenMM LJ non-bonded force
and most importantly, that computing the free energies using both the 'equilibrium' and 'non-equilibrium' methods supported in this framework are in agreement amongst themselves, with values computed using Yank, and with the GROMACS values reported by Loeffler et al 1.
"},{"location":"reproducibility/#regression-results","title":"Regression Results","text":"The results of running the free energy comparisons using the latest version of the framework are shown below:
Here the error bars show the standard deviation computed across three replicas.
-
Hannes H Loeffler, Stefano Bosisio, Guilherme Duarte Ramos Matos, Donghyuk Suh, Benoit Roux, David L Mobley, and Julien Michel. Reproducibility of free energy calculations across different molecular simulation software packages. Journal of chemical theory and computation, 14(11):5567\u20135582, 2018.\u00a0\u21a9\u21a9
"},{"location":"reference/","title":"Index","text":""},{"location":"reference/#absolv","title":"absolv","text":"Absolute solvation free energy calculations using OpenMM
Modules:
-
config \u2013 Configure free energy calculations.
-
fep \u2013 Prepare OpenMM systems for FEP calculations.
-
neq \u2013 Run non-equilibrium forward and reverse sampling.
-
runner \u2013 Run calculations defined by a config.
-
setup \u2013 Setup systems ready for calculations.
-
tests \u2013 -
utils \u2013 Common utils
"},{"location":"reference/SUMMARY/","title":"SUMMARY","text":" - absolv
- config
- fep
- neq
- runner
- setup
- utils
- openmm
- topology
"},{"location":"reference/config/","title":"Configure free energy calculations.
Classes:
-
System \u2013 Define the two solvents that solutes will be transferred between (a -> b),
-
MinimizationProtocol \u2013 Configure how a system should be energy minimized.
-
SimulationProtocol \u2013 Configure how a system should be evolved by molecular simulation.
-
HREMDProtocol \u2013 Configure how a system should be evolved by Hamiltonian replica exchange.
-
EquilibriumProtocol \u2013 Configure how an equilibrium (e.g. TI, MBAR) alchemical free energy
-
SwitchingProtocol \u2013 Configure non-reversibly driving a system between to alchemical states.
-
NonEquilibriumProtocol \u2013 Configure a non-equilibrium alchemical free energy calculation [1, 2].
-
Config \u2013 Configure a transfer free energy calculation.
-
Result \u2013 The result of a free energy calculation.
"},{"location":"reference/config/#absolv.config.System","title":"System pydantic-model","text":" Bases: BaseModel
Define the two solvents that solutes will be transferred between (a -> b), as well as the solutes themselves.
Fields:
-
solutes (dict[str, PositiveInt]) -
solvent_a (dict[str, PositiveInt] | None) -
solvent_b (dict[str, PositiveInt] | None) -
n_solute_molecules (int) -
n_solvent_molecules_a (int) -
n_solvent_molecules_b (int) -
components_a (list[tuple[str, int]]) -
components_b (list[tuple[str, int]])
Validators:
-
_validate_solvent_a \u2192 solvent_a -
_validate_solvent_b \u2192 solvent_b
"},{"location":"reference/config/#absolv.config.System.solutes","title":"solutes pydantic-field","text":"solutes: dict[str, PositiveInt]\n
A dictionary containing the SMILES patterns of each solute in the system as well as how many instances of each there should be.
"},{"location":"reference/config/#absolv.config.System.solvent_a","title":"solvent_a pydantic-field","text":"solvent_a: dict[str, PositiveInt] | None\n
A dictionary containing the SMILES patterns of each component in the first solvent as well as how many instances of each there should be.A value of None should be used to indicate vacuum.
"},{"location":"reference/config/#absolv.config.System.solvent_b","title":"solvent_b pydantic-field","text":"solvent_b: dict[str, PositiveInt] | None\n
A dictionary containing the SMILES patterns of each component in the second solvent as well as how many instances of each there should be. A value of None should be used to indicate vacuum.
"},{"location":"reference/config/#absolv.config.System.n_solute_molecules","title":"n_solute_molecules pydantic-field","text":"n_solute_molecules: int\n
Returns the total number of solute molecules that will be present.
"},{"location":"reference/config/#absolv.config.System.n_solvent_molecules_a","title":"n_solvent_molecules_a pydantic-field","text":"n_solvent_molecules_a: int\n
Returns the total number of solvent molecules that will be present in the first solution.
"},{"location":"reference/config/#absolv.config.System.n_solvent_molecules_b","title":"n_solvent_molecules_b pydantic-field","text":"n_solvent_molecules_b: int\n
Returns the total number of solvent molecules that will be present in the second solution.
"},{"location":"reference/config/#absolv.config.System.components_a","title":"components_a pydantic-field","text":"components_a: list[tuple[str, int]]\n
Returns the identities and counts of the molecules present in the first system.
Returns:
-
list[tuple[str, int]] \u2013 The SMILES representation and count of each component.
"},{"location":"reference/config/#absolv.config.System.components_b","title":"components_b pydantic-field","text":"components_b: list[tuple[str, int]]\n
Returns the identities and counts of the molecules present in the second system.
Returns:
-
list[tuple[str, int]] \u2013 The SMILES representation and count of each component.
"},{"location":"reference/config/#absolv.config.MinimizationProtocol","title":"MinimizationProtocol pydantic-model","text":" Bases: BaseModel
Configure how a system should be energy minimized.
Fields:
-
tolerance (OpenMMQuantity[kilojoule_per_mole / nanometers]) -
max_iterations (NonNegativeInt)
"},{"location":"reference/config/#absolv.config.MinimizationProtocol.tolerance","title":"tolerance pydantic-field","text":"tolerance: OpenMMQuantity[\n kilojoule_per_mole / nanometers\n] = (10.0 * kilojoule_per_mole / nanometers)\n
How precisely the energy minimum must be located [kj / mol / nm]
"},{"location":"reference/config/#absolv.config.MinimizationProtocol.max_iterations","title":"max_iterations pydantic-field","text":"max_iterations: NonNegativeInt = 0\n
The maximum number of iterations to perform. If this is 0, minimization is continued until the results converge.
"},{"location":"reference/config/#absolv.config.SimulationProtocol","title":"SimulationProtocol pydantic-model","text":" Bases: BaseModel
Configure how a system should be evolved by molecular simulation.
Fields:
-
integrator (LangevinIntegrator) -
n_steps (PositiveInt)
"},{"location":"reference/config/#absolv.config.SimulationProtocol.integrator","title":"integrator pydantic-field","text":"integrator: LangevinIntegrator = LangevinIntegrator(\n timestep=2.0 * femtosecond, friction=1.0 / picosecond\n)\n
The integrator to use for the simulation.
"},{"location":"reference/config/#absolv.config.SimulationProtocol.n_steps","title":"n_steps pydantic-field","text":"n_steps: PositiveInt\n
The number of steps to evolve the system by.
"},{"location":"reference/config/#absolv.config.HREMDProtocol","title":"HREMDProtocol pydantic-model","text":" Bases: BaseModel
Configure how a system should be evolved by Hamiltonian replica exchange.
Fields:
-
integrator (LangevinIntegrator) -
n_warmup_steps (int) -
n_steps_per_cycle (int) -
n_cycles (int)
"},{"location":"reference/config/#absolv.config.HREMDProtocol.integrator","title":"integrator pydantic-field","text":"integrator: LangevinIntegrator = LangevinIntegrator(\n timestep=2.0 * femtosecond, friction=1.0 / picosecond\n)\n
The integrator to use for the simulation.
"},{"location":"reference/config/#absolv.config.HREMDProtocol.n_warmup_steps","title":"n_warmup_steps pydantic-field","text":"n_warmup_steps: int = 50000\n
The number of steps to run each replica for before starting hremd trials. All energies gathered during this period will be discarded.
"},{"location":"reference/config/#absolv.config.HREMDProtocol.n_steps_per_cycle","title":"n_steps_per_cycle pydantic-field","text":"n_steps_per_cycle: int = 1000\n
The number of steps to propagate the system by before attempting an exchange.
"},{"location":"reference/config/#absolv.config.HREMDProtocol.n_cycles","title":"n_cycles pydantic-field","text":"n_cycles: int = 2500\n
The number of cycles of 'propagate the system' -> 'exchange replicas' to run.
"},{"location":"reference/config/#absolv.config.EquilibriumProtocol","title":"EquilibriumProtocol pydantic-model","text":" Bases: BaseModel
Configure how an equilibrium (e.g. TI, MBAR) alchemical free energy calculation.
Fields:
-
type (Literal['equilibrium']) -
minimization_protocol (MinimizationProtocol) -
equilibration_protocol (SimulationProtocol) -
production_protocol (HREMDProtocol) -
lambda_sterics (list[confloat(ge=0.0, le=1.0)]) -
lambda_electrostatics (list[confloat(ge=0.0, le=1.0)]) -
n_states (int)
"},{"location":"reference/config/#absolv.config.EquilibriumProtocol.minimization_protocol","title":"minimization_protocol pydantic-field","text":"minimization_protocol: MinimizationProtocol = (\n MinimizationProtocol()\n)\n
Whether to minimize the energy of the system prior to any simulations.
"},{"location":"reference/config/#absolv.config.EquilibriumProtocol.equilibration_protocol","title":"equilibration_protocol pydantic-field","text":"equilibration_protocol: SimulationProtocol = (\n SimulationProtocol(n_steps=100000)\n)\n
The (optional) protocol that describes the equilibration simulation to run prior to the production one.
"},{"location":"reference/config/#absolv.config.EquilibriumProtocol.production_protocol","title":"production_protocol pydantic-field","text":"production_protocol: HREMDProtocol = HREMDProtocol(\n n_steps_per_cycle=6250, n_cycles=160\n)\n
The protocol that describes the production to run.
"},{"location":"reference/config/#absolv.config.EquilibriumProtocol.lambda_sterics","title":"lambda_sterics pydantic-field","text":"lambda_sterics: list[confloat(ge=0.0, le=1.0)]\n
The alchemical pathway to transform the vdW interactions along. A value of 1.0 represents a fully interacting system while a value of 0.0 represents a system with the solute-solute and solute-solvent vdW interactions disabled.
"},{"location":"reference/config/#absolv.config.EquilibriumProtocol.lambda_electrostatics","title":"lambda_electrostatics pydantic-field","text":"lambda_electrostatics: list[confloat(ge=0.0, le=1.0)]\n
The alchemical pathway to transform the electrostatic interactions along. A value of 1.0 represents a fully interacting system while a value of 0.0 represents a system with the solute-solute and solute-solvent electrostatic interactions disabled.
"},{"location":"reference/config/#absolv.config.EquilibriumProtocol.n_states","title":"n_states pydantic-field","text":"n_states: int\n
Returns the number of lambda states that will be sampled at.
"},{"location":"reference/config/#absolv.config.SwitchingProtocol","title":"SwitchingProtocol pydantic-model","text":" Bases: BaseModel
Configure non-reversibly driving a system between to alchemical states.
Fields:
-
n_electrostatic_steps (NonNegativeInt) -
n_steps_per_electrostatic_step (NonNegativeInt) -
n_steric_steps (NonNegativeInt) -
n_steps_per_steric_step (NonNegativeInt)
"},{"location":"reference/config/#absolv.config.SwitchingProtocol.n_electrostatic_steps","title":"n_electrostatic_steps pydantic-field","text":"n_electrostatic_steps: NonNegativeInt\n
The number of steps to annihilate the electrostatics interactions over. The total time needed to annihilate the electrostatics interactions will be n_electrostatic_steps * n_steps_per_electrostatic_step * timestep
"},{"location":"reference/config/#absolv.config.SwitchingProtocol.n_steps_per_electrostatic_step","title":"n_steps_per_electrostatic_step pydantic-field","text":"n_steps_per_electrostatic_step: NonNegativeInt\n
The number of timesteps to evolve the system by each time the electrostatics interactions are modified. A value of 1 will give a 'smooth' transition between the each discrete lambda value whereas a value greater than 1 will yield a stepwise transition as shown in Figure 3 of doi:10.1063/1.4712028.
"},{"location":"reference/config/#absolv.config.SwitchingProtocol.n_steric_steps","title":"n_steric_steps pydantic-field","text":"n_steric_steps: NonNegativeInt\n
The number of steps to decouple the sterics interactions over once the electrostatics interactions have been annihilated. The total time needed to annihilate the sterics interactions will be n_steric_steps * n_steps_per_steric_step * timestep
"},{"location":"reference/config/#absolv.config.SwitchingProtocol.n_steps_per_steric_step","title":"n_steps_per_steric_step pydantic-field","text":"n_steps_per_steric_step: NonNegativeInt\n
The number of timesteps to evolve the system by each time the sterics interactions are modified. A value of 1 will give a 'smooth' transition between the each discrete lambda value whereas a value greater than 1 will yield a stepwise transition as shown in Figure 3 of doi:10.1063/1.4712028.
"},{"location":"reference/config/#absolv.config.NonEquilibriumProtocol","title":"NonEquilibriumProtocol pydantic-model","text":" Bases: BaseModel
Configure a non-equilibrium alchemical free energy calculation [1, 2].
It is expected that first the electrostatics interactions will be annihilated followed by a decoupling of the sterics interactions.
References [1] Ballard, Andrew J., and Christopher Jarzynski. \"Replica exchange with nonequilibrium switches: Enhancing equilibrium sampling by increasing replica overlap.\" The Journal of chemical physics 136.19 (2012): 194101.
[2] Gapsys, Vytautas, et al. \"Large scale relative protein ligand binding affinities using non-equilibrium alchemy.\" Chemical Science 11.4 (2020): 1140-1152.
Fields:
-
type (Literal['non-equilibrium']) -
minimization_protocol (MinimizationProtocol | None) -
equilibration_protocol (SimulationProtocol | None) -
production_protocol (SimulationProtocol) -
production_report_interval (PositiveInt) -
switching_protocol (SwitchingProtocol)
"},{"location":"reference/config/#absolv.config.NonEquilibriumProtocol.minimization_protocol","title":"minimization_protocol pydantic-field","text":"minimization_protocol: MinimizationProtocol | None = (\n MinimizationProtocol()\n)\n
The (optional) protocol to follow when minimizing the system in both the end states prior to running the equilibrium simulations.
"},{"location":"reference/config/#absolv.config.NonEquilibriumProtocol.equilibration_protocol","title":"equilibration_protocol pydantic-field","text":"equilibration_protocol: SimulationProtocol | None = (\n SimulationProtocol(n_steps=100000)\n)\n
The (optional) protocol to follow when equilibrating the system in both the end states prior to running the production equilibrium simulations.
"},{"location":"reference/config/#absolv.config.NonEquilibriumProtocol.production_protocol","title":"production_protocol pydantic-field","text":"production_protocol: SimulationProtocol = (\n SimulationProtocol(n_steps=6250 * 160)\n)\n
The protocol to follow when running the production equilibrium simulation in both the end states. The snapshots generated at the end of each iteration will be used for each non-equilibrium switch.
"},{"location":"reference/config/#absolv.config.NonEquilibriumProtocol.production_report_interval","title":"production_report_interval pydantic-field","text":"production_report_interval: PositiveInt = 6250\n
The interval at which to write out the simulation state during the production simulation.
"},{"location":"reference/config/#absolv.config.NonEquilibriumProtocol.switching_protocol","title":"switching_protocol pydantic-field","text":"switching_protocol: SwitchingProtocol\n
The protocol that describes how each snapshot generated during the production simulation should be driven from state 0 -> 1 and likewise 1 -> 0 in order to compute the non-equilibrium work distributions.
"},{"location":"reference/config/#absolv.config.Config","title":"Config pydantic-model","text":" Bases: BaseModel
Configure a transfer free energy calculation.
Fields:
-
temperature (OpenMMQuantity[kelvin]) -
pressure (OpenMMQuantity[atmosphere] | None) -
alchemical_protocol_a (AlchemicalProtocol) -
alchemical_protocol_b (AlchemicalProtocol)
"},{"location":"reference/config/#absolv.config.Config.temperature","title":"temperature pydantic-field","text":"temperature: OpenMMQuantity[kelvin]\n
The temperature to calculate at [K].
"},{"location":"reference/config/#absolv.config.Config.pressure","title":"pressure pydantic-field","text":"pressure: OpenMMQuantity[atmosphere] | None\n
The pressure to calculate at [atm].
"},{"location":"reference/config/#absolv.config.Config.alchemical_protocol_a","title":"alchemical_protocol_a pydantic-field","text":"alchemical_protocol_a: AlchemicalProtocol\n
The protocol that describes the alchemical pathway to transform the solute along in the first solvent.
"},{"location":"reference/config/#absolv.config.Config.alchemical_protocol_b","title":"alchemical_protocol_b pydantic-field","text":"alchemical_protocol_b: AlchemicalProtocol\n
The protocol that describes the alchemical pathway to transform the solute along in the second solvent.
"},{"location":"reference/config/#absolv.config.Result","title":"Result pydantic-model","text":" Bases: BaseModel
The result of a free energy calculation.
Fields:
-
dg_solvent_a (OpenMMQuantity[KCAL_MOL]) -
dg_std_solvent_a (OpenMMQuantity[KCAL_MOL]) -
dg_solvent_b (OpenMMQuantity[KCAL_MOL]) -
dg_std_solvent_b (OpenMMQuantity[KCAL_MOL]) -
dg (Quantity) -
dg_std (Quantity)
"},{"location":"reference/config/#absolv.config.Result.dg_solvent_a","title":"dg_solvent_a pydantic-field","text":"dg_solvent_a: OpenMMQuantity[KCAL_MOL]\n
The change in free energy of alchemically transforming the solute from an interacting to a non-interacting state in the first solvent.
"},{"location":"reference/config/#absolv.config.Result.dg_std_solvent_a","title":"dg_std_solvent_a pydantic-field","text":"dg_std_solvent_a: OpenMMQuantity[KCAL_MOL]\n
The standard error in dg_solvent_a.
"},{"location":"reference/config/#absolv.config.Result.dg_solvent_b","title":"dg_solvent_b pydantic-field","text":"dg_solvent_b: OpenMMQuantity[KCAL_MOL]\n
The change in free energy of alchemically transforming the solute from an interacting to a non-interacting state in the second solvent.
"},{"location":"reference/config/#absolv.config.Result.dg_std_solvent_b","title":"dg_std_solvent_b pydantic-field","text":"dg_std_solvent_b: OpenMMQuantity[KCAL_MOL]\n
The standard error in dg_solvent_b.
"},{"location":"reference/config/#absolv.config.Result.dg","title":"dg pydantic-field","text":"dg: Quantity\n
The change in free energy of transferring the solute from solvent-a to solvent-b in units of kT.
"},{"location":"reference/config/#absolv.config.Result.dg_std","title":"dg_std pydantic-field","text":"dg_std: Quantity\n
The standard error in ddg.
"},{"location":"reference/fep/","title":"Prepare OpenMM systems for FEP calculations.
Functions:
-
apply_fep \u2013 Generate a system whereby a number of the molecules can be alchemically
-
set_fep_lambdas \u2013 Set the values of the alchemical lambdas on an OpenMM context.
"},{"location":"reference/fep/#absolv.fep.apply_fep","title":"apply_fep","text":"apply_fep(\n system: System,\n alchemical_indices: list[set[int]],\n persistent_indices: list[set[int]],\n custom_alchemical_potential: str | None = None,\n) -> System\n
Generate a system whereby a number of the molecules can be alchemically transformed from a base chemical system.
Notes -
Currently only OpenMM systems that have:
-
vdW + electrostatics in a single built-in non-bonded force
- electrostatics in a built-in non-bonded force, vdW in a custom non-bonded force and vdW 1-4 interactions in a custom bond force
- all of the above sans any electrostatics
are supported. * By default a soft-core version of the LJ potential with a-b-c of 1-1-6 and alpha=0.5 that can be scaled by a global lambda_sterics parameter will be used for alchemical-chemical vdW interactions embedded in an OpenMM NonbondedForce while the energy expression set on a CustomNonbondedForce will be be modified to have the form \"lambda_sterics*({original_expression})\".
Parameters:
-
system (System) \u2013 The chemical system to generate the alchemical system from
-
alchemical_indices (list[set[int]]) \u2013 The atom indices corresponding to each molecule that should be alchemically transformable. The atom indices must correspond to all atoms in each molecule as alchemically transforming part of a molecule is not supported.
-
persistent_indices (list[set[int]]) \u2013 The atom indices corresponding to each molecule that should not be alchemically transformable.
-
custom_alchemical_potential (str | None, default: None ) \u2013 A custom expression to use for the potential energy function that describes the chemical-alchemical intermolecular interactions. See the Notes for information about the default value.
The expression must include \"lambda_sterics\".
Source code in absolv/fep.py def apply_fep(\n system: openmm.System,\n alchemical_indices: list[set[int]],\n persistent_indices: list[set[int]],\n custom_alchemical_potential: str | None = None,\n) -> openmm.System:\n \"\"\"Generate a system whereby a number of the molecules can be alchemically\n transformed from a base chemical system.\n\n Notes:\n * Currently only OpenMM systems that have:\n\n - vdW + electrostatics in a single built-in non-bonded force\n - electrostatics in a built-in non-bonded force, vdW in a custom\n **non-bonded** force and vdW 1-4 interactions in a custom **bond** force\n - all of the above sans any electrostatics\n\n are supported.\n * By default a soft-core version of the LJ potential with a-b-c of 1-1-6\n and alpha=0.5 that can be scaled by a global `lambda_sterics` parameter\n will be used for alchemical-chemical vdW interactions embedded in an\n OpenMM ``NonbondedForce`` while the energy expression set on a\n ``CustomNonbondedForce`` will be be modified to have the form\n ``\"lambda_sterics*({original_expression})\"``.\n\n Args:\n system: The chemical system to generate the alchemical system from\n alchemical_indices: The atom indices corresponding to each molecule that\n should be alchemically transformable. The atom indices **must**\n correspond to **all** atoms in each molecule as alchemically\n transforming part of a molecule is not supported.\n persistent_indices: The atom indices corresponding to each molecule that\n should **not** be alchemically transformable.\n custom_alchemical_potential: A custom expression to use for the potential\n energy function that describes the chemical-alchemical intermolecular\n interactions. See the Notes for information about the default value.\n\n The expression **must** include ``\"lambda_sterics\"``.\n \"\"\"\n\n system = copy.deepcopy(system)\n\n # Make sure we track v-sites attached to any solutes that may be alchemically\n # turned off. We do this as a post-process step as the OpenFF toolkit does not\n # currently expose a clean way to access this information.\n atom_indices = alchemical_indices + persistent_indices\n atom_indices = _find_v_sites(system, atom_indices)\n\n alchemical_indices = atom_indices[: len(alchemical_indices)]\n persistent_indices = atom_indices[len(alchemical_indices) :]\n\n (\n nonbonded_force,\n custom_nonbonded_force,\n custom_bond_force,\n ) = _find_nonbonded_forces(system)\n\n if nonbonded_force is not None:\n _add_electrostatics_lambda(nonbonded_force, alchemical_indices)\n\n if custom_nonbonded_force is not None:\n for alchemical_force in _add_custom_vdw_lambda(\n custom_nonbonded_force,\n alchemical_indices,\n persistent_indices,\n custom_alchemical_potential,\n ):\n system.addForce(alchemical_force)\n\n elif nonbonded_force is not None:\n for alchemical_force in _add_lj_vdw_lambda(\n nonbonded_force,\n alchemical_indices,\n persistent_indices,\n custom_alchemical_potential,\n ):\n system.addForce(alchemical_force)\n\n else:\n raise NotImplementedError\n\n return system\n
"},{"location":"reference/fep/#absolv.fep.set_fep_lambdas","title":"set_fep_lambdas","text":"set_fep_lambdas(\n context: Context,\n lambda_sterics: float | None = None,\n lambda_electrostatics: float | None = None,\n)\n
Set the values of the alchemical lambdas on an OpenMM context.
Parameters:
-
context (Context) \u2013 The context to update.
-
lambda_sterics (float | None, default: None ) \u2013 The (optional) value of the steric lambda.
-
lambda_electrostatics (float | None, default: None ) \u2013 The (optional) value of the electrostatics lambda.
Source code in absolv/fep.py def set_fep_lambdas(\n context: openmm.Context,\n lambda_sterics: float | None = None,\n lambda_electrostatics: float | None = None,\n):\n \"\"\"Set the values of the alchemical lambdas on an OpenMM context.\n\n Args:\n context: The context to update.\n lambda_sterics: The (optional) value of the steric lambda.\n lambda_electrostatics: The (optional) value of the electrostatics lambda.\n \"\"\"\n\n if lambda_sterics is not None:\n assert (\n 0.0 <= lambda_sterics <= 1.0\n ), f\"`{LAMBDA_STERICS}` must be between 0 and 1\"\n context.setParameter(LAMBDA_STERICS, lambda_sterics)\n\n if lambda_electrostatics is not None:\n assert (\n 0.0 <= lambda_electrostatics <= 1.0\n ), f\"`{LAMBDA_ELECTROSTATICS}` must be between 0 and 1\"\n\n context.setParameter(LAMBDA_ELECTROSTATICS, lambda_electrostatics)\n
"},{"location":"reference/neq/","title":"Run non-equilibrium forward and reverse sampling.
Functions:
-
run_neq \u2013 Run a non-equilibrium simulation with OpenMM, whereby a system is non-reversibly
"},{"location":"reference/neq/#absolv.neq.run_neq","title":"run_neq","text":"run_neq(\n simulation: Simulation,\n coords_0: State,\n coords_1: State,\n protocol: SwitchingProtocol,\n) -> tuple[float, float]\n
Run a non-equilibrium simulation with OpenMM, whereby a system is non-reversibly pulled along an alchemical pathway as described by Ballard and Jarzynski [1] (Figure 3) and Gapsys et al [2].
Both the forward and reverse directions will be simulated.
References [1] Ballard, Andrew J., and Christopher Jarzynski. \"Replica exchange with nonequilibrium switches: Enhancing equilibrium sampling by increasing replica overlap.\" The Journal of chemical physics 136.19 (2012): 194101.
[2] Gapsys, Vytautas, et al. \"Large scale relative protein ligand binding affinities using non-equilibrium alchemy.\" Chemical Science 11.4 (2020): 1140-1152.
Returns:
-
tuple[float, float] \u2013 The forward and reverse work values.
Source code in absolv/neq.py def run_neq(\n simulation: openmm.app.Simulation,\n coords_0: openmm.State,\n coords_1: openmm.State,\n protocol: absolv.config.SwitchingProtocol,\n) -> tuple[float, float]:\n \"\"\"Run a non-equilibrium simulation with OpenMM, whereby a system is non-reversibly\n pulled along an alchemical pathway as described by Ballard and Jarzynski [1]\n (Figure 3) and Gapsys et al [2].\n\n Both the forward and reverse directions will be simulated.\n\n References:\n [1] Ballard, Andrew J., and Christopher Jarzynski. \"Replica exchange with\n nonequilibrium switches: Enhancing equilibrium sampling by increasing replica\n overlap.\" The Journal of chemical physics 136.19 (2012): 194101.\n\n [2] Gapsys, Vytautas, et al. \"Large scale relative protein ligand binding\n affinities using non-equilibrium alchemy.\" Chemical Science 11.4 (2020):\n 1140-1152.\n\n Returns:\n The forward and reverse work values.\n \"\"\"\n\n forward_potentials = _simulate(\n simulation, coords_0, protocol, reverse_direction=False\n )\n reverse_potentials = _simulate(\n simulation, coords_1, protocol, reverse_direction=True\n )\n\n forward_work = (forward_potentials[:, 1] - forward_potentials[:, 0]).sum()\n reverse_work = (reverse_potentials[:, 1] - reverse_potentials[:, 0]).sum()\n\n return forward_work, reverse_work\n
"},{"location":"reference/runner/","title":"Run calculations defined by a config.
Classes:
-
PreparedSystem \u2013 A container for the prepared inputs for a particular solvent phase.
Functions:
-
setup \u2013 Prepare each system to be simulated, namely the ligand in each solvent.
-
run_eq \u2013 Perform a simulation at each lambda window and for each solvent.
-
run_neq \u2013 Performs the simulations required to estimate the free energy using a
"},{"location":"reference/runner/#absolv.runner.PreparedSystem","title":"PreparedSystem","text":" Bases: NamedTuple
A container for the prepared inputs for a particular solvent phase.
Attributes:
-
system (System) \u2013 The alchemically modified OpenMM system.
-
topology (Topology) \u2013 The OpenFF topology with any box vectors set.
-
coords (Quantity) \u2013 The coordinates of the system.
"},{"location":"reference/runner/#absolv.runner.PreparedSystem.system","title":"system instance-attribute","text":"system: System\n
The alchemically modified OpenMM system.
"},{"location":"reference/runner/#absolv.runner.PreparedSystem.topology","title":"topology instance-attribute","text":"topology: Topology\n
The OpenFF topology with any box vectors set.
"},{"location":"reference/runner/#absolv.runner.PreparedSystem.coords","title":"coords instance-attribute","text":"coords: Quantity\n
The coordinates of the system.
"},{"location":"reference/runner/#absolv.runner.setup","title":"setup","text":"setup(\n system: System,\n config: Config,\n force_field: ForceField | SystemGenerator,\n custom_alchemical_potential: str | None = None,\n) -> tuple[PreparedSystem, PreparedSystem]\n
Prepare each system to be simulated, namely the ligand in each solvent.
Parameters:
-
system (System) \u2013 The system to prepare.
-
config (Config) \u2013 The config defining the calculation to perform.
-
force_field (ForceField | SystemGenerator) \u2013 The force field, or a callable that transforms an OpenFF topology into an OpenMM system, without any alchemical modifications to run the calculations using.
If a callable is specified, it should take arguments of an OpenFF topology, a unit wrapped numpy array of atom coords, and a string literal with a value of either \"solvent-a\" or \"solvent-b\".
-
custom_alchemical_potential (str | None, default: None ) \u2013 A custom expression to use for the potential energy function that describes the chemical-alchemical intermolecular interactions.
See the absolv.fep.apply_fep function for more details.
Returns:
-
tuple[PreparedSystem, PreparedSystem] \u2013 The two prepared systems, corresponding to solvent-a and solvent-b respectively.
Source code in absolv/runner.py def setup(\n system: absolv.config.System,\n config: absolv.config.Config,\n force_field: openff.toolkit.ForceField | absolv.utils.openmm.SystemGenerator,\n custom_alchemical_potential: str | None = None,\n) -> tuple[PreparedSystem, PreparedSystem]:\n \"\"\"Prepare each system to be simulated, namely the ligand in each solvent.\n\n Args:\n system: The system to prepare.\n config: The config defining the calculation to perform.\n force_field: The force field, or a callable that transforms an OpenFF\n topology into an OpenMM system, **without** any alchemical modifications\n to run the calculations using.\n\n If a callable is specified, it should take arguments of an OpenFF\n topology, a unit wrapped numpy array of atom coords, and a string\n literal with a value of either ``\"solvent-a\"`` or ``\"solvent-b\"``.\n custom_alchemical_potential: A custom expression to use for the potential\n energy function that describes the chemical-alchemical intermolecular\n interactions.\n\n See the ``absolv.fep.apply_fep`` function for more details.\n\n Returns:\n The two prepared systems, corresponding to solvent-a and solvent-b respectively.\n \"\"\"\n\n solvated_a = _setup_solvent(\n \"solvent-a\",\n system.components_a,\n force_field,\n system.n_solute_molecules,\n system.n_solvent_molecules_a,\n custom_alchemical_potential,\n )\n solvated_b = _setup_solvent(\n \"solvent-b\",\n system.components_b,\n force_field,\n system.n_solute_molecules,\n system.n_solvent_molecules_b,\n custom_alchemical_potential,\n )\n\n if system.solvent_a is not None and config.pressure is not None:\n absolv.utils.openmm.add_barostat(\n solvated_a.system, config.temperature, config.pressure\n )\n if system.solvent_b is not None and config.pressure is not None:\n absolv.utils.openmm.add_barostat(\n solvated_b.system, config.temperature, config.pressure\n )\n\n return solvated_a, solvated_b\n
"},{"location":"reference/runner/#absolv.runner.run_eq","title":"run_eq","text":"run_eq(\n config: Config,\n prepared_system_a: PreparedSystem,\n prepared_system_b: PreparedSystem,\n platform: OpenMMPlatform = femto.md.constants.OpenMMPlatform.CUDA,\n output_dir: Path | None = None,\n) -> Result\n
Perform a simulation at each lambda window and for each solvent.
Parameters:
-
config (Config) \u2013 The config defining the calculation to perform.
-
prepared_system_a (PreparedSystem) \u2013 The prepared system a. See setup for more details.
-
prepared_system_b (PreparedSystem) \u2013 The prepared system b. See setup for more details.
-
platform (OpenMMPlatform, default: CUDA ) \u2013 The OpenMM platform to run using.
-
output_dir (Path | None, default: None ) \u2013 The (optional) directory to save HREMD samples to.
Source code in absolv/runner.py def run_eq(\n config: absolv.config.Config,\n prepared_system_a: PreparedSystem,\n prepared_system_b: PreparedSystem,\n platform: femto.md.constants.OpenMMPlatform = femto.md.constants.OpenMMPlatform.CUDA,\n output_dir: pathlib.Path | None = None,\n) -> absolv.config.Result:\n \"\"\"Perform a simulation at each lambda window and for each solvent.\n\n Args:\n config: The config defining the calculation to perform.\n prepared_system_a: The prepared system a. See ``setup`` for more details.\n prepared_system_b: The prepared system b. See ``setup`` for more details.\n platform: The OpenMM platform to run using.\n output_dir: The (optional) directory to save HREMD samples to.\n \"\"\"\n\n results_a, overlap_a = _run_phase_hremd(\n config.alchemical_protocol_a,\n config.temperature,\n prepared_system_a,\n platform,\n None if output_dir is None else output_dir / \"solvant-a\",\n )\n\n dg_a, dg_a_std = results_a[\"ddG_kcal_mol\"], results_a[\"ddG_error_kcal_mol\"]\n # overlap_a = overlap_a[\"overlap_0\"]\n\n results_b, overlap_b = _run_phase_hremd(\n config.alchemical_protocol_b,\n config.temperature,\n prepared_system_b,\n platform,\n None if output_dir is None else output_dir / \"solvant-b\",\n )\n\n dg_b, dg_b_std = results_b[\"ddG_kcal_mol\"], results_b[\"ddG_error_kcal_mol\"]\n # overlap_b = overlap_b[\"overlap_0\"]\n\n return absolv.config.Result(\n dg_solvent_a=dg_a * openmm.unit.kilocalorie_per_mole,\n dg_std_solvent_a=dg_a_std * openmm.unit.kilocalorie_per_mole,\n dg_solvent_b=dg_b * openmm.unit.kilocalorie_per_mole,\n dg_std_solvent_b=dg_b_std * openmm.unit.kilocalorie_per_mole,\n )\n
"},{"location":"reference/runner/#absolv.runner.run_neq","title":"run_neq","text":"run_neq(\n config: Config,\n prepared_system_a: PreparedSystem,\n prepared_system_b: PreparedSystem,\n platform: OpenMMPlatform = femto.md.constants.OpenMMPlatform.CUDA,\n) -> Result\n
Performs the simulations required to estimate the free energy using a non-equilibrium method.
These include equilibrium simulations at the two end states (i.e. fully interacting and fully de-coupled solute) for each solvent followed by non-equilibrium switching simulations between each end state to compute the forward and reverse work values.
Parameters:
-
config (Config) \u2013 The config defining the calculation to perform.
-
prepared_system_a (PreparedSystem) \u2013 The prepared system a. See setup for more details.
-
prepared_system_b (PreparedSystem) \u2013 The prepared system b. See setup for more details.
-
platform (OpenMMPlatform, default: CUDA ) \u2013 The OpenMM platform to run using.
Source code in absolv/runner.py def run_neq(\n config: absolv.config.Config,\n prepared_system_a: PreparedSystem,\n prepared_system_b: PreparedSystem,\n platform: femto.md.constants.OpenMMPlatform = femto.md.constants.OpenMMPlatform.CUDA,\n) -> absolv.config.Result:\n \"\"\"Performs the simulations required to estimate the free energy using a\n non-equilibrium method.\n\n These include **equilibrium** simulations at the two end states (i.e. fully\n interacting and fully de-coupled solute) for each solvent followed by\n non-equilibrium switching simulations between each end state to compute the\n forward and reverse work values.\n\n Args:\n config: The config defining the calculation to perform.\n prepared_system_a: The prepared system a. See ``setup`` for more details.\n prepared_system_b: The prepared system b. See ``setup`` for more details.\n platform: The OpenMM platform to run using.\n \"\"\"\n\n with tempfile.TemporaryDirectory() as tmp_dir:\n tmp_dir = pathlib.Path(tmp_dir)\n\n solvent_a_dir = tmp_dir / \"solvent-a\"\n solvent_b_dir = tmp_dir / \"solvent-b\"\n\n _run_phase_end_states(\n config.alchemical_protocol_a,\n config.temperature,\n prepared_system_a,\n solvent_a_dir,\n platform,\n )\n dg_a, dg_std_a = _run_switching(\n config.alchemical_protocol_a,\n config.temperature,\n prepared_system_a,\n solvent_a_dir,\n platform,\n )\n\n _run_phase_end_states(\n config.alchemical_protocol_b,\n config.temperature,\n prepared_system_b,\n solvent_b_dir,\n platform,\n )\n dg_b, dg_std_b = _run_switching(\n config.alchemical_protocol_b,\n config.temperature,\n prepared_system_b,\n solvent_b_dir,\n platform,\n )\n\n return absolv.config.Result(\n dg_solvent_a=dg_a.in_units_of(openmm.unit.kilocalories_per_mole),\n dg_std_solvent_a=dg_std_a.in_units_of(openmm.unit.kilocalories_per_mole),\n dg_solvent_b=dg_b.in_units_of(openmm.unit.kilocalories_per_mole),\n dg_std_solvent_b=dg_std_b.in_units_of(openmm.unit.kilocalories_per_mole),\n )\n
"},{"location":"reference/setup/","title":"Setup systems ready for calculations.
Functions:
-
setup_system \u2013 Generate a set of molecule coordinate by using the PACKMOL package.
"},{"location":"reference/setup/#absolv.setup.setup_system","title":"setup_system","text":"setup_system(\n components: list[tuple[str, int]],\n box_target_density: Quantity = 0.95 * _G_PER_ML,\n box_scale_factor: float = 1.1,\n box_padding: Quantity = 2.0 * openmm.unit.angstrom,\n tolerance: Quantity = 2.0 * openmm.unit.angstrom,\n) -> tuple[Topology, Quantity]\n
Generate a set of molecule coordinate by using the PACKMOL package.
Parameters:
-
components (list[tuple[str, int]]) \u2013 A list of the form components[i] = (smiles_i, count_i) where smiles_i is the SMILES representation of component i and count_i is the number of corresponding instances of that component to create.
-
box_target_density (Quantity, default: 0.95 * _G_PER_ML ) \u2013 Target mass density when approximating the box size for the final system with units compatible with g / mL.
-
box_scale_factor (float, default: 1.1 ) \u2013 The amount to scale the approximate box size by.
-
box_padding (Quantity, default: 2.0 * angstrom ) \u2013 The amount of extra padding to add to the box size to avoid PBC issues in units compatible with angstroms.
-
tolerance (Quantity, default: 2.0 * angstrom ) \u2013 The minimum spacing between molecules during packing in units compatible with angstroms.
Returns:
-
tuple[Topology, Quantity] \u2013 A topology containing the molecules the coordinates were generated for and a unit [A] wrapped numpy array of coordinates with shape=(n_atoms, 3).
Source code in absolv/setup.py def setup_system(\n components: list[tuple[str, int]],\n box_target_density: openmm.unit.Quantity = 0.95 * _G_PER_ML,\n box_scale_factor: float = 1.1,\n box_padding: openmm.unit.Quantity = 2.0 * openmm.unit.angstrom,\n tolerance: openmm.unit.Quantity = 2.0 * openmm.unit.angstrom,\n) -> tuple[openff.toolkit.Topology, openmm.unit.Quantity]:\n \"\"\"Generate a set of molecule coordinate by using the PACKMOL package.\n\n Args:\n components: A list of the form ``components[i] = (smiles_i, count_i)`` where\n ``smiles_i`` is the SMILES representation of component `i` and\n ``count_i`` is the number of corresponding instances of that component\n to create.\n box_target_density: Target mass density when approximating the box size for the\n final system with units compatible with g / mL.\n box_scale_factor: The amount to scale the approximate box size by.\n box_padding: The amount of extra padding to add to the box size to avoid PBC\n issues in units compatible with angstroms.\n tolerance: The minimum spacing between molecules during packing in units\n compatible with angstroms.\n\n Raises:\n * PACKMOLRuntimeError\n\n Returns:\n A topology containing the molecules the coordinates were generated for and\n a unit [A] wrapped numpy array of coordinates with shape=(n_atoms, 3).\n \"\"\"\n\n packmol_path = shutil.which(\"packmol\")\n\n if packmol_path is None:\n raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), \"packmol\")\n\n box_size = (\n _approximate_box_size_by_density(components, box_target_density)\n * box_scale_factor\n )\n\n molecules = {}\n\n for smiles, _ in components:\n if smiles in molecules:\n continue\n\n molecule = openff.toolkit.Molecule.from_smiles(smiles)\n molecule.generate_conformers(n_conformers=1)\n molecule.name = f\"component-{len(molecules)}.xyz\"\n molecules[smiles] = molecule\n\n with openff.utilities.temporary_cd():\n for molecule in molecules.values():\n molecule.to_file(molecule.name, \"xyz\")\n\n input_file_contents = _generate_input_file(\n [(molecules[smiles].name, count) for smiles, count in components],\n box_size,\n tolerance,\n )\n\n with open(\"input.txt\", \"w\") as file:\n file.write(input_file_contents)\n\n with open(\"input.txt\") as file:\n subprocess.run(packmol_path, stdin=file, check=True, capture_output=True)\n\n with open(\"output.xyz\") as file:\n output_lines = file.read().splitlines(False)\n\n coordinates = (\n numpy.array(\n [\n [float(coordinate) for coordinate in coordinate_line.split()[1:]]\n for coordinate_line in output_lines[2:]\n if len(coordinate_line) > 0\n ]\n )\n * openmm.unit.angstrom\n )\n\n topology = openff.toolkit.Topology.from_molecules(\n [molecules[smiles] for smiles, count in components for _ in range(count)]\n )\n topology.box_vectors = numpy.eye(3) * (box_size + box_padding)\n\n return topology, coordinates\n
"},{"location":"reference/utils/","title":"Index","text":""},{"location":"reference/utils/#absolv.utils","title":"utils","text":"Common utils
Modules:
-
openmm \u2013 Utilities to manipulate OpenMM objects.
-
topology \u2013 Utilities for manipulating OpenFF topology objects.
"},{"location":"reference/utils/openmm/","title":"Utilities to manipulate OpenMM objects.
Functions:
-
add_barostat \u2013 Add a barostat to a system in-place.
-
create_simulation \u2013 Creates an OpenMM simulation object.
-
create_system_generator \u2013 Creates a 'system generator' that can be used when setting up an alchemical
-
extract_frame \u2013 Extracts a frame from a trajectory as an OpenMM state object.
"},{"location":"reference/utils/openmm/#absolv.utils.openmm.add_barostat","title":"add_barostat","text":"add_barostat(\n system: System,\n temperature: Quantity,\n pressure: Quantity,\n frequency: int = 25,\n)\n
Add a barostat to a system in-place.
Parameters:
-
system (System) \u2013 The system to add the barostat to.
-
temperature (Quantity) \u2013 The temperature to simulate at.
-
pressure (Quantity) \u2013 The pressure to simulate at.
-
frequency (int, default: 25 ) \u2013 The frequency at which to apply the barostat.
Source code in absolv/utils/openmm.py def add_barostat(\n system: openmm.System,\n temperature: openmm.unit.Quantity,\n pressure: openmm.unit.Quantity,\n frequency: int = 25,\n):\n \"\"\"Add a barostat to a system in-place.\n\n Args:\n system: The system to add the barostat to.\n temperature: The temperature to simulate at.\n pressure: The pressure to simulate at.\n frequency: The frequency at which to apply the barostat.\n \"\"\"\n\n barostats = [\n force\n for force in system.getForces()\n if isinstance(force, openmm.MonteCarloBarostat)\n ]\n assert len(barostats) == 0, \"the system should not already contain a barostat\"\n\n system.addForce(openmm.MonteCarloBarostat(pressure, temperature, frequency))\n
"},{"location":"reference/utils/openmm/#absolv.utils.openmm.create_simulation","title":"create_simulation","text":"create_simulation(\n system: System,\n topology: Topology,\n coords: Quantity,\n integrator: Integrator,\n platform: OpenMMPlatform,\n) -> Simulation\n
Creates an OpenMM simulation object.
Parameters:
-
system (System) \u2013 The system to simulate
-
topology (Topology) \u2013 The topology being simulated.
-
coords (Quantity) \u2013 The initial coordinates. Box vectors (if any) will be taken from the topology.
-
integrator (Integrator) \u2013 The integrator to evolve the system with.
-
platform (OpenMMPlatform) \u2013 The accelerator to run using.
Returns:
-
Simulation \u2013 The created simulation.
Source code in absolv/utils/openmm.py def create_simulation(\n system: openmm.System,\n topology: openff.toolkit.Topology,\n coords: openmm.unit.Quantity,\n integrator: openmm.Integrator,\n platform: femto.md.constants.OpenMMPlatform,\n) -> openmm.app.Simulation:\n \"\"\"Creates an OpenMM simulation object.\n\n Args:\n system: The system to simulate\n topology: The topology being simulated.\n coords: The initial coordinates. Box vectors (if any) will be taken from the\n topology.\n integrator: The integrator to evolve the system with.\n platform: The accelerator to run using.\n\n Returns:\n The created simulation.\n \"\"\"\n platform_properties = (\n {\"Precision\": \"mixed\"} if platform.upper() in [\"CUDA\", \"OPENCL\"] else {}\n )\n platform = openmm.Platform.getPlatformByName(platform)\n\n if topology.box_vectors is not None:\n system.setDefaultPeriodicBoxVectors(*topology.box_vectors.to_openmm())\n\n simulation = openmm.app.Simulation(\n topology.to_openmm(), system, integrator, platform, platform_properties\n )\n\n if topology.box_vectors is not None:\n simulation.context.setPeriodicBoxVectors(*topology.box_vectors.to_openmm())\n\n simulation.context.setPositions(coords)\n simulation.context.setVelocitiesToTemperature(integrator.getTemperature())\n\n return simulation\n
"},{"location":"reference/utils/openmm/#absolv.utils.openmm.create_system_generator","title":"create_system_generator","text":"create_system_generator(\n force_field: ForceField,\n solvent_a_nonbonded_method: int,\n solvent_b_nonbonded_method: int,\n nonbonded_cutoff: Quantity = 1.0\n * openmm.unit.nanometer,\n constraints: int | None = None,\n rigid_water: bool | None = None,\n remove_cmm_motion: bool = True,\n hydrogen_mass: Quantity | None = None,\n switch_distance: Quantity | None = None,\n) -> SystemGenerator\n
Creates a 'system generator' that can be used when setting up an alchemical free energy calculation from an OpenMM force field.
Parameters:
-
force_field (ForceField) \u2013 The OpenMM force field to parameterize the topology using.
-
solvent_a_nonbonded_method (int) \u2013 The non-bonded method to use in solvent a.
-
solvent_b_nonbonded_method (int) \u2013 The non-bonded method to use in solvent b.
-
nonbonded_cutoff (Quantity, default: 1.0 * nanometer ) \u2013 The non-bonded cutoff to use.
-
constraints (int | None, default: None ) \u2013 The type of constraints to apply to the system.
-
rigid_water (bool | None, default: None ) \u2013 Whether to force rigid water.
-
remove_cmm_motion (bool, default: True ) \u2013 Whether to remove any CMM motion.
-
hydrogen_mass (Quantity | None, default: None ) \u2013 The mass to use for hydrogens.
-
switch_distance (Quantity | None, default: None ) \u2013 The switch distance to use.
Returns:
-
SystemGenerator \u2013 A callable that will create an OpenMM system from an OpenFF topology and the name of the solvent (i.e. \"solvent-a\" or \"solvent-b\") the system will be used for.
Source code in absolv/utils/openmm.py def create_system_generator(\n force_field: openmm.app.ForceField,\n solvent_a_nonbonded_method: int,\n solvent_b_nonbonded_method: int,\n nonbonded_cutoff: openmm.unit.Quantity = 1.0 * openmm.unit.nanometer,\n constraints: int | None = None,\n rigid_water: bool | None = None,\n remove_cmm_motion: bool = True,\n hydrogen_mass: openmm.unit.Quantity | None = None,\n switch_distance: openmm.unit.Quantity | None = None,\n) -> SystemGenerator:\n \"\"\"Creates a 'system generator' that can be used when setting up an alchemical\n free energy calculation from an OpenMM force field.\n\n Args:\n force_field: The OpenMM force field to parameterize the topology using.\n solvent_a_nonbonded_method: The non-bonded method to use in solvent a.\n solvent_b_nonbonded_method: The non-bonded method to use in solvent b.\n nonbonded_cutoff: The non-bonded cutoff to use.\n constraints: The type of constraints to apply to the system.\n rigid_water: Whether to force rigid water.\n remove_cmm_motion: Whether to remove any CMM motion.\n hydrogen_mass: The mass to use for hydrogens.\n switch_distance: The switch distance to use.\n\n Returns:\n A callable that will create an OpenMM system from an OpenFF topology and the\n name of the solvent (i.e. ``\"solvent-a\"`` or ``\"solvent-b\"``) the system will\n be used for.\n \"\"\"\n\n def system_generator(\n topology: openff.toolkit.Topology,\n coordinates: openmm.unit.Quantity,\n solvent_idx: typing.Literal[\"solvent-a\", \"solvent-b\"],\n ) -> openmm.System:\n openmm_topology = topology.to_openmm()\n\n if topology.box_vectors is not None:\n openmm_topology.setPeriodicBoxVectors(topology.box_vectors.to_openmm())\n\n # We need to fix the special case of water in order for OMM to correctly apply\n # a constraint between H atoms.\n for chain in openmm_topology.chains():\n for residue in chain.residues():\n if len(residue) != 3:\n continue\n\n symbols = sorted(atom.element.symbol for atom in residue.atoms())\n\n if symbols == [\"H\", \"H\", \"O\"]:\n residue.name = \"HOH\"\n\n from openmm.app import Modeller\n\n modeller = Modeller(\n openmm_topology,\n [\n openmm.Vec3(coordinate[0], coordinate[1], coordinate[2])\n for coordinate in coordinates.value_in_unit(openmm.unit.nanometers)\n ]\n * openmm.unit.nanometers,\n )\n modeller.addExtraParticles(force_field)\n\n system = force_field.createSystem(\n modeller.getTopology(),\n nonbondedMethod=(\n solvent_a_nonbonded_method\n if solvent_idx == \"solvent-a\"\n else solvent_b_nonbonded_method\n ),\n nonbondedCutoff=nonbonded_cutoff,\n constraints=constraints,\n rigidWater=rigid_water,\n removeCMMotion=remove_cmm_motion,\n hydrogenMass=hydrogen_mass,\n switchDistance=switch_distance,\n )\n\n return system\n\n return system_generator\n
"},{"location":"reference/utils/openmm/#absolv.utils.openmm.extract_frame","title":"extract_frame","text":"extract_frame(trajectory: Trajectory, idx: int) -> State\n
Extracts a frame from a trajectory as an OpenMM state object.
Parameters:
-
trajectory (Trajectory) \u2013 The trajectory to extract the frame from.
-
idx (int) \u2013 The index of the frame to extract.
Returns:
-
State \u2013 The extracted frame.
Source code in absolv/utils/openmm.py def extract_frame(trajectory: mdtraj.Trajectory, idx: int) -> openmm.State:\n \"\"\"Extracts a frame from a trajectory as an OpenMM state object.\n\n Args:\n trajectory: The trajectory to extract the frame from.\n idx: The index of the frame to extract.\n\n Returns:\n The extracted frame.\n \"\"\"\n\n system = openmm.System()\n\n for _ in range(trajectory.n_atoms):\n system.addParticle(1.0 * openmm.unit.dalton)\n\n context = openmm.Context(system, openmm.VerletIntegrator(0.0001))\n context.setPositions(trajectory.openmm_positions(idx))\n\n if trajectory.unitcell_vectors is not None:\n context.setPeriodicBoxVectors(*trajectory.openmm_boxes(idx))\n\n return context.getState(getPositions=True)\n
"},{"location":"reference/utils/topology/","title":"Utilities for manipulating OpenFF topology objects.
Functions:
-
topology_to_components \u2013 A helper method for condensing a topology down to a list of components
-
topology_to_atom_indices \u2013 A helper method for extracting the sets of atom indices associated with each
"},{"location":"reference/utils/topology/#absolv.utils.topology.topology_to_components","title":"topology_to_components","text":"topology_to_components(\n topology: Topology,\n) -> list[tuple[str, int]]\n
A helper method for condensing a topology down to a list of components and their counts.
Notes If the topology is not contiguous then the returned list may contain multiple tuples with the same smiles but different counts.
Parameters:
-
topology (Topology) \u2013 The topology to condense.
Returns:
-
list[tuple[str, int]] \u2013 A list of the form components[i] = (smiles_i, count_i) where smiles_i is the SMILES representation of component i and count_i is the number of corresponding instances of that component to create.
Source code in absolv/utils/topology.py def topology_to_components(topology: openff.toolkit.Topology) -> list[tuple[str, int]]:\n \"\"\"A helper method for condensing a topology down to a list of components\n and their counts.\n\n Notes:\n If the topology is not contiguous then the returned list may contain multiple\n tuples with the same smiles but different counts.\n\n Args:\n topology: The topology to condense.\n\n Returns:\n A list of the form ``components[i] = (smiles_i, count_i)`` where\n ``smiles_i`` is the SMILES representation of component `i` and\n ``count_i`` is the number of corresponding instances of that component\n to create.\n \"\"\"\n components = []\n\n current_smiles = None\n current_count = 0\n\n for molecule in topology.molecules:\n smiles = molecule.to_smiles()\n\n if smiles == current_smiles:\n current_count += 1\n continue\n\n if current_count > 0:\n components.append((current_smiles, current_count))\n\n current_smiles = smiles\n current_count = 1\n\n if current_count > 0:\n components.append((current_smiles, current_count))\n\n return components\n
"},{"location":"reference/utils/topology/#absolv.utils.topology.topology_to_atom_indices","title":"topology_to_atom_indices","text":"topology_to_atom_indices(\n topology: Topology,\n) -> list[set[int]]\n
A helper method for extracting the sets of atom indices associated with each molecule in a topology.
Parameters:
-
topology (Topology) \u2013 The topology to extract the atom indices from.
Returns:
-
list[set[int]] \u2013 The set of atoms indices associated with each molecule in the topology.
Source code in absolv/utils/topology.py def topology_to_atom_indices(topology: openff.toolkit.Topology) -> list[set[int]]:\n \"\"\"A helper method for extracting the sets of atom indices associated with each\n molecule in a topology.\n\n Args:\n topology: The topology to extract the atom indices from.\n\n Returns:\n The set of atoms indices associated with each molecule in the topology.\n \"\"\"\n\n atom_indices: list[set[int]] = []\n current_atom_idx = 0\n\n for molecule in topology.molecules:\n atom_indices.append({i + current_atom_idx for i in range(molecule.n_atoms)})\n current_atom_idx += molecule.n_atoms\n\n return atom_indices\n
"},{"location":"user-guide/overview/","title":"Overview","text":"To begin with, we define a configuration that encodes the entirety of the absolute free energy calculation.
This includes specifying the solutes and the two solvents that they will be transferred between:
import absolv.config\n\nsystem=absolv.config.System(\n solutes={\"CCO\": 1}, solvent_a=None, solvent_b={\"O\": 895}\n)\n
Each key should be a SMILES representation of a molecule, and the value should be the number of copies of that molecule to include in the system. There may be multiple in the case of, e.g., ion pairs like Na+ and Cl-. None may be used to specify vacuum, e.g., in the above case that the solute will be transferred from vacuum into bulk water.
The temperature and pressure that the calculation will be performed at must also be specified:
import openmm.unit\n\ntemperature=298.15 * openmm.unit.kelvin,\npressure=1.0 * openmm.unit.atmosphere\n
Finally, the alchemical pathway to transform the solute along in each solvent must be specified. This can either be a more traditional 'equilibrium' pathway, or a 'non-equilibrium' pathway:
EquilibriumNon-equilibrium import absolv.config\n\nalchemical_protocol_a=absolv.config.EquilibriumProtocol(\n lambda_sterics=[1.0, 1.0, 1.0, 1.0, 1.0],\n lambda_electrostatics=[1.0, 0.75, 0.5, 0.25, 0.0]\n)\nalchemical_protocol_b=absolv.config.EquilibriumProtocol(\n lambda_sterics=[\n 1.00, 1.00, 1.00, 1.00, 1.00, 0.95, 0.90, 0.80, 0.70, 0.60, 0.50, 0.40,\n 0.35, 0.30, 0.25, 0.20, 0.15, 0.10, 0.05, 0.00,\n ],\n lambda_electrostatics=[\n 1.00, 0.75, 0.50, 0.25, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00,\n 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00,\n ]\n)\n
Here the default lambda schedule from FreeSolv is used. A lambda of 1.0 indicates a fully interacting system, and a lambda of 0.0 indicates a fully decoupled system.
import femto.md.config\nimport openmm.unit\n\nimport absolv.config\n\nintegrator = femto.md.config.LangevinIntegrator(\n timestep=2.0 * openmm.unit.femtoseconds\n)\n\nalchemical_protocol_a=absolv.config.NonEquilibriumProtocol(\n # the protocol to use for the production run at each end state\n production_protocol=absolv.config.SimulationProtocol(\n integrator=integrator, n_steps=6250 * 160\n ),\n # the interval with which to store frames that NEQ switching\n # simulations will be launched from\n production_report_interval=6250,\n # define how the NEQ switching will be performed\n switching_protocol=absolv.config.SwitchingProtocol(\n n_electrostatic_steps=60,\n n_steps_per_electrostatic_step=100,\n # intra-molecular vdW interactions are not decoupled by default, so we\n # don't need to do any vdW decoupling in vacuum when there's only one solvent\n n_steric_steps=0,\n n_steps_per_steric_step=0\n )\n)\nalchemical_protocol_b=absolv.config.NonEquilibriumProtocol(\n production_protocol=absolv.config.SimulationProtocol(\n integrator=integrator, n_steps=6250 * 160,\n ),\n switching_protocol=absolv.config.SwitchingProtocol(\n # Annihilate the electrostatic interactions over the first 12 ps\n n_electrostatic_steps=60,\n n_steps_per_electrostatic_step=100,\n # followed by decoupling the vdW interactions over the next 38 ps\n n_steric_steps=190,\n n_steps_per_steric_step=100,\n )\n)\n
These individual components are then combined into a single configuration object:
import absolv.config\n\nconfig = absolv.config.Config(\n temperature=temperature,\n pressure=pressure,\n alchemical_protocol_a=alchemical_protocol_a,\n alchemical_protocol_b=alchemical_protocol_b,\n)\n
which can be used to trivially set up
import openff.toolkit\nforce_field = openff.toolkit.ForceField(\"openff-2.1.0.offxml\")\n\nimport absolv.runner\nprepared_system_a, prepared_system_b = absolv.runner.setup(system, config, force_field)\n
and run the calculation:
EquilibriumNon-equilibrium result = absolv.runner.run_eq(\n config, prepared_system_a, prepared_system_b, \"CUDA\"\n)\nprint(result)\n
result = absolv.runner.run_neq(\n config, prepared_system_a, prepared_system_b, \"CUDA\"\n)\nprint(result)\n
where the result will be a Result object.
"},{"location":"user-guide/theory/","title":"Theory","text":"The change in free energy of transferring a solute from one solvent to another \\(\\Delta G_{A->B}\\), can be readily computed by constructing a thermodynamic cycle in which the solute in each solvent is alchemically transformed into a non-interacting state as shown below:
Namely,
\\[ \\Delta G_{A->B} = \\Delta G_1 - \\Delta G_3 \\] where here \\(\\Delta G_1\\) and \\(\\Delta G_3\\) are the change in free energy of alchemically transforming a solute (or perhaps multiple solutes in the case of charged molecules with counter ions) so that it no longer interacts with the surrounding corresponding 'solvent'.
Computing the solvation (/ hydration) free energy is a special case of this cycle when the first 'solvent' is vacuum.
This framework currently offers two routes to computing such a change in free energy as the solute is being alchemically transformed, a more commonly used 'equilibrium route' and a 'non-equilibrium' route.
Note For more information about how the different inter- and intramolecular interactions are alchemically modified see the transformations page.
"},{"location":"user-guide/theory/#equilibrium-calculations","title":"Equilibrium Calculations","text":"Within this framework we refer to free energy calculations that involve dividing the alchemical pathway into discrete windows at each value of the coupling parameter \\(\\lambda\\), collecting equilibrium samples for each such discrete state, and using these samples to compute the free energy using an approach such as thermodynamic integration (TI) 1, BAR 2, and MBAR 3 as 'equilibrium' free energy calculations.
At present absolv does not offer functionality for computing the derivatives with respect to lambda required for TI, and only supports MBAR and technically BAR although this estimator is not recommended.
See the overview for more information on running equilibrium free energy calculations using absolv.
"},{"location":"user-guide/theory/#non-equilibrium-calculations","title":"Non-equilibrium Calculations","text":"Within this framework we refer to free energy calculations that proceed by generating equilibrium configurations at both end states (with the solute fully interacting with the solute and with the solute-solvent interactions fully decoupled / annihilated), and then driving each configuration non-reverisbly along the alchemical pathway by scaling the coupling factor \\(\\lambda\\) as a function of time 4.
From a practical perspective it is computationally more efficient and convenient to proceed along the alchemical pathway in a stepwise, rather than continuous fashion. More specifically, the protocol proceeds by making a perturbation to \\(\\lambda\\), followed by several relaxation steps, and repeating these two sub-stebs until the alchemical transformation is complete.
in this way the required to transform the solute from the interacting to the non-interacting state and likewise from the non-interacting to the interacting state can be computed according to
\\[W = \\sum_{i=0}^{n-1} \\left[u_{\\lambda_{i+1}}\\left(x_i\\right) - u_{\\lambda_{i}}\\left(x_i\\right)\\right]\\] where here \\(u_{\\lambda_i}\\left(x_i\\right)\\) is the reduced potential evaluated at configuration \\(i\\) and \\(\\lambda_i\\).
The free energy is then estimated by solving
\\[\\sum^N_{i=1}\\dfrac{1}{1+\\exp(\\beta(W^f_i-\\Delta F))} = \\sum^N_{i=1}\\dfrac{1}{1+\\exp(\\beta(W^r_i+\\Delta F))}\\] self consistently where \\(W^f_i\\) corresponds to work computed along the forward pathway going from the interacting to the non-interacting state and \\(W^r_i\\) to work computed along the reverse pathway going from the non-interacting to the interacting state. \\(N\\) refers to the total number of equilibrium snapshots that were generated at each end state.
-
TP Straatsma and JA McCammon. Computational alchemy. Annual Review of Physical Chemistry, 43(1):407\u2013435, 1992.\u00a0\u21a9
-
Charles H Bennett. Efficient estimation of free energy differences from monte carlo data. Journal of Computational Physics, 22(2):245\u2013268, 1976.\u00a0\u21a9
-
Michael R Shirts and John D Chodera. Statistically optimal analysis of samples from multiple equilibrium states. The Journal of chemical physics, 129(12):124105, 2008.\u00a0\u21a9
-
Andrew J Ballard and Christopher Jarzynski. Replica exchange with nonequilibrium switches: enhancing equilibrium sampling by increasing replica overlap. The Journal of chemical physics, 136(19):194101, 2012.\u00a0\u21a9
"},{"location":"user-guide/transformations/","title":"Transformations","text":"The absolv framework supports alchemically transforming the electrostatic and vdW interactions within an OpenMM System object via the absolv.fep.apply_fep function.
Due to the huge flexibility of the OpenMM system, allowing completely custom electrostatic and vdW forces to be implemented, it would be impossible to support all possible combinations of non-bonded forces. As such, the framework currently only systems that contain
- vdW + electrostatic interactions in a single built-in
NonbondedForce object - electrostatic interactions in a built-in
NonbondedForce object, vdW interactions in a CustomNonbondedForce object and vdW 1-4 interactions in a CustomBondForce - the above combinations sans any electrostatics interactions
are supported.
"},{"location":"user-guide/transformations/#electrostatics","title":"Electrostatics","text":"The electrostatic interactions in the system are alchemically transformed by linearly scaling all partial charges on particles in solute molecules by \\(\\lambda_{elec}\\), corresponding to a \"lambda_electrostatics\" context variable that will be added by the factory, such that
\\[q^{sol}_i \\left(\\lambda_{elec}\\right) = \\lambda_{elec} \\times q^{sol}_i\\] All intramolecular interactions will be switched off during the alchemical transformation. This is referred to as 'annihilating' the electrostatic interactions in other packages and some literature.
Because the charges are scaled directly, the energy contributions of the alchemically scaled electrostatic interactions will be
\\[U^E = \\lambda_{elec} U^E_{sol-solv} + \\lambda_{elec}^2 U^E_{sol-sol} + U^E_{solv-solv}\\] where \\(U^E_{sol-sol}\\), \\(U^E_{sol-solv}\\) and \\(U^E_{solv-solv}\\) are the un-scaled electrostatic contributions to the energies of the solute-solute, solute-solvent and solvent-solvent interactions respectively.
"},{"location":"user-guide/transformations/#vdw","title":"vdW","text":"Currently vdW interactions can only be transformed if they are stored in a standard NonbondedForce OR if they are split between a CustomBondForce (1-2, 1-3, and 1-4 interactions) and a CustomNonbondedForce (the remaining 1-n and intermolecular interactions).
The interactions will be transformed according to \\(\\lambda_{vdW}\\) which corresponds to a \"lambda_sterics\" context variable that will be added by the factory.
Only intermolecular vdW interactions will be alchemically scaled, while all intramolecular interactions will be left un-scaled. This is is referred to as 'decoupling' the vdW interactions in other packages and some literature.
"},{"location":"user-guide/transformations/#lennard-jones","title":"Lennard--Jones","text":"If the vdW interactions are stored in a standard NonbondedForce, then the alchemical factory will split them so that
-
the NonbondedForce force retains all interactions between solvent particles
-
all intermolecular alchemical (both solute-solvent and solute-solute) interactions are moved into a new CustomNonbondedForce
-
all solute intramolecular interactions are moved into a new CustomNonbondedForce
Note The intramolecular solute-solute interactions won't use any periodic boundary corrections such that the the decoupled state of the solute corresponds to the proper vacuum state without periodicity effects.
The CustomNonbondedForce will copy over all settings (including cut-off, long-range correction, etc) from the original NonbondedForce, but will replace the normal LJ energy expression with a soft-core version. By default, this takes the form:
\\[U^{vdW} \\left( \\lambda_{vdW} \\right) = \\lambda_{vdW} \\times 4 \\varepsilon \\left[ \\left( \\dfrac{\\sigma}{r_{eff}}\\right)^{12} - \\left( \\dfrac{\\sigma}{r_{eff}}\\right)^{6} \\right]\\] where
\\[r_{eff} = \\sigma \\left( \\dfrac{1}{2} \\left(1 - \\lambda_{vdW}\\right) + \\left( \\dfrac{r}{\\sigma} \\right) ^ 6 \\right) ^ \\frac{1}{6}\\]"},{"location":"user-guide/transformations/#custom-vdw-forms","title":"Custom vdW Forms","text":"If the vdW interactions are split across a CustomNonbondedForce and a CustomBondForce then the alchemical factory will further split them so that
-
the original CustomNonbondedForce force will retain all interactions between solvent particles
-
all solute intramolecular interactions are moved into a new CustomNonbondedForce
Note The intramolecular solute-solute interactions won't use any periodic boundary corrections such that the the decoupled state of the solute corresponds to the proper vacuum state without periodicity effects.
-
all intermolecular alchemical (both solute-solvent and solute-solute) interactions are moved into another new CustomNonbondedForce with a modified energy expression such that
\\[U^{vdW}_{sol-solv} \\left( \\lambda_{vdW} \\right) = \\lambda_{vdW} \\times U^{vdW}_{sol-solv}\\]
"}]}
\ No newline at end of file
diff --git a/latest/sitemap.xml.gz b/latest/sitemap.xml.gz
index 2015dd8..bc9ced2 100644
Binary files a/latest/sitemap.xml.gz and b/latest/sitemap.xml.gz differ