Skip to content

Commit 0562cfb

Browse files
jthortonjthorton
committed
add trajectory notebook
1 parent 5eb2a12 commit 0562cfb

1 file changed

Lines changed: 255 additions & 0 deletions

File tree

‎simulation_visualisation/simulation_visualisation.ipynb‎

Lines changed: 255 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,255 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"id": "20a8d98a",
6+
"metadata": {},
7+
"source": [
8+
"# Extracting and visualising a free energy simulation\n",
9+
"\n",
10+
"This notebook provides a step-by-step guide to extract and visualise a free energy simulation trajectory from a ``simulation.nc`` file using [openfe-analysis](https://github.com/OpenFreeEnergy/openfe_analysis) and [MDAnalysis](https://github.com/MDAnalysis/mdanalysis). By the end, you should understand how to:\n",
11+
"\n",
12+
"1. Extract the trajectory of a ``replica`` or ``single lambda state`` from a ``simulation.nc`` file\n",
13+
"2. Centre the ligand in the simulation box using `MDAnalysis`\n",
14+
"3. Write out the trajectorie(s) using `MDAnalysis`\n",
15+
"\n",
16+
"## Downloading the example data\n",
17+
"\n",
18+
"First, download some example trajectory data. This may take a few minutes due to the size of the simulation file. Please skip this section if you have already done this!"
19+
]
20+
},
21+
{
22+
"cell_type": "code",
23+
"execution_count": 1,
24+
"id": "be6317291bbf3804",
25+
"metadata": {},
26+
"outputs": [
27+
{
28+
"name": "stdout",
29+
"output_type": "stream",
30+
"text": [
31+
"--2025-09-05 10:31:50-- https://zenodo.org/records/15375081/files/simulation.nc\n",
32+
"Resolving zenodo.org (zenodo.org)... 2001:1458:d00:25::100:372, 2001:1458:d00:24::100:f6, 2001:1458:d00:61::100:2f3, ...\n",
33+
"Connecting to zenodo.org (zenodo.org)|2001:1458:d00:25::100:372|:443... connected.\n",
34+
"HTTP request sent, awaiting response... 200 OK\n",
35+
"Length: 516886878 (493M) [application/octet-stream]\n",
36+
"Saving to: ‘simulation.nc.1’\n",
37+
"\n",
38+
"simulation.nc.1 100%[===================>] 492.94M 12.2MB/s in 62s \n",
39+
"\n",
40+
"2025-09-05 10:32:53 (7.96 MB/s) - ‘simulation.nc.1’ saved [516886878/516886878]\n",
41+
"\n",
42+
"--2025-09-05 10:32:53-- https://zenodo.org/records/15375081/files/hybrid_system.pdb\n",
43+
"Resolving zenodo.org (zenodo.org)... 2001:1458:d00:24::100:f6, 2001:1458:d00:61::100:2f3, 2001:1458:d00:25::100:372, ...\n",
44+
"Connecting to zenodo.org (zenodo.org)|2001:1458:d00:24::100:f6|:443... connected.\n",
45+
"HTTP request sent, awaiting response... 200 OK\n",
46+
"Length: 388547 (379K) [application/octet-stream]\n",
47+
"Saving to: ‘hybrid_system.pdb.1’\n",
48+
"\n",
49+
"hybrid_system.pdb.1 100%[===================>] 379.44K 1.65MB/s in 0.2s \n",
50+
"\n",
51+
"2025-09-05 10:32:54 (1.65 MB/s) - ‘hybrid_system.pdb.1’ saved [388547/388547]\n",
52+
"\n"
53+
]
54+
}
55+
],
56+
"source": [
57+
"! wget https://zenodo.org/records/15375081/files/simulation.nc\n",
58+
"! wget https://zenodo.org/records/15375081/files/hybrid_system.pdb"
59+
]
60+
},
61+
{
62+
"cell_type": "markdown",
63+
"id": "5d066221-0b8e-4b1b-a047-c989a5957cf7",
64+
"metadata": {},
65+
"source": [
66+
"## Extracting the trajectory with `MDAnalysis`\n",
67+
"\n",
68+
"The `openfe-analysis` package provides an `MDAnalysis` reader to help extract the trajectory data from the `simulation.nc` file. As the file contains multipule replicas simulated at different lambda states, we must choose which of these to load as a single trajectory. We have two options available to construct the trajectory:\n",
69+
"- `state_id`: will construct a trajectory which follows a single Hamiltonian lambda state at the specified value.\n",
70+
"- `recplica_id`: will construct a trajectory which follows a single replica at the specified value.\n",
71+
"\n",
72+
"In this example which uses a trajectory from a relative binding free energy calculation we will load the trajectory at `lambda=0` or the end state corresponding to Ligand A and visulaise the trajectory with `nglview`."
73+
]
74+
},
75+
{
76+
"cell_type": "code",
77+
"execution_count": 6,
78+
"id": "05ba7dc7",
79+
"metadata": {},
80+
"outputs": [
81+
{
82+
"name": "stderr",
83+
"output_type": "stream",
84+
"text": [
85+
"/Users/joshua/mambaforge/envs/openfe_fast_set/lib/python3.12/site-packages/openfe_analysis/utils/multistate.py:41: UserWarning: This is an older NetCDF file that does not yet contain information about the write frequency of positions and velocities. We will assume that positions and velocities were written out at every iteration. \n",
86+
" warnings.warn(wmsg)\n"
87+
]
88+
},
89+
{
90+
"data": {
91+
"application/vnd.jupyter.widget-view+json": {
92+
"model_id": "ca7edee8791e482fa36a171ee3edb895",
93+
"version_major": 2,
94+
"version_minor": 0
95+
},
96+
"text/plain": [
97+
"NGLWidget(max_frame=500)"
98+
]
99+
},
100+
"metadata": {},
101+
"output_type": "display_data"
102+
}
103+
],
104+
"source": [
105+
"import MDAnalysis as mda\n",
106+
"from openfe_analysis import FEReader\n",
107+
"import nglview as nv\n",
108+
"\n",
109+
"u_0 = mda.Universe(\"hybrid_system.pdb\", \"simulation.nc\", format=FEReader, state_id=0)\n",
110+
"\n",
111+
"w = nv.show_mdanalysis(u_0)\n",
112+
"w"
113+
]
114+
},
115+
{
116+
"cell_type": "markdown",
117+
"id": "f0e7b0e7-04b5-42c0-ab65-2202275288d4",
118+
"metadata": {},
119+
"source": [
120+
"<div class=\\\"alert alert-block alert-info\\\"> <b>Note:</b> The OpenFE relative binding free energy protocol does not save water positions by default, this can be changed via the <a href=\"https://docs.openfree.energy/en/latest/reference/api/openmm_protocol_settings.html#openfe.protocols.openmm_utils.omm_settings.MultiStateOutputSettings.output_indices\">output_indices</a> protocol setting. </div>\n",
121+
"\n",
122+
"\n",
123+
"To view the final state at `lambda=1` we can use negative indexing if we don't know the total number of lambda states."
124+
]
125+
},
126+
{
127+
"cell_type": "code",
128+
"execution_count": 7,
129+
"id": "009a8ac7-0ca9-494c-bbd6-177e752b0d4f",
130+
"metadata": {},
131+
"outputs": [
132+
{
133+
"name": "stderr",
134+
"output_type": "stream",
135+
"text": [
136+
"/Users/joshua/mambaforge/envs/openfe_fast_set/lib/python3.12/site-packages/openfe_analysis/utils/multistate.py:41: UserWarning: This is an older NetCDF file that does not yet contain information about the write frequency of positions and velocities. We will assume that positions and velocities were written out at every iteration. \n",
137+
" warnings.warn(wmsg)\n",
138+
"DEPRECATED: Please use 'center' method\n"
139+
]
140+
},
141+
{
142+
"data": {
143+
"application/vnd.jupyter.widget-view+json": {
144+
"model_id": "18e16cf1a1144b228e6ab22a06d42002",
145+
"version_major": 2,
146+
"version_minor": 0
147+
},
148+
"text/plain": [
149+
"NGLWidget(max_frame=500)"
150+
]
151+
},
152+
"metadata": {},
153+
"output_type": "display_data"
154+
}
155+
],
156+
"source": [
157+
"u_1 = mda.Universe(\"hybrid_system.pdb\", \"simulation.nc\", format=FEReader, state_id=-1)\n",
158+
"\n",
159+
"w = nv.show_mdanalysis(u_1)\n",
160+
"w.center_view()\n",
161+
"w"
162+
]
163+
},
164+
{
165+
"cell_type": "markdown",
166+
"id": "75602d58-6709-4d9b-8a52-9896fdc6a0af",
167+
"metadata": {},
168+
"source": [
169+
"## Centring the Ligand\n",
170+
"\n",
171+
"You may have noticed in the view above that the ligand seems to have drifted away from the protein, this is a visualisation artifact caused by the use of periodic boundary conditions and the way in which OpenMM tries to ensure that all particle positions are written into a single periodic box. We can fix this however, using transformations in `MDAnalysis` to centre the ligand and wrap the protein around it like so:"
172+
]
173+
},
174+
{
175+
"cell_type": "code",
176+
"execution_count": 8,
177+
"id": "a5f16795-1b9c-4198-8694-6568eaba06c7",
178+
"metadata": {},
179+
"outputs": [
180+
{
181+
"data": {
182+
"application/vnd.jupyter.widget-view+json": {
183+
"model_id": "2f9fe1601d5a4c22954bdc2dd69ad482",
184+
"version_major": 2,
185+
"version_minor": 0
186+
},
187+
"text/plain": [
188+
"NGLWidget(max_frame=500)"
189+
]
190+
},
191+
"metadata": {},
192+
"output_type": "display_data"
193+
}
194+
],
195+
"source": [
196+
"import MDAnalysis.transformations as trans\n",
197+
"\n",
198+
"ligand = u_0.select_atoms(\"resname UNK\")\n",
199+
"protein = u_0.select_atoms(\"protein\")\n",
200+
"transforms = [\n",
201+
" trans.center_in_box(ligand, wrap=True),\n",
202+
" trans.wrap(protein),\n",
203+
" trans.nojump.NoJump()\n",
204+
"]\n",
205+
"u_0.trajectory.add_transformations(*transforms)\n",
206+
"\n",
207+
"w = nv.show_mdanalysis(u_0)\n",
208+
"w"
209+
]
210+
},
211+
{
212+
"cell_type": "markdown",
213+
"id": "7026738c-e4fe-461e-85ca-f726e9768a5e",
214+
"metadata": {},
215+
"source": [
216+
"## Saving the trajectory to file\n",
217+
"\n",
218+
"We can now use `MDAnalysis` to save our trajectory to a common file format:"
219+
]
220+
},
221+
{
222+
"cell_type": "code",
223+
"execution_count": 9,
224+
"id": "6e410679",
225+
"metadata": {},
226+
"outputs": [],
227+
"source": [
228+
"with mda.Writer('out.xtc', n_atoms=len(u_0.atoms)) as w:\n",
229+
" for ts in u_0.trajectory:\n",
230+
" w.write(u_0.atoms)"
231+
]
232+
}
233+
],
234+
"metadata": {
235+
"kernelspec": {
236+
"display_name": "Python 3 (ipykernel)",
237+
"language": "python",
238+
"name": "python3"
239+
},
240+
"language_info": {
241+
"codemirror_mode": {
242+
"name": "ipython",
243+
"version": 3
244+
},
245+
"file_extension": ".py",
246+
"mimetype": "text/x-python",
247+
"name": "python",
248+
"nbconvert_exporter": "python",
249+
"pygments_lexer": "ipython3",
250+
"version": "3.12.11"
251+
}
252+
},
253+
"nbformat": 4,
254+
"nbformat_minor": 5
255+
}

0 commit comments

Comments
 (0)