Skip to content
New issue

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

Sign up for GitHub

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

Already on GitHub? Sign in to your account

Re-write avg. ensemble op to use analytical gradient #78

Merged
merged 4 commits into from
Oct 19, 2023
Merged

Re-write avg. ensemble op to use analytical gradient #78

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
ffcf3e6
Re-write avg. ensemble op to use analytical gradient
SimonBoothroyd Oct 19, 2023
642fb16
Add tests and update notebook
SimonBoothroyd Oct 19, 2023
7989ef6
Fix broadcasting
SimonBoothroyd Oct 19, 2023
fdd2a35
Pin interchange
SimonBoothroyd Oct 19, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion devtools/envs/base.yaml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@ dependencies:
# Core packages
- openff-units
- openff-toolkit-base >=0.9.2
- openff-interchange-base >=0.3.15
- openff-interchange-base ==0.3.15

- pytorch-cpu
- pydantic
Expand All @@ -24,6 +24,7 @@ dependencies:
- rdkit
- packmol
- numpy
- msgpack-python

# Examples
- jupyter
Expand Down
176 changes: 107 additions & 69 deletions examples/md-simulations.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@
"application/vnd.jupyter.widget-view+json": {
"version_major": 2,
"version_minor": 0,
"model_id": "cdbe68b205874fd99fc47a70449ef621"
"model_id": "517ec81345264ae8b2a042fd56b568a6"
}
},
"metadata": {},
Expand All @@ -42,16 +42,16 @@
" openff.toolkit.ForceField(\"openff-2.0.0.offxml\"),\n",
" openff.toolkit.Molecule.from_smiles(smiles).to_topology(),\n",
" )\n",
" for smiles in (\"O\", \"CO\")\n",
" for smiles in (\"CCO\", \"CO\")\n",
"]\n",
"\n",
"tensor_ff, topologies = smee.converters.convert_interchange(interchanges)"
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
"end_time": "2023-10-11T23:33:50.423674Z",
"start_time": "2023-10-11T23:33:45.643146Z"
"end_time": "2023-10-19T21:50:43.299571Z",
"start_time": "2023-10-19T21:50:38.880450Z"
}
},
"id": "e178932166969df5"
Expand All @@ -77,8 +77,8 @@
"metadata": {
"collapsed": false,
"ExecuteTime": {
"end_time": "2023-10-11T23:33:50.440248Z",
"start_time": "2023-10-11T23:33:50.423221Z"
"end_time": "2023-10-19T21:50:43.308877Z",
"start_time": "2023-10-19T21:50:43.302488Z"
}
},
"id": "38700a4d251b1ab9"
Expand All @@ -100,18 +100,18 @@
"source": [
"import smee\n",
"\n",
"# define a periodic box containing 216 water molecules\n",
"system_water = smee.TensorSystem([topologies[0]], [216], is_periodic=True)\n",
"# define a periodic box containing 216 ethanol molecules\n",
"system_ethanol = smee.TensorSystem([topologies[0]], [216], is_periodic=True)\n",
"# define a periodic box containing 216 methanol molecules\n",
"system_methanol = smee.TensorSystem([topologies[1]], [216], True)\n",
"# define a periodic box containing 128 water molecules and 128 methanol molecules\n",
"# define a periodic box containing 128 ethanol molecules and 128 methanol molecules\n",
"system_mixture = smee.TensorSystem(topologies, [128, 128], True)"
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
"end_time": "2023-10-11T23:33:50.440423Z",
"start_time": "2023-10-11T23:33:50.429216Z"
"end_time": "2023-10-19T21:50:43.343163Z",
"start_time": "2023-10-19T21:50:43.307306Z"
}
},
"id": "58a584bf7997e194"
Expand All @@ -133,53 +133,86 @@
"execution_count": 4,
"outputs": [],
"source": [
"import tempfile\n",
"\n",
"import openmm.unit\n",
"\n",
"import smee.mm\n",
"\n",
"temperature = 298.15 * openmm.unit.kelvin\n",
"pressure = 1.0 * openmm.unit.atmosphere\n",
"\n",
"def simulation_protocol():\n",
" temperature = 298.15 * openmm.unit.kelvin\n",
" pressure = 1.0 * openmm.unit.atmosphere\n",
"\n",
" coords_config = smee.mm.GenerateCoordsConfig()\n",
" # we can run an arbitrary number of equilibration simulations / minimizations.\n",
" # all generated data will be discarded, but the final coordinates will be used\n",
" # to initialize the production simulation\n",
" equilibrate_config = [\n",
" smee.mm.MinimizationConfig(),\n",
" # short NVT equilibration simulation\n",
" smee.mm.SimulationConfig(\n",
" temperature=temperature,\n",
" pressure=None,\n",
" n_steps=50000,\n",
" timestep=1.0 * openmm.unit.femtosecond,\n",
" ),\n",
" # short NPT equilibration simulation\n",
" smee.mm.SimulationConfig(\n",
" temperature=temperature,\n",
" pressure=pressure,\n",
" n_steps=50000,\n",
" timestep=1.0 * openmm.unit.femtosecond,\n",
" ),\n",
" ]\n",
" # long NPT production simulation\n",
" production_config = smee.mm.SimulationConfig(\n",
"# we can run an arbitrary number of equilibration simulations / minimizations.\n",
"# all generated data will be discarded, but the final coordinates will be used\n",
"# to initialize the production simulation\n",
"equilibrate_config = [\n",
" smee.mm.MinimizationConfig(),\n",
" # short NVT equilibration simulation\n",
" smee.mm.SimulationConfig(\n",
" temperature=temperature,\n",
" pressure=None,\n",
" n_steps=50000,\n",
" timestep=1.0 * openmm.unit.femtosecond,\n",
" ),\n",
" # short NPT equilibration simulation\n",
" smee.mm.SimulationConfig(\n",
" temperature=temperature,\n",
" pressure=pressure,\n",
" n_steps=500000,\n",
" timestep=2.0 * openmm.unit.femtosecond,\n",
" )\n",
" # store coords and values every 2 ps\n",
" report_interval = 1000\n",
" n_steps=50000,\n",
" timestep=1.0 * openmm.unit.femtosecond,\n",
" ),\n",
"]\n",
"# long NPT production simulation\n",
"production_config = smee.mm.SimulationConfig(\n",
" temperature=temperature,\n",
" pressure=pressure,\n",
" n_steps=500000,\n",
" timestep=2.0 * openmm.unit.femtosecond,\n",
")\n",
"\n",
"\n",
"def compute_ensemble_averages(system: smee.TensorSystem):\n",
" # computing the ensemble averages is a two step process - we first need to run\n",
" # an MD simulation using the force field making sure to store the coordinates,\n",
" # box vectors and kinetic energies\n",
" coords, box_vectors = smee.mm.generate_system_coords(system)\n",
"\n",
" with tempfile.NamedTemporaryFile() as tmp_path:\n",
" # save the simulation output every 1000th frame (2 ps) to a temporary file.\n",
" # we could also save the trajectory more permanently, but as we do nothing\n",
" # with it after computing the averages in this example, we simply want to\n",
" # discard it.\n",
"\n",
" reporter_interval = 1\n",
"\n",
" with open(tmp_path.name, \"wb\") as tmp_file:\n",
" reporter = smee.mm.TensorReporter(tmp_file, reporter_interval)\n",
"\n",
" smee.mm.simulate(\n",
" system,\n",
" tensor_ff,\n",
" coords,\n",
" box_vectors,\n",
" equilibrate_config,\n",
" production_config,\n",
" [reporter],\n",
" )\n",
"\n",
" return coords_config, equilibrate_config, production_config, report_interval"
" # we can then compute the ensemble averages from the trajectory. generating\n",
" # the trajectory separately from computing the ensemble averages allows us\n",
" # to run the simulation in parallel with other simulations more easily, without\n",
" # having to worry about copying gradients between workers / processes.\n",
" import pathlib\n",
"\n",
" return smee.mm.compute_ensemble_averages(\n",
" system, tensor_ff, pathlib.Path(tmp_path.name), temperature, pressure\n",
" )"
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
"end_time": "2023-10-11T23:33:50.476680Z",
"start_time": "2023-10-11T23:33:50.435362Z"
"end_time": "2023-10-19T21:50:43.376891Z",
"start_time": "2023-10-19T21:50:43.315296Z"
}
},
"id": "ccba93245cf83ff7"
Expand All @@ -197,25 +230,28 @@
{
"cell_type": "code",
"execution_count": 5,
"outputs": [],
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/simon/mambaforge/envs/smee/lib/python3.11/site-packages/torch/autograd/function.py:506: UserWarning: at::frobenius_norm is deprecated and it is just left for JIT compatibility. It will be removed in a future PyTorch release. Please use `linalg.vector_norm(A, 2., dim, keepdim)` instead (Triggered internally at /Users/runner/miniforge3/conda-bld/pytorch-recipe_1696859382176/work/aten/src/ATen/native/LinearAlgebra.cpp:2789.)\n",
" return super().apply(*args, **kwargs) # type: ignore[misc]\n"
]
}
],
"source": [
"# run simulations of each system and compute ensemble averages over the trajectories\n",
"# of the potential energy, volume, and density\n",
"water_avgs = smee.mm.compute_ensemble_averages(\n",
" system_water, tensor_ff, *simulation_protocol()\n",
")\n",
"methanol_avgs = smee.mm.compute_ensemble_averages(\n",
" system_methanol, tensor_ff, *simulation_protocol()\n",
")\n",
"mixture_avgs = smee.mm.compute_ensemble_averages(\n",
" system_mixture, tensor_ff, *simulation_protocol()\n",
")"
"ethanol_avgs = compute_ensemble_averages(system_ethanol)\n",
"methanol_avgs = compute_ensemble_averages(system_methanol)\n",
"mixture_avgs = compute_ensemble_averages(system_mixture)"
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
"end_time": "2023-10-11T23:33:54.809457Z",
"start_time": "2023-10-11T23:33:50.474518Z"
"end_time": "2023-10-19T21:50:53.873961Z",
"start_time": "2023-10-19T21:50:43.370825Z"
}
},
"id": "3156bcfc509380f7"
Expand All @@ -238,27 +274,29 @@
"outputs": [],
"source": [
"# define some MOCK data and loss function\n",
"mock_water_density = 1.0 # g/mL\n",
"mock_ethanol_density = 0.789 # g/mL\n",
"mock_methanol_density = 0.791 # g/mL\n",
"\n",
"mock_enthalpy_of_mixing = 0.891 # kcal/mol\n",
"\n",
"loss = (water_avgs[\"density\"] - mock_water_density) ** 2\n",
"loss = (ethanol_avgs[\"density\"] - mock_ethanol_density) ** 2\n",
"loss += (methanol_avgs[\"density\"] - mock_methanol_density) ** 2\n",
"\n",
"mixture_enthalpy = mixture_avgs[\"enthalpy\"] / 256\n",
"\n",
"water_enthalpy = water_avgs[\"enthalpy\"] / 128\n",
"ethanol_enthalpy = ethanol_avgs[\"enthalpy\"] / 128\n",
"methanol_enthalpy = methanol_avgs[\"enthalpy\"] / 128\n",
"\n",
"enthalpy_of_mixing = mixture_enthalpy - (0.5 * water_enthalpy + 0.5 * methanol_enthalpy)\n",
"enthalpy_of_mixing = mixture_enthalpy - (\n",
" 0.5 * ethanol_enthalpy + 0.5 * methanol_enthalpy\n",
")\n",
"loss += (enthalpy_of_mixing - mock_enthalpy_of_mixing) ** 2"
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
"end_time": "2023-10-11T23:33:54.815759Z",
"start_time": "2023-10-11T23:33:54.810831Z"
"end_time": "2023-10-19T21:50:53.879337Z",
"start_time": "2023-10-19T21:50:53.876688Z"
}
},
"id": "38b9a27d7cd06c1a"
Expand All @@ -281,10 +319,10 @@
"name": "stdout",
"output_type": "stream",
"text": [
"VdW Ɛ Gradients tensor([ 14.8422, -57.7803, -5.4171, -1.8189, 74.6179, -53.0473],\n",
"VdW Ɛ Gradients tensor([ -88.9419, -15.8891, 265.0663, 96.5979, -1092.2346],\n",
" dtype=torch.float64)\n",
"VdW σ Gradients tensor([ 1.5717e+01, 4.1386e-04, 1.0395e+01, 1.3762e+01, 1.6823e+01,\n",
" -1.1552e-02], dtype=torch.float64)\n"
"VdW σ Gradients tensor([ 3.1174e+01, 4.4822e+01, 4.4113e+01, 4.4252e+01, -3.8771e-02],\n",
" dtype=torch.float64)\n"
]
}
],
Expand All @@ -300,8 +338,8 @@
"metadata": {
"collapsed": false,
"ExecuteTime": {
"end_time": "2023-10-11T23:34:06.458777Z",
"start_time": "2023-10-11T23:33:54.817906Z"
"end_time": "2023-10-19T21:50:53.889296Z",
"start_time": "2023-10-19T21:50:53.881786Z"
}
},
"id": "dd3ccdfe61a0cd09"
Expand Down
11 changes: 4 additions & 7 deletions smee/mm/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,19 +2,16 @@

