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
|
<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01//EN" "http://www.w3.org/TR/html4/strict.dtd">
<html><head>
<meta content="text/html; charset=ISO-8859-1" http-equiv="content-type"><title>Tutorial paral Molecular Dynamics</title></head><body><br>
<hr><h1>ABINIT,
lesson Molecular Dynamics :</h1>
<h2>How to perform Molecular Dynamics calculations using parallelism </h2>
<hr><p>This lesson aims at showing how to perform molecular dynamics with ABINIT using a parallel computer. </p>
You will learn how to launch molecular dynamics calculation and what are
the main
input
variables that govern convergency and numerical efficiency.<br>
<p>You are supposed to know already some basics of parallelism in ABINIT, explained in the tutorial
<a href="lesson_basepar.html">A first introduction to ABINIT in parallel</a>, and
<a href="lesson_paral_gspw.html">ground state with plane waves</a>.
<p>This lesson should take about 1.5 hour to be done and requires
to have at least a 200 CPU core parallel computer.<i> </i> </p>
<h5>Copyright (C) 2000-2014 ABINIT group (JB) <br>
This file is distributed under the terms of the GNU General Public
License, see ~ABINIT/COPYING or <a href="http://www.gnu.org/copyleft/gpl.txt" target="_blank">
http://www.gnu.org/copyleft/<wbr>gpl.txt </a>. <br>
For the initials of contributors, see ~ABINIT/Infos/contributors . </h5>
<script type="text/javascript" src="list_internal_links.js"> </script>
<h3><b><b>Contents of lesson PARAL_MOLDYN :</b></b></h3>
<ul><li><a href="lesson_paw1.html#0" target="_blank">0.</a> Summary of the molecular dynamics method</li><li><a href="lesson_paw1.html#1" target="_blank">1.</a> Performing molecular dynamics with ABINIT</li><li><a href="lesson_paw1.html#2" target="_blank">2.</a> The convergence on K-points and supercell size</li><li><a href="lesson_paw1.html#3" target="_blank">3.</a> Compute the melting temperature of Al at a given pressure</li></ul>
<hr>
<h3><a name="1334abc7eee58317_0"></a><b><br>0. Summary of </b>the molecular dynamics<b> method</b></h3>
<p style="text-align: justify;">The basic idea underlying <span style="font-style: italic;">Ab Initio Molecular Dynamics</span> (AIMD) is to
compute
the forces acting on the nuclei from electronic structure calculations
that are performed as the molecular dynamics trajectory is
generated. An AIMD calculation assumes only that the system is
composed of <em></em>nuclei and<em></em>
electrons, that the Born–Oppenheimer approximation is valid, and that
the dynamics of the nuclei can be treated classically on the
ground-state electronic surface. It
allows both equilibrium thermodynamic and dynamical properties of a
system at finite temperature to be computed. For exemple melting
temperatures, phase transitions, atomic vibrations, structure factor...
but also XANES or IR spectrum can be obtained with this technique. AIMD
deals with supercells of hundred to thousand of atoms (usually, the
larger, the better!). In addition Molecular Dynamics simulations can be
performed for days, weeks or even months! They are therefore very time
consuming and can not be done without the help of high speed and
massively parallel computing.</p>
<p><br>
</p>
<hr style="width: 100%; height: 2px;">
<p><i>Before continuing, you might
consider to work in a different
subdirectory as for the other lessons. Why not "Work_paral_moldyn" ?
</i><i>In what follows, the name of files are
mentioned as if
you were in this subdirectory.<br>
All the input files can be found in the </i><i><span style="color: rgb(0, 102, 0);">~abinit/tests/tutorial/Input</span>
directory.</i></p>
<p style="font-style: italic;">In the following, when <span style="font-weight: bold;">"</span><span style="color: rgb(0, 153, 0);">run ABINIT
over nn CPU cores</span>" appears, you have
to use a specific command line according to the operating system
and
architecture of the computer you are using. This can be for
instance:
<span style="color: rgb(0, 102, 0);">mpirun -n nn abinit < abinit.files </span>or the use of a specific
<span style="color: rgb(0, 153, 0);">submission file</span>. </p>
<p><i>You can compare your results with several reference
output files located in </i><i><span style="color: rgb(0, 102, 0);">~abinit/tests/tutorial/Refs</span>
and </i><i><span style="color: rgb(0, 102, 0);">~abinit/tests/tutorial/Refs/<wbr>
</span></i><i>directories (for the present tutorial they are named tmoldyn_*.out). <br>
</i></p>
<h3><hr> </h3>
<h3><a name="1334abc7eee58317_1"></a> <b><b><br>
1. </b></b>Performing molecular dynamics<b><b> with ABINIT</b></b></h3>
<br>
<p style="text-align: justify;">There are different algorithms to do molecular dynamics.
See the input variable "<a href="../input_variables/varrlx.html#ionmov" target="kwimg" onclick="return (false);">ionmov</a>",
with values 1, 6, 7, 8, 9, 12, 13 and 14.<span style="text-decoration: underline;"></span><a href="../input_variables/varrlx.html#dtion" onclick="return (false);"><br>
</a></p>
<a href="../input_variables/varrlx.html#dtion" onclick="return (false);">dtion</a> controls the ion time steps in
atomic units of time (one atomic time unit is
2.418884e-17 seconds, which is the value of
Planck's constant in hartree*sec). The default value is 100. You should try several values
for <b>dtion</b> in order to establish the stable and efficient
choice.
For example this value should decrease at high pressure.<br>
<br>
Except
for the <span style="font-style: italic;">isothermal</span>/<span style="font-style: italic;">isenthalpic</span> (<a href="../input_variables/varrlx.html#ionmov" onclick="return (false);">ionmov</a> 13) ensemble the input variable "<a href="../input_variables/varrlx.html#optcell">optcell</a>"
must be set to 0. <br>
You have also to define the maximal number of
timesteps of the molecular dynamics.<br>
<br>
Usually you can set the input
variable "<a href="../input_variables/varrlx.html#ntime" target="kwimg">ntime</a>"
to a large value, 5000, since there is no "end" to a molecular dynamics
simulation. You can always stop or restart the calculation at your
convenience by using the input variable <a href="../input_variables/varrlx.html#restartxf" onclick="return (false);">restartxf</a>.
<p style="text-align: justify;">The input file<span style="color: red;"> tmoldyn_01.in</span>
is an example of a file
that contains data for a molecular dynamics simulation using the isokinetic ensemble for aluminum.
Open the <span style="color: red;">tmoldyn_01.in</span><span style="color: rgb(204, 0, 0);"></span>
file and look at it carefully. The
unit cell is defined at the end. It is a <span style="font-weight: bold;">2x2x2 fcc supercell</span> containing 32 atoms of Al. <a href="../input_variables/varrlx.html#ionmov" target="kwimg" onclick="return (false);">ionmov</a> is set to 12 for the isokinetic ensemble, and since <a href="../input_variables/varrlx.html#ntime" target="kwimg">ntime</a> is set to 50, ABINIT
will carry on 50 time steps of molecular dynamics. The calculation will
be performed for a temperature of 3000 K, see the key variable <a href="../input_variables/varrlx.html#mdtemp" onclick="return (false);">mdtemp</a>.<span style="font-weight: bold;"> </span>It
gives the initial and final temperature in Kelvin of the
simulation. The temperature will change linearly from the initial
temperature
<b>mdtemp(1)</b> at itime=1 to
the final temperature <b>mdtemp(2)</b> at the end of the
<a href="../input_variables/varrlx.html#ntime" onclick="return (false);">ntime</a> timesteps.
Here the temperature will stay constant during the whole simulation.<br>
Note that we use the same temperature for the ions and the electrons : <a href="../input_variables/varbas.html#occopt" onclick="return (false);">occopt</a> has been set to 3 for a Fermi-Dirac smearing and <a href="../input_variables/vargs.html#tsmear" onclick="return (false);">tsmear</a><a href="../input_variables/vargs.html#tsmear" onclick="return (false);">r</a>
has been set to 3000 Kelvin. Nothing prevent you to use different
electronic and ionic temperature, you just have to know why you are
doing so!</p>
<p style="text-align: justify;">Molecular
dynamics simulations are always large calculations, dealing with
supercells of hundreds to thousands of atoms. Therefore they are always
performed in parallel. In <span style="color: red;">tmoldyn_01.in<span style="color: black;">,</span> <span style="color: black;"><a href="../input_variables/varpar.html#paral_kgb" onclick="return (false);">paral_kgb</a> has been set to 1 </span></span>to activate the parallelisation over K-points, G-vectors and bands. The three following keywords give
the number of processors for each level of parallelisation. Since we
have only one K-point in the simulation (the gamma point), <a href="../input_variables/varpar.html#npkpt" onclick="return (false);">nkpt</a> has been set to 1. <a href="../input_variables/varpar.html#npfft" onclick="return (false);">npfft</a> is set to 3 and <a href="../input_variables/varpar.html#npband" onclick="return (false);">npband</a> to 10, for a total number of 3x10=30 processors. You might use the <span style="color: rgb(204, 0, 0);">tmoldyn.files</span>
file.
Edit it and adapt it with the appropriate file names.</p>
<p style="text-align: justify;"><span style="color: black;">Then</span> run the
calculation in parallel over <span style="font-weight: bold;">30
CPU cores</span>.
You can change the distritbution of processors over the level of
parallelisation to try to find the most efficient one. Set for
exemple <a href="../input_variables/varpar.html#npfft" onclick="return (false);">npfft</a> to 1 and <a href="../input_variables/varpar.html#npband" onclick="return (false);">npband</a>
to 40. You can make other choices and compare the individual cpu time.
Since
molecular dynamics can last for weeks, it is crucial to find the
appropriate distribution to reduce the computational time at the
maximum. Look at the ouput file. For each iteration you will see the
coordinates, the forces, the velocities and the kinetic and the total
energy.</p>
<p style="text-align: justify;">In addition, ABINIT should have
generated a <span style="font-weight: bold;">HIST file</span>, which contains the whole history of the
molecular dynamics simulation : atomic positions, velocities, primitive
translations,
stress tensor, energies... at each time step. This file will be used to
restart the calculation if you want to perform more time steps or to
extract the necessary informations to make use of the molecular
dynamics simulation. In <span style="color: red;">tmoldyn_01.in </span>add the keyword <a href="../input_variables/varrlx.html#restartxf" onclick="return (false);">restartxf</a>
and set it to -1. Run the calculation again, in the same directory.
Look at the new output file. The number of each time step are indicated
over the total number of steps :</p>
<p><span style="color: rgb(0, 102, 0);">--- Iteration: ( 1/100) Internal Cycle: (1/1)</span></p><p style="text-align: justify;">Since
we already performed 50 steps of molecular dynamics, the total number
of time steps are now 100. So the first 50 iterations are from the
previous calculation. You can check that by comparing <span style="color: red;">tmoldyn_01.out <span style="color: black;">and</span> </span><span style="color: red;">tmoldyn_01.outA<span style="color: black;">.
There is only one <span style="font-weight: bold;">HIST file</span> and it contains the history of the two
calculations.</span></span></p>
<p style="text-align: justify;"><span style="color: red;"><span style="color: black;">Now we can calculate and plot several quantities. We need
for that the <span style="font-style: italic; font-weight: bold;">diag_moldyn.py</span> python script. You can find it in the </span></span><i><span style="color: rgb(0, 102, 0);">~abinit/doc/tutorial/lesson_paral_moldyn</span></i><span style="color: red;"><span style="color: black;"> directory (link <a href="./lesson_paral_moldyn/diag_moldyn.py">here</a>).<br>
Run the diag_modyn.py script (type: "python diag_moldyn.py").<br>
</span></span></p>
<p style="text-align: justify;"><span style="color: red;"><span style="color: black;">You can read (on standardoutput) the average value and the standard
deviation of the total energy, the temperature and the pressure. You have also generated
several files which contain pressures, energies, stresses, positions
and temperatures. You can plot this files to observe the behavior of
the quantities during the molecular dynamics. Note that 100 time steps
is far from being sufficient to equilibrate physical quantities as
pressure. 2000 or 3000 are more common numbers to reach this goal but it
would exceed the time alloted to this tutorial.</span></span></p>
<h3><hr> <a name="1334abc7eee58317_2"></a> <b><b><br>
2. The convergence on K-points and supercell size</b></b></h3>In the previous section you have learned to perform molecular dynamics
with abinit. We used a fcc supercell of 32 atoms with only one K-point.
In this section we will make convergence studies
with respect to these parameters.
<p><b><br>
</b></p>
<p><b>1.a</b><b> Computing the convergence in K-points</b></p>
<p style="text-align: justify;">The files <span style="color: red;">tmoldyn_02.in <span style="color: black;">and</span> </span><span style="color: red;">tmoldyn_03.in <span style="color: black;">are input files for 2x2x2 and 3x3x3 K-points grid respectively, or 4 and 14 K-points in the </span></span>irreducible
Brillouin zone<span style="color: red;"><span style="color: black;">. Since the parallelisation is the most efficient </span></span>over the k-point level you should always put <a href="../input_variables/varpar.html#npkpt" onclick="return (false);">nkpt</a> to the largest possible value before increasing <a href="../input_variables/varpar.html#npfft" onclick="return (false);">npfft</a> and <a href="../input_variables/varpar.html#npband" onclick="return (false);">npband</a><span style="color: black;">.
We have followed this rule in the input files. Change the name of the
previous file PRESS to PRESS01 to save it. Run now ABINIT </span>in parallel over <span style="font-weight: bold;">120
CPU cores</span> with <span style="color: red;">tmoldyn_02.in </span><span style="color: black;"> and</span> over <span style="font-weight: bold;">140
CPU cores</span> with <span style="color: red;">tmoldyn_03.in<span style="color: black;">. At the end of each calculation use</span></span><span style="color: red;"><span style="color: black;"> the <span style="font-style: italic;">diag_moldyn.py</span></span></span><span style="color: red;"><span style="color: black;"> script and save the results in PRESS02 and PRESS03. You can now plot
the pressures in term of the K-points grids and compare the average
values:</span></span></p>
<p style="text-align: center;"><span style="color: red;"><span style="color: black;"><img style="width: 443px; height: 342px;" alt="Kpoints" src="./lesson_paral_moldyn/Kpoints.png"><br>
</span></span></p>
<p style="text-align: justify;"><span style="color: red;"><span style="color: black;">
As said previously our simulations are too short to be completely
convincing but you can see that you need at least a 2x2x2 K-points grid
for a 32 atoms cell. If you have some time, increase ntime to 300 and run again ABINIT.</span></span></p>
<p><b><br>
</b></p>
<p><b>1.b</b><b> Computing the convergence in cell size</b></p>
<p style="text-align: justify;">We
also have to check if our cell is sufficiently large to give reliable
physical quantities. In the previous section we used a 2x2x2 <span style="font-style: italic;">fcc</span>
supercell. <span style="color: red;">tmoldyn_04.in <span style="color: black;">is an input file for a 3x3x3 <span style="font-style: italic;">fcc</span> supercell and therefore contains 108 atoms. <a href="../input_variables/varbas.html#nband" onclick="return (false);">nband</a> and <a href="../input_variables/varbas.html#acell" onclick="return (false);">acell</a> has been scaled accordingly to take into account the new size of the cell. </span></span><span style="color: black;">Run now ABINIT </span>in parallel over <span style="font-weight: bold;">45
CPU cores</span> and then <span style="color: red;"><span style="color: black;"><span style="font-style: italic;">diag_moldyn.py</span></span></span>
(note that the output file is very big, and no reference has been provided for comparison).
Save the pressure to PRESS04. <span style="color: red;">tmoldyn_05.in <span style="color: black;">has
the same cell but a 2x2x2 K-points grid (note that the output file is very big, and no reference has been provided for comparison).
Run it over the adequate
number of cores and save the pressure to PRESS05. Plot now PRESS04 and
PRESS05 and compare the average values. You'll see that for this size
of cell, the gamma point is sufficient.</span></span></p>
<p style="text-align: justify;"><span style="color: red;"><span style="color: black;">We are now going to increase
again the cell size. With a 4x4x4 fcc cell, the file </span></span><span style="color: red;">tmoldyn_06.in <span style="color: black;">has 256 atoms. Of course, </span></span><span style="color: red;"><span style="color: black;"><a href="../input_variables/varbas.html#nband" onclick="return (false);">nband</a> and <a href="../input_variables/varbas.html#acell" onclick="return (false);">acell</a>
has been scaled. This calculation should last for 30 min over <span style="font-weight: bold;">60 CPU cores</span>
(note that the output file is very big, and no reference has been provided for comparison). Run
it and save the pressure to PRESS06. Plot now PRESS02, PRESS04 and
PRESS06, remove the first steps and </span></span><span style="color: red;"><span style="color: black;">compare
the pressure average values:</span></span></p>
<div style="text-align: center;"><span style="color: red;"><span style="color: black;"><img style="width: 545px; height: 422px;" alt="cells" src="lesson_paral_moldyn/cell.png"></span></span><br>
</div>
<span style="color: red;"><span style="color: black;"></span></span><span style="color: red;"><span style="color: black;"><br>
</span></span><span style="color: red;"><span style="color: black;"><br>
</span></span>
<div style="text-align: justify;"><span style="color: red;"><span style="color: black;">You can see that even if the pressure was
converge in term of K-points, the 32 atoms supercell was not sufficient to
give a reliable pressure. A 3x3x3 supercell with the gamma point seems to be the adequate
size. You see also that with a bigger cell, the pressure fluctuations are considerably reduced.</span></span><br>
</div>
<span style="color: red;"><span style="color: black;"><span style="font-style: italic;"><br>
Note that here we made the convergence studies using the</span> <span style="font-style: italic;">pressure as a criteria. The results can depend on the physical quantity you are looking at, pressure, temperature, energy, or</span></span></span><span style="font-style: italic;">
dynamical matrices by observing the displacement fluctuations....
Always check if your cell is large enough and give the corresponding
uncertainty.<br>
Also, to reduce the time necessary to do this tutorial we set the value
of ecut to 3Ha. This is too small, for Al, it should be closer to 8 Ha.<br>
<br>
</span><hr>
<h3><a name="1334abc7eee58317_3"> </a><b><b><br>
3. Compute the melting temperature of Al at a given pressure
</b></b></h3>
<div style="text-align: justify;">As an exemple of what can be done in molecular dynamics,
we are going to calculate the melting temperature of aluminum using the
so-called <span style="font-style: italic;">Heat Until it Melts</span> (HUM) method. In this method the solid
phase is heated gradually until melting occurs. Let us start with a
temperature of 5500K.<br>
An example of file is given with <span style="color: red;">tmoldyn_07.in<span style="color: black;">. To work fast, we use a 32 atoms supercell and the gamma point (note that the output file is very big, and no reference has been provided for comparison).<br>
<br>
</span></span><span style="color: black;">Run ABINIT </span>in parallel over <span style="font-weight: bold;">30
CPU cores</span> and then <span style="color: red;"><span style="color: black;"><span style="font-style: italic;">diag_moldyn.py</span></span></span>.<br>
Save the pressure to PRESS71. Plot the atomic positions, you see
that at this temperature, the cell is solid. Increase the temperature
to 6000K (do not forget to also change the electronic temperature
with tsmear) and run ABINIT. Save the pressure and again, look at the
positions. The cell is still solid.<br>
Set the temperature to 6500K. Look
at the positions: the ions are moving across the cell and do not come
back to their equilibrium positions. The cell has melted and is now
liquid. Plot the pressures for the three simulations: <br>
</div>
<br>
<div style="text-align: center;"><img style="width: 453px; height: 349px;" alt="melting" src="./lesson_paral_moldyn/melting.png"><br>
</div>
<div style="text-align: center;"><img style="width: 361px; height: 278px;" alt="xcart" src="./lesson_paral_moldyn/xcart.png"><br>
</div>
<br>
<div style="text-align: justify;">You can clearly
observe a discontinuous change in pressure due to the volume
difference between the solid and liquid phases. This give a melting
temperature of 6250 K at 149 GPa. With more sophisticated techniques
the melting temperature at this pressure is around 5500 K. Indeed,
in addition to the crude parameters we used (ecut, natom...), the HUM method has some intrinsec
drawbacks. In HUM the crystal is heated homogeneously, the melting
initiates in the bulk and this results in an overheating effect.<br>
<br>
<br>
</div>
<script type="text/javascript" src="list_internal_links.js"> </script>
</body></html>
|
