forked from lammpstutorials/lammpstutorials.github.io
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathtutorial05.html
More file actions
378 lines (308 loc) · 23.5 KB
/
Copy pathtutorial05.html
File metadata and controls
378 lines (308 loc) · 23.5 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
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
<!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>GCMC water silica</title>
<meta name="description" content="GCMC water silica">
<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="#Generation">System generation</a></li>
<li class="nav-item"><a class="nav-link" href="#Cracking">Cracking</a></li>
<li class="nav-item"><a class="nav-link" href="#Water">Adding water</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>Water adsorption in silica</h1>
<h3>Adsorption of water in a silica crack using the grand canonical Monte Carlo method</h3>
<br>
<p><img src="Figures/tutorial05/3silica.png" class="img-fluid" alt="Responsive image">
<small style="color:#6B6B6B">Figure: From left to right, amorphous silica, amorphous silica with crack, cracked amorphous silica with adsorbed water molecule. Silicon (Si) and oxygen (O) atoms of the silica are in yellow and red, respectively, and represented at a network of bonds. Oxygen (O) and hydrogen (O) atoms of the water molecules are in red and white, respectively, and represented at spheres.
</p>
<p><b>The objective of this tutorial</b> is to combine molecular dynamics and grand canonical Monte Carlo to simulate the adsorption of water molecules in a cracked silica, see the figure above. There are three main parts to this tutorial:
<ul>
<li> <a href="#Generation">System generation - </a> First, amorphous silica is generated by temperature annealing.</li>
<li> <a href="#Cracking">Cracking - </a>Seconds, the silica is deformed by dilatation of the box until a crack forms.</li>
<li> <a href="#Water">Adding water - </a> And third, the system is equilibrated with a water reservoir of given chemical potential using grand canonical Monte Carlo.</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>
<br><br><br><section id="Generation">
<h2>Generation of the silica block</h2>
<p>Let us generate a block of amorphous silica (SiO2). To do so, we are going to replicate a building block containing 3 Si and 6 O atoms. The data file for the SiO atoms can be downloaded by clicking <a href="Data/tutorial05/dataSiO.data" target="_blank" class="removelinkdefault">here</a>. This data file contains the coordinates of the atoms, their masses, and their charges, and can be directly read by LAMMPS using the read_file command. Let us replicate these atoms using LAMMPS, and apply an annealing procedure to obtain a block of amorphous silica.
<br><br>
Create a new input file in the same folder as the downloaded dataSiO.data, and copy the following lines in it.
<pre><code># Initialization
units metal
boundary p p p
atom_style full
pair_style vashishta
neighbor 1.0 bin
neigh_modify delay 1
# System definition
read_data data.SiO
replicate 4 4 4
</code></pre></p><p>
Metalic units are used as required by the Vashishta potential. The <a href="https://pubmed.ncbi.nlm.nih.gov/9993674/" target="_blank" class="removelinkdefault">Vashishta</a> potential is a bond-angle energy based potential, it deduces the bonds between atoms from their relative positions. Therefore, there is no need to provide bond and angle information as we do with classic force field. These potentials are more computationally heavy than classical force fields and require the use of a smaller timestep, but they allow for the modelling of bond formation and breaking, which is what we need here as we want to create a crack in the silica. You can download it by clicking <a href="Data/tutorial05/SiO.1990.vashishta" target="_blank" class="removelinkdefault">here</a>. The system is then replicated four times in all three directions of space.
<br><br>
Then, let us add the pair coefficient: we simply indicate that the first atom is Si, and the second is O, and add a dump for printing atoms trajectories:
<pre><code># Simulation settings
pair_coeff * * SiO.1990.vashishta Si O
dump dmp all atom 5000 dump.lammpstrj
</code></pre></p><p>
Finally, let us create the last part of our script. The annealing procedure is the following: we first start with a small phase at 6000 K, then cool down the system to 4000 K using a pressure of 100 atm. Then we cool down the system further while also reducing the pressure, then perform a small equilibration step at the final desired condition, 300 K and 1 atm.
<p class="alert alert-info"><b>Disclaimer.</b> I created this procedure by intuition and not from proper calibration, do not copy it without making your own tests if you intend to publish your results.</p>
<pre><code># Run
velocity all create 6000 4928459 rot yes dist gaussian
fix nvt1 all nvt temp 6000 6000 0.1
timestep 0.001
thermo 1000
run 5000
unfix nvt1
fix npt1 all npt temp 6000 4000 0.1 aniso 100 100 1
run 50000
fix npt1 all npt temp 4000 300 0.1 aniso 100 1 1
run 200000
fix npt1 all npt temp 300 300 0.1 aniso 1 1 1
run 4000
write_data data.amorphousSiO
</code></pre></p><p>
Notice the use of an anisotropic barostat for which all three directions of space are managed independently. An anisotropic barostat is better for a solid phase like here, but for a liquid of a gas, use an isotropic one instead. After running the simulation, the final configuration data.amorphousSiO will be located in the same folder as your input file. Alternatively, if you are only interested in the next steps of this tutorial, you can download it by clicking <a href="Data/tutorial05/data.amorphousSiO " target="_blank" class="removelinkdefault">here</a>. The final system resembles the image below, where the oxygen atoms are in red and the silicon atoms in yellow:
<p><img src="Figures/tutorial05/silica.jpg" class="img-fluid" alt="Responsive image"></p>
<p class="alert alert-info">Ideally you want to test the validity of the generated structure, for example by measuring the radial distribution function and/or the Young modulus, and compare them to experimental data. This is beyond the scope of this tutorial.</p>
</section>
<br><br><br><section id="Cracking">
<h2>Cracking the silica</h2>
We are now going to dilate the block of silica to create a crack. Create a new folder, copy data.amorphousSiO and SiO.1990.vashishta in it, and create a new input.lammps file starting with:
<pre><code># Initialization
units metal
boundary p p p
atom_style full
neighbor 1.0 bin
neigh_modify delay 1
# System definition
read_data data.amorphousSiO
# Simulation settings
pair_style vashishta
pair_coeff * * SiO.1990.vashishta Si O
dump dmp all atom 1000 dump.lammpstrj
</code></pre></p><p>
Then, we are going to progressively increase the size of the box over z, thus forcing the silica to crack. To do so, we are going to make a loop using the jump command. At every step of the loop, the box dimension over z will be multiplied by a factor 1.005. For this step, we a NVT thermostat because we want to impose a deformation of the volume (i.e. NPT would be inappropriate). Add the following lines to the input script:
<pre><code># Run
fix nvt1 all nvt temp 300 300 0.1
timestep 0.001
thermo 1000
variable var loop 35
label loop
change_box all z scale 1.005 remap
run 2000
next var
jump input.lammps loop
run 20000
write_data data.dilatedSiO
</code></pre></p><p>
After the dilatation, a final equilibration step of 20 picoseconds is performed. If you look at the dump file produced after executing this script, or at <a href="https://www.youtube.com/watch?v=8rBqYIcTgno&ab_channel=SimonGravelle" target="_blank" class="removelinkdefault">this video</a>, you can see the dilatation occurring step-by-step and the atoms adjusting to the box size. At first, the deformations are reversible (elastic regime), but at some point, bonds start breaking and dislocations appear (plastic regime). You can download the final state directly by clicking <a href="Data/tutorial05/data.dilatedSiO " target="_blank" class="removelinkdefault">here</a>. The final system, with the crack, resembles:
<p><img src="Figures/tutorial05/crackedsilica.jpg" class="img-fluid" alt="Responsive image"></p>
<p class="alert alert-info">In ambient conditions, silicon (Si) atoms are chemically passivated by forming covalent bonds with hydrogen (H) atoms. For the sake of simplicity, we are not going to decorate our silica with hydrogen atoms in this tutorial. It would require a procedure to properly insert hydrogen atoms at the right place. Also, it would require the use of another potential.</p>
</section>
<br><br><br><section id="Water">
<h2>Adding water</h2>
In order to add the water molecules to the silica, we are going to use the Monte Carlo method in the grand canonical ensemble (GCMC). In short, the system is put into contact with a virtual reservoir of given chemical potential \(\mu\), and multiple attempts to insert water molecules at random positions are made. Attempts are either accepted or rejected based on a Monte Carlo an acceptance rule.
<br><br>
In a different folder, copy the data.dilatedSiO file previously generated as well as this <a href="Data/tutorial05/TIP4P2005.txt" target="_blank" class="removelinkdefault">TIP4P2005.txt</a> file for the water molecule, and the SiO.1990.vashishta file. Create a new input file, and copy the following lines into it:
<pre><code># Initialization
units metal
boundary p p p
atom_style full
neighbor 1.0 bin
neigh_modify delay 1
pair_style hybrid/overlay vashishta lj/cut/tip4p/long 3 4 1 1 0.1546 10
kspace_style pppm/tip4p 1.0e-4
bond_style harmonic
angle_style harmonic
</code></pre></p><p>
There are several difference with the previous input files used in this tutorial. All these differences are here to include the water to our simulation. First, we have to combine two force fields, vashishta for SiO, and lj/cut/tip4p/long for TIP4P water model. Combining force fields in LAMMPS can be done using the hybrid/overlay pair style. We also need a kspace solver to solve the long range Coulomb interaction associated with tip4p/long. Finally, we need to define the style for the bond and angle of the water molecules. Some of these features have been seen in a <a href="tutorial03.html">previous and simpler tutorial</a>. Before going further, we need to make a few change to our data file. Currently, data.dilatedSiO only includes two type of atoms, but we need four:
<pre><code> LAMMPS data file via write_data, version 30 Jul 2021, timestep = 90000
576 atoms
2 atom types
0.910777522101565 19.67480018949893 xlo xhi
2.1092682236518137 18.476309487947546 ylo yhi
-4.1701120819606885 24.75568979356097 zlo zhi
Masses
1 28.0855
2 15.9994
Atoms # full
(...)
</code></pre></p><p>
We need to make some changes for the addition of water molecules. Modify the file so that it looks like that (the changes are highlighted in bold):
<pre><code>LAMMPS data file via write_data, version 30 Jul 2021, timestep = 90000
576 atoms
<b>4</b> atom types
<b>1 bond types</b>
<b>1 angle types</b>
<b>2 extra bond per atom
1 extra angle per atom
2 extra special per atom</b>
0.910777522101565 19.67480018949893 xlo xhi
2.1092682236518137 18.476309487947546 ylo yhi
-4.1701120819606885 24.75568979356097 zlo zhi
Masses
1 28.0855
2 15.9994
<b>3 15.9994
4 1.008</b>
Atoms # full
(...)
</code></pre></p><p>
Doing so, we anticipate that there will be 4 atoms types in the simulations, with O and H of H2O being indexes 3 and 4, respectively. There will also be 1 bond type and 1 angle type. The extra bond, extra angle, and extra special lines are for memory allocation. Now we can continue to fill the input.lammps file, by adding the system definition:
<pre><code># System definition
read_data data.dilatedSiO
molecule h2omol TIP4P2005.txt
create_atoms 0 single 19.5 10 10 mol h2omol 45585
group SiO type 1 2
group H2O type 3 4
</code></pre></p><p>
After reading the data file and defining the h2omol molecule from the txt file, the create atoms command is used to include one single molecule in the system at the location 19.5 10 10. I've chosen these 3 values so that the water molecule is initially inside the crack, and not overlapping with SiO atoms, you may have to choose different values. You can estimate these values by looking at your silica structure using VMD. There are cleaner way to do it, such as detecting the free space before inserting the molecule, but for the sake of simplicity, here we are going to do the dirty way here. Not adding a molecule before starting the GCMC steps usually lead to failure, I am not sure why. Then, add the following settings of the simulation:
<pre><code># Simulation settings
pair_coeff * * vashishta SiO.1990.vashishta Si O NULL NULL
pair_coeff * * lj/cut/tip4p/long 0 0
pair_coeff 1 3 lj/cut/tip4p/long 0.0057 4.42 # epsilonSi = 0.00403, sigmaSi = 3.69
pair_coeff 2 3 lj/cut/tip4p/long 0.0043 3.12 # epsilonO = 0.0023, sigmaO = 3.091
pair_coeff 3 3 lj/cut/tip4p/long 0.008 3.1589
pair_coeff 4 4 lj/cut/tip4p/long 0.0 0.0
bond_coeff 1 0 0.9572
angle_coeff 1 0 104.52
variable oxygen atom "type==3"
group oxygen dynamic all var oxygen
variable nO equal count(oxygen)
fix myat1 all ave/time 100 10 1000 v_nO file numbermolecule.dat
fix shak H2O shake 1.0e-4 200 0 b 1 a 1 mol h2omol
dump dmp all atom 2000 dump.lammpstrj
</code></pre></p><p>
The force field vashishta apply only to Si and O of SiO, and not to the O and H of H2O thanks to the NULL parameters. Pair coefficients for lj/cut/tip4p/long are defined between O atoms, as well as between O(SiO)-O(H2O) and Si(SiO)-O(H2O). Finally, the number of oxygen atoms will be printed in the file numbermolecule.dat, and the shake algorithm is used to maintain the shape of the water molecule over time. Some of these features have been seen previously, such as in <a href="tutorial03.html">this tutorial</a>.
Let us make a first equilibration step:
<pre><code># Run
compute_modify thermo_temp dynamic yes
compute ctH2O H2O temp
compute_modify ctH2O dynamic yes
fix mynvt1 H2O nvt temp 300 300 0.1
fix_modify mynvt1 temp ctH2O
compute ctSiO SiO temp
fix mynvt2 SiO nvt temp 300 300 0.1
fix_modify mynvt2 temp ctSiO
timestep 0.001
thermo 1000
run 5000
</code></pre></p><p>
We use to different thermostats for SiO and H2O, which is better when you have two species, such as one solid and one liquid. It is particularly important to use two thermostats here as the number of water molecules will fluctuate. We use a compute_modify 'dynamic yes' for water to specify that the number of molecules is not constant.
Finally, let us use the gcmc fix and perform the grand canonical Monte Carlo step:
<pre><code>variable tfac equal 5.0/3.0
fix fgcmc H2O gcmc 100 100 0 0 65899 300 -0.5 0.1 mol h2omol tfac_insert ${tfac} group H2O shake shak full_energy
run 250000
write_data SiOwithwater.data
</code></pre></p><p>
The tfac_insert option ensures that the correct estimate is made for the temperature of the inserted water molecules by taking into account the internal degrees of freedom. Running this simulation, you should see the number of molecule increasing progressively. Using \(\mu = -0.5\) eV is a reasonable value for the chemical potential that corresponds roughly to ambient conditions (i.e. RH approx 50%).
<br><br>
After 250000 steps, I've got 12 water molecules. The number will vary from one simulation to another. Here the final number of molecule will depend a lot on the space available in the crack. It also depends on the imposed chemical potential, temperature, and interaction between water and silica. The final state looks like this:
<p><img src="Figures/tutorial05/SiOwater.jpg" class="img-fluid" alt="Responsive image"></p>
<p>You can download the input scripts that have been written in this tutorial by clicking <a href="Inputs/tutorial05.zip" target="_blank" class="removelinkdefault">here</a>.</p>
</section>
<br><br><br><section id="further">
<h2>Going further</h2>
<p>This tutorial is already quite complicated, but if you want to go further, there are several interesting things that can be done with this system:
<br><br>
<b>Relative humidity.</b> You can link the imposed chemical potential with the value of relative humidity (RH). For that, you have to calibrate your simulation by measuring the equilibrium amount of water in an empty box for varying imposed chemical potential, see <a href="https://aip.scitation.org/doi/full/10.1063/1.5126481" target="_blank" class="removelinkdefault">one of my paper</a> for example.
<br><br>
<b>Isotherm.</b> You can perform a full adsorption isotherm by varying the chemical potential and extracting the equilibrium water content as a function of the imposed RH. Isotherms can be compared to experimental data, and are used sometimes to calibrate force field as it contains a lot of information about the fluid-solid interactions.
<br><br>
<b>Isosteric heat of adsorption.</b> Using the GCMC procedure, you can measure the heat of adsorption by measuring the fluctuations in water molecule and total energy of the system.
<br><br>
<b>Hybrid MD/GCMC.</b> In the grand canonical ensemble, the volume of the box is fixed, so its not possible to capture the swelling of a material with its water content (most material swells with water, like sponges). If you want to model the swelling while also performing a GCMC adsorption simulation, you can alternate between GCMC steps and molecular dynamics steps in the NPT ensemble. This method is called hybrid MD/GCMC.
</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>
</div>
</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>