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
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
|
<html>
<head>
<title>Tutorial #3</title>
</head>
<body bgcolor="#ffffff">
<hr>
<h1>ABINIT, third lesson of the tutorial: </h1>
<h2>Crystalline silicon. </h2>
<hr>
<p>This lesson aims at showing you how to get
the following physical properties, for an insulator :
<ul>
<li>the total energy
<li>the lattice parameter
<li>the band structure (actually, the Kohn-Sham band structure)
</ul>
You will learn about the use of k-points, as well
as the smearing of the plane-wave kinetic energy cut-off.
<p>This lesson should take about 1 hour.
<h5>Copyright (C) 2000-2014 ABINIT group (XG,RC)
<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">
http://www.gnu.org/copyleft/gpl.txt </a>.
<br> For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
</h5>
<script type="text/javascript" src="list_internal_links.js"> </script>
<br>
<h3><b>Content of lesson 3</b></h3>
<ul>
<li><a href="lesson_base3.html#31">3.1</a> Computing the total energy of silicon at
fixed number of k points.</li>
<li><a href="lesson_base3.html#32">3.2</a> Starting the convergence study with respect
to k points</li>
<li><a href="lesson_base3.html#33">3.3</a> Actually performing the convergence study
with respect to k points</li>
<li><a href="lesson_base3.html#34">3.4</a> Determination of the lattice parameters</li>
<li><a href="lesson_base3.html#35">3.5</a> Computing the band structure</li>
</ul>
<br>
<hr>
<a name="31"> </a>
<h3><b>
3.1.
Computing the total energy of silicon at fixed number of k points.
</b></h3>
<p><i>Before beginning, you might consider to work in a different subdirectory
as for lesson_base1 or lesson_base2 . Why not "Work3" ?</i>
<p>The file ~abinit/tests/tutorial/Input/tbase3_x.files lists the file
names and root names. You can copy it in the Work3 directory and change it
as you did for the tbase1_x.files and tbase2_x.files files. You can also copy
the file ~abinit/tests/tutorial/Input/tbase3_1.in in Work3. This is your
input file. You should edit it, read it carefully, have a look at the following
"new" input variables, and their explanation :
<ul>
<li><a href="../input_variables/varbas.html#rprim" target="kwimg">rprim</a>
<li><a href="../input_variables/varbas.html#xred" target="kwimg">xred</a> (used instead of xcart)
<li><a href="../input_variables/varbas.html#kptopt" target="kwimg">kptopt</a>,
<a href="../input_variables/varbas.html#ngkpt" target="kwimg">ngkpt</a>,
<a href="../input_variables/varbas.html#nshiftk" target="kwimg">nshiftk</a>,
<a href="../input_variables/varbas.html#shiftk" target="kwimg">shiftk</a>,
<a href="../input_variables/vargs.html#kptrlatt" target="kwimg">kptrlatt</a>,
(not easy ... take your time !)
<li><a href="../input_variables/vargs.html#diemac" target="kwimg">diemac</a> (compared to isolated
molecules, another value is used, while <a href="../input_variables/vargs.html#diemix" target="kwimg">diemix</a>
has been suppressed).
</ul>
<p>Note also the following :
you will work at fixed <a href="../input_variables/varbas.html#ecut" target="kwimg">ecut</a>
(=8Ha). It is implicit that in "real life",
you should do a convergence test with respect to
<a href="../input_variables/varbas.html#ecut" target="kwimg">ecut</a> ...
<br>
Here, a suitable <a href="../input_variables/varbas.html#ecut" target="kwimg">ecut</a> is given
to you. It will allow to obtain 0.2% relative accuracy on lattice parameters.
<p>When you have read the input file, you can run the code, as usual
(it will run for a few seconds).
<br>Then, read the output file, and note the total energy.
<pre>
etotal -8.8662238960E+00
</pre>
<hr>
<a name="32"> </a>
<h3><b>
3.2.
Starting the convergence study with respect to k points
</b></h3>
<p>There is of course a convergence study associated with the sampling
of the Brillouin zone. You should examine different grids,
of increasing resolution. You might try the following series
of grids :
<pre>
ngkpt1 2 2 2
ngkpt2 4 4 4
ngkpt3 6 6 6
ngkpt4 8 8 8
</pre>
However, the associated number of k points in the irreducible
Brillouin zone grows very fast. It is
<pre>
nkpt1 2
nkpt2 10
nkpt3 28
nkpt4 60
</pre>
ABINIT computes automatically this number of k points, from the definition of the grid
and the symmetries.
<br>You might nevertheless define an input nkpt value in the input file,
in which case ABINIT will compare its computed value (from the grid) with this input
value. We take this opportunity to examine the behaviour of ABINIT
when a problem is detected. Let's suppose that with <code>ngkpt1 4 4 4</code> ,
one mentions <code>nkpt1 2</code> .
The input file ~abinit/tests/tutorial/Input/tbase3_2.in is an example.
Do not forget to change tbase3_x.files,
if you are using that file name . The message that you get a few dozen of lines before the end of the
log file is :
<pre>
--- !BUG
message: |
The argument nkpt= 2, does not match
the number of k points generated by kptopt, kptrlatt, shiftk,
and the eventual symmetries, that is, nkpt= 10.
However, note that it might be due to the user,
if nkpt is explicitely defined in the input file.
In this case, please check your input file.
src_file: getkgrid.F90
src_line: 415
...
Action : contact ABINIT group.
</pre>
<p>This is a typical ABINIT error message. It states what is the problem that causes the stop of ABINIT,
then suggests that it might be due to an error in the input file, namely, an erroneous value of nkpt.
The expected value,
<a href="../input_variables/varbas.html#nkpt" target="kwimg">nkpt</a> 10.
is mentioned before the notice that the input file might be erroneous.
Then, the file at which the problem occured is mentioned, as well as the number of the line in that file. </p>
<p>As the computation of <a href="../input_variables/varbas.html#nkpt" target="kwimg">nkpt</a>
for specific grids of k points is not an easy task, while the even more
important selection of specific economical grids (the best ratio
between the accuracy of the integration in the Brillouin zone
and the number of k-points) is more difficult, some help to the
user is provided by ABINIT. ABINIT is able to examine automatically
different k point grids, and to propose the best grids for integration.
This is described in the abinit_help file, see the input
variable <a href="../input_variables/varfil.html#prtkpt" target="kwimg">prtkpt</a>,
and the associated characterisation of the integral accuracy, described
in <a href="../input_variables/vargs.html#kptrlen" target="kwimg">kptrlen</a>.
The generation of lists of k-point sets is done in different test cases,
in the directory ~abinit/tests/v2 . You can directly have a look at the output
files in ~abinit/testis/v2/Refs , the output files for the tests 61 to 73.
<p>When one begins the study of a new material, it is strongly advised
to examine first the list of k points grids, and select (at least)
three efficient ones, for the k point convergence study.
Do not forget that the CPU time will be linearly proportional to the
number of k points to be treated : using 10 k points will take five
more time than using 2 k points. Even for a similar accuracy of the
Brillouin zone integration (about the same value of
<a href="../input_variables/vargs.html#kptrlen" target="kwimg">kptrlen</a>),
it might be easy to generate a grid
that will fold to 10 k points in the irreducible Brillouin zone, as well as
one that will fold to 2 k points in the irreducible Brillouin zone.
The latter is clearly to be preferred !
<p>
<hr>
<p> <a name="33"> </a> </p>
<h3><b>
3.3.
Actually performing the convergence study with respect to k points
</b></h3>
<p>In order to understand k-point grids, you should read the Monkhorst and Pack
paper, <cite>Phys. Rev. B 13, 5188 (1976)</cite> ... Well, maybe not immediately
... In the meantime, you can try the above-mentioned convergence study.
<p>The input file ~abinit/tests/tutorial/Input/tbase3_3.in is an example,
while ~abinit/tests/tutorial/Refs/tbase3_3.out is a reference output
file. In this output file, you should have a look at the echo of input variables.
As you know, these are preprocessed, and, in particular, <a href="../input_variables/varbas.html#ngkpt" target="kwimg">ngkpt</a>
and <a href="../input_variables/varbas.html#shiftk" target="kwimg">shiftk</a> are used to generate
the list of k points (<a href="../input_variables/varbas.html#kpt" target="kwimg">kpt</a>) and
their weights (<a href="../input_variables/varbas.html#wtk" target="kwimg">wtk</a>). You should
read the information about <a href="../input_variables/varbas.html#kpt" target="kwimg">kpt</a>
and <a href="../input_variables/varbas.html#wtk" target="kwimg">wtk</a>.
<p>From the output file, here is the evolution of total energy per unit cell :
<pre>
etotal1 -8.8662238960E+00
etotal2 -8.8724909739E+00
etotal3 -8.8726017432E+00
etotal4 -8.8726056405E+00
</pre>
<p>The difference between dataset 3 and dataset 4 is rather small. Even the dataset
2 gives an accuracy of about 0.0001 Ha <br>
So, our converged value for the total energy, at fixed <a href="../input_variables/varbas.html#acell" target="kwimg">acell</a>,
fixed <a href="../input_variables/varbas.html#ecut" target="kwimg">ecut</a>, is -8.8726 Ha .
<p>
<hr>
<p> <a name="34"> </a>
<h3><b>
3.4.
Determination of the lattice parameters.
</b></h3>
<p>The input variable "<a href="../input_variables/varrlx.html#optcell" target="kwimg">optcell</a>"
governs the automatic optimisation of cell shape and volume. <br>
For the automatic optimisation of cell volume, use :
<pre>
optcell 1
ionmov 3
ntime 10
dilatmx 1.05
ecutsm 0.5
</pre>
<p>You should read the indications about <a href="../input_variables/varrlx.html#dilatmx" target="kwimg">dilatmx</a>
and <a href="../input_variables/varrlx.html#ecutsm" target="kwimg">ecutsm</a>. <br>
Do not test all the k point grids, only those with nkpt 2 and 10.
<p>The input file ~abinit/tests/tutorial/Input/tbase3_4.in is an example,
while ~abinit/tests/tutorial/Refs/tbase3_4.out is a reference output
file. <br>
<p>You should obtain the following evolution
of the lattice parameters :
<pre>
acell1 1.0233363682E+01 1.0233363682E+01 1.0233363682E+01 Bohr
acell2 1.0216447241E+01 1.0216447241E+01 1.0216447241E+01 Bohr
</pre>
with the following very small residual stresses :
<pre>
strten1 1.8591719160E-07 1.8591719160E-07 1.8591719160E-07
0.0000000000E+00 0.0000000000E+00 0.0000000000E+00
strten2 -2.8279720007E-08 -2.8279720007E-08 -2.8279720007E-08
0.0000000000E+00 0.0000000000E+00 0.0000000000E+00
</pre>
<p>The stress tensor is given in Hartree/Bohr^3, and the order
of the components is
<pre>
11 22 33
23 13 12
</pre>
<p>There is only a 0.13% relative difference between <a href="../input_variables/varbas.html#acell" target="kwimg">acell</a>1
and <a href="../input_variables/varbas.html#acell" target="kwimg">acell</a>2 . <br>
So, our converged LDA value for Silicon, with the 14si.pspnc pseudopotential
(see the tbase3_x.files file) is 10.216
Bohr (actually 10.21644...), that is 5.406 Angstrom. The experimental value is <code>5.431 Angstrom
at 25 degree Celsius</code>, see <cite>R.W.G. Wyckoff, Crystal structures
Ed. Wiley and sons, New-York (1963).</cite>
<p>
<hr>
<a name="35"> </a>
<h3><b>3.5.
Computing the band structure</b></h3>
<p>We fix the parameters <a href="../input_variables/varbas.html#acell" target="kwimg">acell</a>
to the theoretical value of 3*10.216, and we fix also the grid of k points (the
4x4x4 FCC grid, equivalent to a 8x8x8 Monkhorst-pack grid) <br>
We will ask for 8 bands (4 valence and 4 conduction).
<p>A band structure can be computed by solving the Kohn-Sham equation for
many different k points, along different lines of the Brillouin zone.
<br>The potential that enters the Kohn-Sham must be derived from a previous
self-consistent calculation, and will not vary during the scan of different
k-point lines.
<br>Suppose that you want to make a L-Gamma-X-(U-)Gamma circuit, with 10,
12 and 17 divisions for each line (each segment has a different length
in reciprocal space, and these divisions give approximately the same
distance between points along a line).
<br>The circuit will be obtained easily by the following choice
of segment end points :
<pre>
L (1/2 0 0)
Gamma (0 0 0)
X (0 1/2 1/2)
Gamma (1 1 1)
</pre>
<i>Note :
<ol>
<li>the last Gamma point is in another cell of the reciprocal
space than the first one, this choice allows to construct
the X-U-Gamma line easily ;</li>
<li>the k-points are specified
using reduced coordinates - in agreement with the input setting of the
primitive 2-atom unit cell - in standard textbooks, you will often find
the L, Gamma or X point given in coordinates of the conventional
8-atom cell : the above-mentioned circuit is then
(1/2 1/2 1/2)-(0 0 0)-(1 0 0)-(1 1 1), but such (conventional) coordinates
cannot be used with the 2-atom (non-conventional) cell.</li>
</ol>
</i>
<p>So, you should set up in your input file, for the first dataset, a usual SCF
calculation in which you output the density
(<a href="../input_variables/varfil.html#prtden" target="kwimg">prtden</a>
1), and, for the second dataset :
<ul>
<li>fix <a href="../input_variables/varbas.html#iscf" target="kwimg">iscf</a> to -2, to make a
non-self-consistent calculation ;
<li>define <a href="../input_variables/varfil.html#getden" target="kwimg">getden</a> -1, to take
the output density of dataset 1 ;
<li>set <a href="../input_variables/varbas.html#nband" target="kwimg">nband</a> to 8 ;
<li>set <a href="../input_variables/varbas.html#kptopt" target="kwimg">kptopt</a> to -3, to define
three segments in the brillouin Zone ;
<li>set <a href="../input_variables/vargs.html#ndivk" target="kwimg">ndivk</a> to 10 12 17
(this means a circuit defined by 4 points, with 10 divisions of the first
segment, 12 divisions of the second, 17 divisions of the third)
<li>set <a href="../input_variables/vargs.html#kptbounds" target="kwimg">kptbounds</a> to
<pre>
0.5 0.0 0.0 # L point
0.0 0.0 0.0 # Gamma point
0.0 0.5 0.5 # X point
1.0 1.0 1.0 # Gamma point in another cell.
</pre>
<li>set <a href="../input_variables/vargs.html#enunit" target="kwimg">enunit</a> to 1, in order
to have eigenenergies in eV
<li>the only tolerance criterion admitted for non-self-consistent calculations
is <a href="../input_variables/varbas.html#tolwfr" target="kwimg">tolwfr</a>. You should set
it to 1.0d-10 (or so), and suppress <a href="../input_variables/varbas.html#toldfe" target="kwimg">toldfe</a>.
</ul>
<p>The input file ~abinit/tests/tutorial/Input/tbase3_5.in is an example,
while ~abinit/tests/tutorial/Refs/tbase3_5.out is a reference output
file.
<p>You should find the band structure starting at (second dataset):
<pre>
Eigenvalues ( eV ) for nkpt= 40 k points:
Eigenvalues ( eV ) for nkpt= 40 k points:
kpt# 1, nband= 8, wtk= 1.00000, kpt= 0.5000 0.0000 0.0000 (reduced coord)
-3.78815 -1.15872 4.69668 4.69668 7.38795 9.23867 9.23867 13.45707
kpt# 2, nband= 8, wtk= 1.00000, kpt= 0.4500 0.0000 0.0000 (reduced coord)
-3.92759 -0.95774 4.71292 4.71292 7.40692 9.25561 9.25561 13.48927
kpt# 3, nband= 8, wtk= 1.00000, kpt= 0.4000 0.0000 0.0000 (reduced coord)
-4.25432 -0.44393 4.76726 4.76726 7.46846 9.31193 9.31193 13.57737
kpt# 4, nband= 8, wtk= 1.00000, kpt= 0.3500 0.0000 0.0000 (reduced coord)
-4.64019 0.24941 4.85732 4.85732 7.56855 9.38323 9.38323 13.64601
....
</pre>
<p>One needs a graphical tool to represent all these data ...
(For the MAPR 2451 lecture : try with MATLAB).
In a separate file (_EIG), you will find the list of k-points and eigenenergies
(the input variable <a href="../input_variables/varfil.html#prteig" target="kwimg">prteig</a> is set by default to 1).
<p>Even without a graphical tool we will have a quick look
at the values at L, Gamma, X and Gamma again :
<pre>
kpt# 1, nband= 8, wtk= 1.00000, kpt= 0.5000 0.0000 0.0000 (reduced coord)
-3.78815 -1.15872 4.69668 4.69668 7.38795 9.23867 9.23867 13.45707
kpt# 11, nband= 8, wtk= 1.00000, kpt= 0.0000 0.0000 0.0000 (reduced coord)
-6.17005 5.91814 5.91814 5.91814 8.44836 8.44836 8.44836 9.17755
kpt# 23, nband= 8, wtk= 1.00000, kpt= 0.0000 0.5000 0.5000 (reduced coord)
-1.96393 -1.96393 3.00569 3.00569 6.51173 6.51173 15.95524 15.95524
kpt# 40, nband= 8, wtk= 1.00000, kpt= 1.0000 1.0000 1.0000 (reduced coord)
-6.17005 5.91814 5.91814 5.91814 8.44836 8.44836 8.44836 9.17755
</pre>
The last gamma is exactly equivalent to the first gamma.
It can be checked that the top of the valence band is obtained
at Gamma (=5.91814 eV). The width of the valence band is
12.09 eV, the lowest unoccupied state at X is 0.594 eV higher
than the top of the valence band, at Gamma. The Si is described
as an indirect band gap material (this is correct), with a band-gap of
about 0.594 eV (this is quantitatively quite wrong : the
experimental value 1.17 eV is at 25 degree
Celsius). The minimum of the conduction
band is even slightly displaced with respect to X, see kpt # 21 .
This underestimation of the band gap is well-known (the famous DFT band-gap problem).
In order to obtain correct band gaps, you need to go beyond the
Kohn-Sham Density Functional Theory : use the GW approximation.
This is described in <a href="lesson_gw1.html">the first lesson on the GW approximation</a>.
<p>For experimental data and band structure representation, see
<br>M.L. Cohen and J.R. Chelikowski
<br><cite>Electronic structure and optical properties of semiconductors</cite>
<br>Springer-Verlag New-York (1988).
<p>There is a subtlety that is worth to comment about. In non-self-consistent
calculations, like those performed in the present band structure calculation,
with <a href="../input_variables/varbas.html#iscf" target="kwimg">iscf</a>=-2, not all bands
are converged within the tolerance <a href="../input_variables/varbas.html#tolwfr" target="kwimg">tolwfr</a>.
Indeed, the two upper bands (by default) have not been taken into account
to apply this convergence criterion: they constitute a "buffer".
The number of such
"buffer" bands is governed by the input variable
<a href="../input_variables/vargs.html#nbdbuf" target="kwimg">nbdbuf</a>.
<p> It can happen that the highest (or two highest) band(s), if not separated by a gap from
non-treated bands, can exhibit a very slow convergence rate. This buffer
allows to achieve convergence of "important", non-buffer bands.
In the present case, 6 bands have been converged with a residual better
than <a href="../input_variables/varbas.html#tolwfr" target="kwimg">tolwfr</a>, while the
two upper bands are less converged (still sufficiently for graphical representation
of the band structure). In order to achieve the same convergence
for all 8 bands, it is advised to use
<a href="../input_variables/varbas.html#nband" target="kwimg">nband</a>=10 (that is, 8 + 2).
<script type="text/javascript" src="list_internal_links.js"> </script>
</body>
</html>
|
