Skip to content

Commit c140587

Browse files
p-zubietap-zubieta
authored
Pair latest notebooks to md (#183)
1 parent 0cdf384 commit c140587

3 files changed

Lines changed: 385 additions & 2 deletions

File tree

examples/hoomd-blue/spectral_abf/butane_spectralABF.ipynb renamed to examples/hoomd-blue/spectral_abf/Butane-SpectralABF.ipynb

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -705,8 +705,7 @@
705705
"provenance": []
706706
},
707707
"jupytext": {
708-
"formats": "ipynb,md",
709-
"main_language": "python"
708+
"formats": "ipynb,md"
710709
},
711710
"kernelspec": {
712711
"display_name": "Python 3",

‎examples/hoomd-blue/spectral_abf/Butane-SpectralABF.md‎

Lines changed: 384 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,384 @@
1+
---
2+
jupyter:
3+
jupytext:
4+
formats: ipynb,md
5+
text_representation:
6+
extension: .md
7+
format_name: markdown
8+
format_version: '1.3'
9+
jupytext_version: 1.14.0
10+
kernelspec:
11+
display_name: Python 3
12+
language: python
13+
name: python3
14+
---
15+
16+
<!-- #region id="T-Qkg9C9n7Cc" -->
17+
# Setting up the environment
18+
19+
First, we are setting up our environment. We use an already compiled and packaged installation of HOOMD-blue and the DLExt plugin.
20+
We copy it from Google Drive and install PySAGES for it.
21+
<!-- #endregion -->
22+
23+
```bash id="3eTbKklCnyd_"
24+
25+
BASE_URL="https://drive.google.com/u/0/uc?id=1hsKkKtdxZTVfHKgqVF6qV2e-4SShmhr7&export=download"
26+
wget -q --load-cookies /tmp/cookies.txt "$BASE_URL&confirm=$(wget -q --save-cookies /tmp/cookies.txt --keep-session-cookies --no-check-certificate $BASE_URL -O- | sed -rn 's/.*confirm=(\w+).*/\1\n/p')" -O pysages-env.zip
27+
rm -rf /tmp/cookies.txt
28+
```
29+
30+
```python colab={"base_uri": "https://localhost:8080/"} id="KRPmkpd9n_NG" outputId="acc3a92c-182f-415b-d8dc-b5af076b3d01"
31+
%env PYSAGES_ENV=/env/pysages
32+
```
33+
34+
```bash id="J7OY5K9VoBBh"
35+
36+
mkdir -p $PYSAGES_ENV .
37+
unzip -qquo pysages-env.zip -d $PYSAGES_ENV
38+
```
39+
40+
```python id="LlVSU_-FoD4w"
41+
!update-alternatives --auto libcudnn &> /dev/null
42+
```
43+
44+
```python id="EMAWp8VloIk4"
45+
import os
46+
import sys
47+
48+
ver = sys.version_info
49+
sys.path.append(os.environ["PYSAGES_ENV"] + "/lib/python" + str(ver.major) + "." + str(ver.minor) + "/site-packages/")
50+
51+
os.environ["XLA_FLAGS"] = "--xla_gpu_strict_conv_algorithm_picker=false"
52+
os.environ["LD_LIBRARY_PATH"] = "/usr/lib/x86_64-linux-gnu:" + os.environ["LD_LIBRARY_PATH"]
53+
```
54+
55+
<!-- #region id="we_mTkFioS6R" -->
56+
## PySAGES
57+
58+
The next step is to install PySAGES.
59+
First, we install the jaxlib version that matches the CUDA installation of this Colab setup. See the JAX documentation [here](https://github.com/google/jax) for more details.
60+
<!-- #endregion -->
61+
62+
```bash id="vK0RZtbroQWe"
63+
64+
pip install -q --upgrade pip
65+
# Installs the wheel compatible with CUDA 11 and cuDNN 8.0.5.
66+
pip install -q --upgrade "jax[cuda11_cudnn805]" -f https://storage.googleapis.com/jax-releases/jax_cuda_releases.html &> /dev/null
67+
```
68+
69+
<!-- #region id="wAtjM-IroYX8" -->
70+
Now we can finally install PySAGES. We clone the newest version from [here](https://github.com/SSAGESLabs/PySAGES) and build the remaining pure python dependencies and PySAGES itself.
71+
<!-- #endregion -->
72+
73+
```bash id="B-HB9CzioV5j"
74+
75+
rm -rf PySAGES
76+
git clone https://github.com/SSAGESLabs/PySAGES.git &> /dev/null
77+
cd PySAGES
78+
pip install -q . &> /dev/null
79+
```
80+
81+
```bash id="ppTzMmyyobHB"
82+
83+
mkdir /content/cff
84+
cd /content/cff
85+
```
86+
87+
<!-- #region id="KBFVcG1FoeMq" -->
88+
# SpectralABF-biased simulations
89+
<!-- #endregion -->
90+
91+
<!-- #region id="0W2ukJuuojAl" -->
92+
SpectralABF gradually learns a better approximation to the coefficients of a basis functions expansion of the free energy of a system, from the generalized mean forces in a similar fashion to the ABF sampling method.
93+
94+
For this Colab, we are using butane as the example molecule.
95+
<!-- #endregion -->
96+
97+
```python id="BBvC7Spoog82"
98+
import hoomd
99+
import hoomd.md
100+
101+
import numpy
102+
103+
104+
pi = numpy.pi
105+
kT = 0.596161
106+
dt = 0.02045
107+
mode = "--mode=gpu"
108+
109+
110+
def generate_context(kT = kT, dt = dt, mode = mode):
111+
"""
112+
Generates a simulation context, we pass this function to the attribute
113+
`run` of our sampling method.
114+
"""
115+
hoomd.context.initialize(mode)
116+
117+
### System Definition
118+
snapshot = hoomd.data.make_snapshot(
119+
N = 14,
120+
box = hoomd.data.boxdim(Lx = 41, Ly = 41, Lz = 41),
121+
particle_types = ['C', 'H'],
122+
bond_types = ["CC", "CH"],
123+
angle_types = ["CCC", "CCH", "HCH"],
124+
dihedral_types = ["CCCC", "HCCC", "HCCH"],
125+
pair_types = ["CCCC", "HCCC", "HCCH"],
126+
dtype = "double"
127+
)
128+
129+
snapshot.particles.typeid[0] = 0
130+
snapshot.particles.typeid[1:4] = 1
131+
snapshot.particles.typeid[4] = 0
132+
snapshot.particles.typeid[5:7] = 1
133+
snapshot.particles.typeid[7] = 0
134+
snapshot.particles.typeid[8:10] = 1
135+
snapshot.particles.typeid[10] = 0
136+
snapshot.particles.typeid[11:14] = 1
137+
138+
positions = numpy.array([
139+
[-2.990196, 0.097881, 0.000091],
140+
[-2.634894, -0.911406, 0.001002],
141+
[-2.632173, 0.601251, -0.873601],
142+
[-4.060195, 0.099327, -0.000736],
143+
[-2.476854, 0.823942, 1.257436],
144+
[-2.832157, 1.833228, 1.256526],
145+
[-2.834877, 0.320572, 2.131128],
146+
[-0.936856, 0.821861, 1.258628],
147+
[-0.578833, 1.325231, 0.384935],
148+
[-0.581553, -0.187426, 1.259538],
149+
[-0.423514, 1.547922, 2.515972],
150+
[-0.781537, 1.044552, 3.389664],
151+
[ 0.646485, 1.546476, 2.516800],
152+
[-0.778816, 2.557208, 2.515062]
153+
])
154+
155+
reference_box_low_coords = numpy.array([-22.206855, -19.677099, -19.241968])
156+
box_low_coords = numpy.array([
157+
-snapshot.box.Lx / 2,
158+
-snapshot.box.Ly / 2,
159+
-snapshot.box.Lz / 2
160+
])
161+
positions += (box_low_coords - reference_box_low_coords)
162+
163+
snapshot.particles.position[:] = positions[:]
164+
165+
mC = 12.00
166+
mH = 1.008
167+
snapshot.particles.mass[:] = [
168+
mC, mH, mH, mH,
169+
mC, mH, mH,
170+
mC, mH, mH,
171+
mC, mH, mH, mH
172+
]
173+
174+
reference_charges = numpy.array([
175+
-0.180000, 0.060000, 0.060000, 0.060000,
176+
-0.120000, 0.060000, 0.060000,
177+
-0.120000, 0.060000, 0.060000,
178+
-0.180000, 0.060000, 0.060000, 0.060000]
179+
)
180+
charge_conversion = 18.22262
181+
snapshot.particles.charge[:] = charge_conversion * reference_charges[:]
182+
183+
snapshot.bonds.resize(13)
184+
snapshot.bonds.typeid[0:3] = 1
185+
snapshot.bonds.typeid[3] = 0
186+
snapshot.bonds.typeid[4:6] = 1
187+
snapshot.bonds.typeid[6] = 0
188+
snapshot.bonds.typeid[7:9] = 1
189+
snapshot.bonds.typeid[9] = 0
190+
snapshot.bonds.typeid[10:13] = 1
191+
192+
snapshot.bonds.group[:] = [
193+
[0, 2], [0, 1], [0, 3], [0, 4],
194+
[4, 5], [4, 6], [4, 7],
195+
[7, 8], [7, 9], [7, 10],
196+
[10, 11], [10, 12], [10, 13]
197+
]
198+
199+
snapshot.angles.resize(24)
200+
snapshot.angles.typeid[0:2] = 2
201+
snapshot.angles.typeid[2] = 1
202+
snapshot.angles.typeid[3] = 2
203+
snapshot.angles.typeid[4:8] = 1
204+
snapshot.angles.typeid[8] = 0
205+
snapshot.angles.typeid[9] = 2
206+
snapshot.angles.typeid[10:14] = 1
207+
snapshot.angles.typeid[14] = 0
208+
snapshot.angles.typeid[15] = 2
209+
snapshot.angles.typeid[16:21] = 1
210+
snapshot.angles.typeid[21:24] = 2
211+
212+
snapshot.angles.group[:] = [
213+
[1, 0, 2], [2, 0, 3], [2, 0, 4],
214+
[1, 0, 3], [1, 0, 4], [3, 0, 4],
215+
[0, 4, 5], [0, 4, 6], [0, 4, 7],
216+
[5, 4, 6], [5, 4, 7], [6, 4, 7],
217+
[4, 7, 8], [4, 7, 9], [4, 7, 10],
218+
[8, 7, 9], [8, 7, 10], [9, 7, 10],
219+
[7, 10, 11], [7, 10, 12], [7, 10, 13],
220+
[11, 10, 12], [11, 10, 13], [12, 10, 13]
221+
]
222+
223+
snapshot.dihedrals.resize(27)
224+
snapshot.dihedrals.typeid[0:2] = 2
225+
snapshot.dihedrals.typeid[2] = 1
226+
snapshot.dihedrals.typeid[3:5] = 2
227+
snapshot.dihedrals.typeid[5] = 1
228+
snapshot.dihedrals.typeid[6:8] = 2
229+
snapshot.dihedrals.typeid[8:11] = 1
230+
snapshot.dihedrals.typeid[11] = 0
231+
snapshot.dihedrals.typeid[12:14] = 2
232+
snapshot.dihedrals.typeid[14] = 1
233+
snapshot.dihedrals.typeid[15:17] = 2
234+
snapshot.dihedrals.typeid[17:21] = 1
235+
snapshot.dihedrals.typeid[21:27] = 2
236+
237+
snapshot.dihedrals.group[:] = [
238+
[2, 0, 4, 5], [2, 0, 4, 6], [2, 0, 4, 7],
239+
[1, 0, 4, 5], [1, 0, 4, 6], [1, 0, 4, 7],
240+
[3, 0, 4, 5], [3, 0, 4, 6], [3, 0, 4, 7],
241+
[0, 4, 7, 8], [0, 4, 7, 9], [0, 4, 7, 10],
242+
[5, 4, 7, 8], [5, 4, 7, 9], [5, 4, 7, 10],
243+
[6, 4, 7, 8], [6, 4, 7, 9], [6, 4, 7, 10],
244+
[4, 7, 10, 11], [4, 7, 10, 12], [4, 7, 10, 13],
245+
[8, 7, 10, 11], [8, 7, 10, 12], [8, 7, 10, 13],
246+
[9, 7, 10, 11], [9, 7, 10, 12], [9, 7, 10, 13]
247+
]
248+
249+
snapshot.pairs.resize(27)
250+
snapshot.pairs.typeid[0:1] = 0
251+
snapshot.pairs.typeid[1:11] = 1
252+
snapshot.pairs.typeid[11:27] = 2
253+
snapshot.pairs.group[:] = [
254+
# CCCC
255+
[0, 10],
256+
# HCCC
257+
[0, 8], [0, 9], [5, 10], [6, 10],
258+
[1, 7], [2, 7], [3, 7],
259+
[11, 4], [12, 4], [13, 4],
260+
# HCCH
261+
[1, 5], [1, 6], [2, 5], [2, 6], [3, 5], [3, 6],
262+
[5, 8], [6, 8], [5, 9], [6, 9],
263+
[8, 11], [8, 12], [8, 13], [9, 11], [9, 12], [9, 13]
264+
]
265+
266+
hoomd.init.read_snapshot(snapshot)
267+
268+
### Set interactions
269+
nl_ex = hoomd.md.nlist.cell()
270+
nl_ex.reset_exclusions(exclusions = ["1-2", "1-3", "1-4"])
271+
272+
lj = hoomd.md.pair.lj(r_cut = 12.0, nlist = nl_ex)
273+
lj.pair_coeff.set('C', 'C', epsilon = 0.07, sigma = 3.55)
274+
lj.pair_coeff.set('H', 'H', epsilon = 0.03, sigma = 2.42)
275+
lj.pair_coeff.set('C', 'H', epsilon = numpy.sqrt(0.07*0.03), sigma = numpy.sqrt(3.55*2.42))
276+
277+
coulomb = hoomd.md.charge.pppm(hoomd.group.charged(), nlist = nl_ex)
278+
coulomb.set_params(Nx = 64, Ny = 64, Nz = 64, order = 6, rcut = 12.0)
279+
280+
harmonic = hoomd.md.bond.harmonic()
281+
harmonic.bond_coeff.set("CC", k = 2*268.0, r0 = 1.529)
282+
harmonic.bond_coeff.set("CH", k = 2*340.0, r0 = 1.09)
283+
284+
angle = hoomd.md.angle.harmonic()
285+
angle.angle_coeff.set("CCC", k = 2*58.35, t0 = 112.7 * pi / 180)
286+
angle.angle_coeff.set("CCH", k = 2*37.5, t0 = 110.7 * pi / 180)
287+
angle.angle_coeff.set("HCH", k = 2*33.0, t0 = 107.8 * pi / 180)
288+
289+
290+
dihedral = hoomd.md.dihedral.opls()
291+
dihedral.dihedral_coeff.set("CCCC", k1 = 1.3, k2 = -0.05, k3 = 0.2, k4 = 0.0)
292+
dihedral.dihedral_coeff.set("HCCC", k1 = 0.0, k2 = 0.0, k3 = 0.3, k4 = 0.0)
293+
dihedral.dihedral_coeff.set("HCCH", k1 = 0.0, k2 = 0.0, k3 = 0.3, k4 = 0.0)
294+
295+
lj_special_pairs = hoomd.md.special_pair.lj()
296+
lj_special_pairs.pair_coeff.set("CCCC", epsilon = 0.07, sigma = 3.55, r_cut = 12.0)
297+
lj_special_pairs.pair_coeff.set("HCCH", epsilon = 0.03, sigma = 2.42, r_cut = 12.0)
298+
lj_special_pairs.pair_coeff.set("HCCC",
299+
epsilon = numpy.sqrt(0.07 * 0.03), sigma = numpy.sqrt(3.55 * 2.42), r_cut = 12.0
300+
)
301+
302+
coulomb_special_pairs = hoomd.md.special_pair.coulomb()
303+
coulomb_special_pairs.pair_coeff.set("CCCC", alpha = 0.5, r_cut = 12.0)
304+
coulomb_special_pairs.pair_coeff.set("HCCC", alpha = 0.5, r_cut = 12.0)
305+
coulomb_special_pairs.pair_coeff.set("HCCH", alpha = 0.5, r_cut = 12.0)
306+
307+
hoomd.md.integrate.mode_standard(dt = dt)
308+
integrator = hoomd.md.integrate.nvt(group = hoomd.group.all(), kT = kT, tau = 100*dt)
309+
integrator.randomize_velocities(seed = 42)
310+
311+
return hoomd.context.current
312+
```
313+
314+
<!-- #region id="3UrzENm_oo6U" -->
315+
Next, we load PySAGES and the relevant classes and methods for our problem
316+
<!-- #endregion -->
317+
318+
```python id="fpMg-o8WomAA"
319+
from pysages.grids import Grid
320+
from pysages.colvars import DihedralAngle
321+
from pysages.methods import SpectralABF
322+
323+
import pysages
324+
```
325+
326+
<!-- #region id="LknkRvo1o4av" -->
327+
The next step is to define the collective variable (CV). In this case, we choose the central dihedral angle.
328+
329+
We define a grid, which will be used to indicate how we want to bin the forces that will be used to approximate the biasing potential and its gradient.
330+
<!-- #endregion -->
331+
332+
```python id="B1Z8FWz0o7u_"
333+
cvs = [DihedralAngle([0, 4, 7, 10])]
334+
grid = Grid(lower=(-pi,), upper=(pi,), shape=(64,), periodic=True)
335+
timesteps = int(5e5)
336+
337+
method = SpectralABF(cvs, grid)
338+
```
339+
340+
<!-- #region id="Fz8BfU34pA_N" -->
341+
We now simulate $5\times10^5$ time steps.
342+
Make sure to run with GPU support, otherwise, it can take a very long time.
343+
<!-- #endregion -->
344+
345+
```python colab={"base_uri": "https://localhost:8080/"} id="K951m4BbpUar" outputId="8005b8a9-2967-4eb9-f9db-e0dc0d523835"
346+
run_result = pysages.run(method, generate_context, timesteps)
347+
```
348+
349+
<!-- #region id="26zdu6yAht5Y" -->
350+
## Analysis
351+
352+
PySAGES provides an `analyze` method that makes it easier to get the free energy of different simulation runs.
353+
<!-- #endregion -->
354+
355+
```python id="2NWmahlfhoj8"
356+
result = pysages.analyze(run_result)
357+
```
358+
359+
<!-- #region id="PXBKUfK0p9T2" -->
360+
Let's plot now the free energy!
361+
<!-- #endregion -->
362+
363+
```python id="X69d1R7OpW4P"
364+
import matplotlib.pyplot as plt
365+
```
366+
367+
```python id="pTIGVSSqKdbs"
368+
mesh = result["mesh"]
369+
A = result["free_energy"]
370+
# Alternatively:
371+
# fes_fn = result["fes_fn"]
372+
# A = fes_fn(mesh)
373+
A = A.max() - A
374+
```
375+
376+
```python colab={"base_uri": "https://localhost:8080/", "height": 302} id="7_d_XfVLLkbI" outputId="e35db259-31f8-4a3b-b1fa-7e91a8a5c88a"
377+
fig, ax = plt.subplots()
378+
379+
ax.set_xlabel(r"Dihedral Angle, $\xi$")
380+
ax.set_ylabel(r"$A(\xi)$")
381+
382+
ax.plot(mesh, A)
383+
plt.gca()
384+
```

examples/hoomd-blue/spectral_abf/butane_spectralABF.py renamed to examples/hoomd-blue/spectral_abf/butane.py

File renamed without changes.

0 commit comments

Comments
 (0)