forked from lammpstutorials/lammpstutorials.github.io
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathtutorial02.html
More file actions
331 lines (270 loc) · 17.4 KB
/
Copy pathtutorial02.html
File metadata and controls
331 lines (270 loc) · 17.4 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
<!DOCTYPE html>
<html lang="en">
<head>
<meta charset="UTF-8" />
<meta http-equiv="X-UA-Compatible" content="IE=edge">
<meta name="viewport" content="width=device-width, initial-scale=1, minimum-scale=1.0, shrink-to-fit=no">
<link href="favicon-32x32.png" rel="icon" />
<title>Simple free energy calculation</title>
<meta name="description" content="Simple free energy calculation tutorial">
<meta name="author" content="Simon Gravelle">
<link rel="stylesheet" type="text/css" href="assets_tutorials/vendor/bootstrap/css/bootstrap.min.css" />
<link rel="stylesheet" type="text/css" href="assets_tutorials/vendor/font-awesome/css/all.min.css" /> <!-- Font Awesome Icon -->
<link rel="stylesheet" type="text/css" href="assets_tutorials/css/stylesheet.css" />
<script async src="https://www.googletagmanager.com/gtag/js?id=G-W1WGEC5GQ8"></script>
<script>
window.dataLayer = window.dataLayer || [];
function gtag(){dataLayer.push(arguments);}
gtag('js', new Date());
gtag('config', 'G-W1WGEC5GQ8');
</script>
<!-- Latex formula -->
<script id="MathJax-script" async src="https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js"></script>
</head>
<body data-spy="scroll" data-target=".idocs-navigation" data-offset="125">
<div id="main-wrapper">
<header id="header" class="sticky-top">
<nav class="primary-menu navbar navbar-expand-lg navbar-dropdown-dark">
<div class="container-fluid">
<a class="logo ml-md-3" href="index.html" title="LAMMPS tutorials"> <img src="Figures/logo.png" alt="LAMMPS tutorials" height="60"/> </a>
<button class="navbar-toggler ml-auto" type="button" data-toggle="collapse" data-target="#header-nav"><span></span><span></span><span></span></button>
<div id="header-nav" class="collapse navbar-collapse justify-content-end">
<ul class="navbar-nav">
<li class="dropdown"> <a class="dropdown-toggle" href="#">All tutorials</a>
<ul class="dropdown-menu">
<li><a class="dropdown-item" href="tutorial01.html">Tutorial 01</a></li>
<li><a class="dropdown-item" href="tutorial02.html">Tutorial 02</a></li>
<li><a class="dropdown-item" href="tutorial03.html">Tutorial 03</a></li>
<li><a class="dropdown-item" href="tutorial04.html">Tutorial 04</a></li>
<li><a class="dropdown-item" href="tutorial05.html">Tutorial 05</a></li>
</ul>
</li>
<li><a target="_blank" href="https://simongravelle.github.io/">About me</a></li>
<li><a target="_blank" href="https://www.patreon.com/molecularsimulations">Get help with your script</a></li>
</ul>
</div>
<ul class="social-icons social-icons-sm ml-lg-2 mr-2">
<li class="social-icons-twitter"><a data-toggle="tooltip" href="https://twitter.com/GravelleSimon" target="_blank" title="" data-original-title="Twitter"><i class="fab fa-twitter"></i></a></li>
<li class="social-icons-twitter"><a data-toggle="tooltip" href="https://gitlab.com/sgravelle" target="_blank" title="" data-original-title="Gitlab"><i class="fab fa-gitlab"></i></i></a></li>
<li class="social-icons-twitter"><a data-toggle="tooltip" href="https://github.com/simongravelle" target="_blank" title="" data-original-title="Github"><i class="fab fa-github"></i></i></a></li>
<li class="social-icons-twitter"><a data-toggle="tooltip" href="https://www.youtube.com/channel/UCLmK_9wpyLVpcP7BPgN6BIw?view_as=subscriber" target="_blank" title="" data-original-title="Youtube"><i class="fab fa-youtube"></i></i></a></li>
</ul>
</div>
</nav>
</header>
<div class="idocs-navigation bg-light">
<ul class="nav flex-column ">
<li class="nav-item"><a class="nav-link active" href="#intro">Getting Started</a></li>
<li class="nav-item"><a class="nav-link active" href="#freesampling">Free sampling</a></li>
<li class="nav-item"><a class="nav-link" href="#Umbrellasampling">Umbrella sampling</a></li>
<li class="nav-item"><a class="nav-link" href="#further">Going further</a></li>
</ul>
</div>
<div class="idocs-content">
<div class="container">
<section id="intro">
<h1>Free energy calculation</h1>
<h3>A simple free energy calculation using LAMMPS and WHAM.</h3>
<br>
<p><img src="Figures/tutorial02/largerRes.jpg" width="300" class="img-fluid" alt="Responsive image"></p>
<br>
<p><b>The objective of this tutorial</b> is to measure the free energy profile across a barrier potential using two methods: free sampling and umbrella sampling. For the sake of simplicity and in order to reduce computation time, the barrier potential will be imposed artificially to the atoms, but the procedure is valid for more complex systems. There are two main parts to this tutorial:
<ul>
<li> <a href="#freesampling">Free sampling - </a>First, the free energy will be calculated using the relatively intuitive free sampling method.</li>
<li> <a href="#Umbrellasampling">Umbrella sampling - </a> Second, the more advanced method called umbrella sampling wil be use.</li>
</ul>
</p>
<p class="alert alert-info">Support the creation of material for LAMMPS by subscribing to my <a href="https://www.youtube.com/channel/UCLmK_9wpyLVpcP7BPgN6BIw?view_as=subscriber" target="_blank">youtube</a> channel, or making a small <a href="https://en.tipeee.com/lammps-tutorials" target="_blank">donation here</a>, any amount would be really appreciated and would highly encourage me to improve this page.</p>
<p>If you are new to LAMMPS, I recommend you to follow <a href="tutorial01.html">this simpler tutorial</a> first.
If you have any suggestion about these tutorials, please contact me by email at simon.gravelle at live.fr.</p>
<section id="freesampling">
<div class="container">
<h2>Method 1: Free sampling </h2>
<p>
One way to calculate the free energy profile is to extract the partition function from a classic (unbiased) molecular dynamics simulation, and then to estimate the Gibbs free energy using
\[\Delta G = -RT \ln(p),\]
where \(\Delta G\) is the free energy difference, \(R\) the gas constant, \(T\) the temperature, and \(p\) the partition function.
<br><br>
As an illustration, let us apply this method to an extremely simple configuration that consists in a few particles diffusing in a box in presence of a position-dependent repealing force that make the centre of the box a relatively unfavourable area to explore. Create an input script, and copy the following lines:
<pre><code># Initialization
variable sigma equal 3.405 # Angstrom
variable epsilon equal 0.238 # Kcal/mol
variable U0 equal ${epsilon} # Kcal/mol
variable dlt equal 1.0 # Angstrom
variable x0 equal 30 # Angstrom
units real
atom_style atomic
pair_style lj/cut 10
boundary p p p
# System definition
region myreg block -100 100 -20 20 -20 20
create_box 1 myreg
create_atoms 1 random 6 341341 myreg
# Simulation settings
mass * 39.95
pair_coeff * * ${epsilon} ${sigma}
neigh_modify every 1 delay 4 check yes
</code></pre></p><p>
If you followed <a href="tutorial01.html">tutorial 1</a>, you must be familiar with these commands. The main novelty here is the use of the system of unit 'real', for which energy in kcal/mol, distance in Angstroms, time in femtosecond. I made this choice for practical reason, as the WHAM algorithm we are going to use in the second part of the tutorial automatically assumes the energy to be in kcal/mol. I choose Argon as the gas of interest, which explains the values of the Lennard-Jones parameter \(\sigma\) and \(\epsilon\), as well as the mass \(m = 39.95\) grams/mole. The variables \(U_0\), \(\delta\), and \(x_0\) are used to create the potential. I have chosen it to be of the form:
\[U(x) / U_0 = \arctan \left( \dfrac{x+x_0}{\delta} \right)- \arctan \left( \dfrac{x-x_0}{\delta} \right),\]
which looks like that:
</p>
<p><img src="Figures/tutorial02/potential.jpg" class="img-fluid" alt="Responsive image"></p>
<p>
From the derivative of the potential with respect to \(x\), we obtain the expression for the force that we are going to impose to the particles in the simulation,
\[F(x)=U_0/((x-x_0)^2/\delta^2+1)/\delta-U_0/((x+x_0)^2/\delta^2+1)/\delta,\]
which looks like that:
</p>
<p><img src="Figures/tutorial02/force.jpg" class="img-fluid" alt="Responsive image"></p>
<p>
In order to impose \(F(x)\) to all of the atoms in the simulation, let us use the 'addforce' command, and define an 'atom' variable as follows:</code></pre></p><p>
<pre><code># Run
variable F atom ${U0}/((x-${x0})^2/${dlt}^2+1)/${dlt}-${U0}/((x+${x0})^2/${dlt}^2+1)/${dlt}
fix myadf all addforce v_F 0.0 0.0
</code></pre></p><p>
Finally, let us use the NVE command (constant number of atoms \(N\), constant volume \(V\), and constant energy \(E\)) with a Langevin thermostat. With these two commands, the simulation will be performed in the NVT ensemble (constant number of atoms \(N\), constant volume \(V\), and constant temperature \(T\)). After an equilibration step of 1000000 timestep, the density profile of the atoms along the \(x\) axis is recorded using an 'ave/chunk' command.
<pre><code>fix mynve all nve
fix mylgv all langevin 119.8 119.8 100 1530917
timestep 2.0
thermo 100000
run 1000000
reset_timestep 0
compute cc1 all chunk/atom bin/1d x 0.0 2.0
fix myac all ave/chunk 100 10000000 1000000000 cc1 density/number file density.dat
dump mydmp all atom 100000000 dump.lammpstrj
run 1000000000
</code></pre></p><p>
We can extract the partition function from the density.dat file,
</p>
<p><img src="Figures/tutorial02/partitionfunction.jpg" class="img-fluid" alt="Responsive image"></p>
<p>
from which it is possible to calculate the free energy profile using the formula given at the top of this page (black disk) and compare it to the imposed potential (red line):
</p>
<p><img src="Figures/tutorial02/potential2.jpg" class="img-fluid" alt="Responsive image"><
<p>
The results confirm that we indeed recover the imposed potential \(U(x)\). I did run the simulation five times with different seeds to improve the quality of the curve.
<br><br>
If we increase the value of \(U_0\), the number of times the particles will explore the central region during the simulation will decrease, making it difficult to obtain a good resolution for the free energy profile. This is why it is sometimes necessary to use slightly more evolved methods, such as umbrella sampling, to extract free energy profiles. This is what we are going to do next.
</p>
</div>
</section>
<section id="Umbrellasampling">
<div class="container">
<h2>Method 2: Umbrella sampling </h2>
<br>
<p>Umbrella sampling is a 'biased molecular dynamics' method, i.e a method in which the interactions of the system are modified in order to make the 'unfavourable states' more likely to be explored by the particle. Starting from the same system as previously, we are going to add a potential \(V\) to one of the particle, and force it to move along the axe \(x\). The chosen path is called the axe of reaction. The final simulation will be analysed using the weighted histogram analysis method (WHAM), which allows to remove the effect of the bias and eventually deduce the unbiased free energy profile. In a different folder, create a new script and copy the following lines:
<pre><code># Initialization
variable sigma equal 3.405 # Angstrom
variable epsilon equal 0.238 # Kcal/mol
variable U0 equal ${epsilon} # Kcal/mol
variable dlt equal 1.0 # Angstrom
variable x0 equal 30 # Angstrom
variable k equal 0.0205 # Kcal/mol/Angstrom^2
units real
atom_style atomic
pair_style lj/cut 10
boundary p p p
# System definition
region myreg block -100 100 -20 20 -20 20
create_box 2 myreg
create_atoms 2 single 0 0 0
create_atoms 1 random 5 341341 myreg
# Simulation settings
mass * 39.948
pair_coeff * * ${epsilon} ${sigma}
neigh_modify every 1 delay 4 check yes
group topull type 2
# Run
variable F atom ${U0}/((x-${x0})^2/${dlt}^2+1)/${dlt}-${U0}/((x+${x0})^2/${dlt}^2+1)/${dlt}
fix pot all addforce v_F 0.0 0.0
fix mynve all nve
fix mylgv all langevin 119.8 119.8 100 1530917
timestep 2.0
thermo 100000
run 200000
reset_timestep 0
dump mydmp all atom 1000000 dump.lammpstrj
</code></pre></p><p>
This code resembles the one of Method 1, except that we created one single particle of type 2. This particle is identical to the particles of type 1, and will be the only one to feel the biasing potential. Let us create a loop with 67 steps, aand move progressively the centre of the bias potential by increment of 0.3 nm.
<pre><code>variable a loop 67
label loop
variable xdes equal ${a}*3-102
variable xave equal xcm(topull,x)
fix mytth topull spring tether ${k} ${xdes} 0 0 0
run 200000
fix myat1 all ave/time 10 10 100 v_xave v_xdes file position.${a}.dat
run 50000000
unfix myat1
next a
jump input.lammps loop
</code></pre></p><p>
The spring command serves to impose the additional harmonic potential with spring constant \(k\). The centre of the harmonic potential \(x_\text{des}\) goes from -99 to 99. For each value of \(x_\text{des}\), an equilibration step of 400 ps is performed, followed by a step of 100 ns during which the position along \(x\) of the particle is saved in data files (one data file per value of \(x_\text{des}\)).
<br><br>
In order to treat the data, we are going to use the WHAM algorithm. You can download and compile the version of <a href="http://membrane.urmc.rochester.edu/?page_id=126" target="_blank" class="removelinkdefault">Alan Grossfield</a>. In order to apply the WHAM algorithm to our simulation, we first need to create a metadata file. This file simply contain the paths of the data file, the value of \(x_\text{des}\), and the values of \(k\). To generate the file more easily, you can run this script using Octave or Matlab (assuming that the wham algorithm is located in the same folder as the LAMMPS simulations)
<pre><code>file=fopen('metadata.dat','wt');
for a=1:67
X=['./position.',num2str(a),'.dat ',num2str(a*3-102),' 0.0205'];
fprintf(file,X);
fprintf(file,'\n');
end
</code></pre></p><p>
The generated file named metadata.dat looks like that (alternatively you can download it <a href="Data/tutorial02/metadata.dat" target="_blank" class="removelinkdefault">here</a>):
<pre><code>./position.1.dat -99 0.0205
./position.2.dat -96 0.0205
./position.3.dat -93 0.0205
./position.4.dat -90 0.0205
./position.5.dat -87 0.0205
(...)
./position.66.dat 96 0.0205
./position.67.dat 99 0.0205
</code></pre></p><p>
Then, simply run the following command
<pre><code>./wham -100 100 160 1e-8 119.8 0 metadata.dat PMF.dat
</code></pre></p><p>
where -100 and 100 are the boundaries, 160 the number of bins, 1e-8 the tolerance, and 119.8 the temperature. A file named PMF.dat has been created, and contains the free energy profile in Kcal/mol. We can compare the PMF we the imposed potential, and the agreement in again quite good:
</p>
<p><img src="Figures/tutorial02/PMF.jpg" class="img-fluid" alt="Responsive image"></p>
<p>
Again, five independent simulations where run to improve the curve. </p>
</div>
</section>
<br><br><br><section id="further">
<h2>Going further</h2>
<p><b>Realistic system.</b> You can apply the same two methods to a more realistic and interesting system. For example, you can measure the free energy profile normal to a water-vapour interface, or normal to a solid liquid interface.
<br><br>
<p><b>Channel selectivity.</b> Umbrella sampling method can be used to evaluate the barrier potential that a channel, for example a protein, oppose to a molecule or an ion. In this case, the procedure is to force the molecule/ion to cross the channel, and then reconstruct the free energy profile as we did here.
<br><br>
<p class="alert alert-info">Support the creation of material for LAMMPS by subscribing to my <a href="https://www.youtube.com/channel/UCLmK_9wpyLVpcP7BPgN6BIw?view_as=subscriber" target="_blank">youtube</a> channel, or making a small <a href="https://en.tipeee.com/lammps-tutorials" target="_blank">donation here</a>, any amount would be really appreciated and would highly encourage me to improve this page.</p>
</div>
</section>
</section>
</div>
</div>
</div>
<!-- Content end -->
<!-- Footer
============================ -->
<footer id="footer" class="section">
<div class="container">
<p class="text-2 text-center mb-0">This template has been adapted from <a class="btn-link" target="_blank" href="http://www.harnishdesign.net/">HarnishDesign</a>.</p>
</div>
</footer>
<!-- Footer end -->
</div>
<!-- Document Wrapper end -->
<!-- Back To Top -->
<a id="back-to-top" data-toggle="tooltip" title="Back to Top" href="javascript:void(0)"><i class="fa fa-chevron-up"></i></a>
<!-- JavaScript
============================ -->
<script src="assets_tutorials/vendor/jquery/jquery.min.js"></script>
<script src="assets_tutorials/vendor/bootstrap/js/bootstrap.bundle.min.js"></script>
<!-- Highlight JS -->
<script src="assets_tutorials/vendor/highlight.js/highlight.min.js"></script>
<!-- Easing -->
<script src="assets_tutorials/vendor/jquery.easing/jquery.easing.min.js"></script>
<!-- Magnific Popup -->
<script src="assets_tutorials/vendor/magnific-popup/jquery.magnific-popup.min.js"></script>
<!-- Custom Script -->
<script src="assets_tutorials/js/theme.js"></script>
</body>
</html>