from smee.mm._config import GenerateCoordsConfig, MinimizationConfig, SimulationConfig
from smee.mm._mm import generate_system_coords, simulate
from smee.mm._ops import (
GRADIENT_DELTA,
GRADIENT_EXCLUDED_ATTRIBUTES,
compute_ensemble_averages,
)
from smee.mm._ops import compute_ensemble_averages
from smee.mm._reporters import TensorReporter, unpack_frames

__all__ = [
"GRADIENT_EXCLUDED_ATTRIBUTES",
"GRADIENT_DELTA",
"compute_ensemble_averages",
"generate_system_coords",
"simulate",
"GenerateCoordsConfig",
"MinimizationConfig",
"SimulationConfig",
"TensorReporter",
"unpack_frames",
]
9 changes: 5 additions & 4 deletions smee/mm/_mm.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -272,7 +272,8 @@ def _run_simulation(
def simulate(
system: smee.TensorSystem | smee.TensorTopology,
force_field: smee.TensorForceField,
coords_config: "smee.mm.GenerateCoordsConfig",
coords: openmm.unit.Quantity,
box_vectors: openmm.unit.Quantity | None,
equilibrate_configs: list[
typing.Union["smee.mm.MinimizationConfig", "smee.mm.SimulationConfig"]
],
Expand All @@ -284,8 +285,8 @@ def simulate(
Args:
system: The system / topology to simulate.
force_field: The force field to simulate with.
coords_config: The configuration defining how to generate the system
coordinates.
coords: The coordinates [Å] to use for the simulation.
box_vectors: The box vectors [Å] to use for the simulation if periodic.
equilibrate_configs: A list of configurations defining the steps to run for
equilibration. No data will be stored from these simulations.
production_config: The configuration defining the production simulation to run.
Expand All @@ -310,7 +311,7 @@ def simulate(

platform = _get_platform(system.is_periodic)

omm_state = generate_system_coords(system, coords_config)
omm_state = coords, box_vectors

omm_system = smee.converters.convert_to_openmm_system(force_field, system)
omm_topology = smee.converters.convert_to_openmm_topology(system)
Expand Down
Loading