Elastic and piezoelectric properties
<html> <head>
<title>Tutorial "elastic"</title>
</head> <body bgcolor=“#ffffff”> <hr> <p>This lesson shows how to calculate physical properties related to strain,
for an insulator and a metal : </p>
<ul>
<li>the rigid-atom elastic tensor </li>
<li>the rigid-atom piezoelectric tensor (insulators only) </li>
<li>the internal strain tensor </li>
<li>the atomic relaxation corrections to the elastic and piezoelectric
tensor </li>
</ul>
You should complete lessons <a href="lesson_rf1.html">RF1</a>
and <a href="lesson_rf2.html">RF2</a>
to introduce the response-function features of ABINIT before
starting this lesson. You will learn to use additional response-function
features of ABINIT, and to use relevant parts of the associated
codes Mrgddb and Anaddb.
<p>This lesson should take about two hours. </p>
<h5>Copyright (C) 2000-2013 ABINIT group (DRH) <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>
<h3><b>Content of lesson “elastic”</b></h3>
<ul>
<li><a href="#1">1</a>
The ground-state geometry of (hypothetical)
wurtzite AlAs.</li>
<li><a href="#2">2</a>
Response-function calculation of several second
derivatives of the total energy.</li>
<li><a href="#3">3</a>
anaddb calculation to incorporate atom-relaxation
effects.</li>
<li><a href="#4">4</a>
Finite-difference calculation of elastic and piezoelectric
constants.</li>
<li><a href="#5">5</a>
Alternative response-function calculation of some piezoelectric
constants.</li>
<li><a href="#6">6</a>
Response-function calculation of the elastic
constants for Al metal.</li>
</ul>
<hr> <p><a name=“1”></a>
<br>
</p>
<h3><b>1. The ground-state geometry of (hypothetical) wurtzite AlAs.</b></h3>
<p><i>Before beginning, you might consider working in a different subdirectory
as for the other lessons.
Why not create “Work_elast” in ~abinit/tests/tutorespfn/Input ? </i><br>
</p>
<p>You should copy the files ~abinit/tests/tutorespfn/Input/telast_1.files and telast_1.in
into Work_elast. You may wish to start the calculation (less than one minute on a standard 3GHz machine) before you read the following. You should open your input file telast_1.in with an editor and examine it
as you read this discussion. </p>
<p>The hypothetical wurtzite structure for AlAs retains the tetrahedral coordination
of the atoms of the actual zincblende structure of AlAs,
but has a hexagonal lattice. It was chosen for this lesson because the atomic positions are not completely determined by symmetry. Both the atomic positions and the lattice constants should be optimized before beginning response-function calculations, especially those related to strain properties. While GS structural
optimization was treated in lessons 1-3, we are introducing a few
new features here, and you should look at the following new input variables which will be discussed below:<br>
</p>
<ul>
<li><a href="../input_variables/varrlx.html#getxred" target="kwimg"> getxred</a>
</li>
<li><a href="../input_variables/varrlx.html#iatfix" target="kwimg">iatfix</a>
</li>
<li><a href="../input_variables/varrlx.html#natfix" target="kwimg"> natfix</a>
</li>
<li><a href="../input_variables/varrlx.html#strfact" target="kwimg"> strfact</a>
<br>
</li>
</ul>
<p>There are two datasets specified in telast_1.in. First, let us examine
the common input data. We specify a starting guess for <a href="../input_variables/varbas.html#acell" target="kwimg">
acell</a>
, and give an accurate decimal specification
for <a href="../input_variables/varbas.html#rprim" target="kwimg">
rprim</a>
. The definition of the atom types and
atoms follows <a href=“lesson_rf1.html”>lesson RF1</a>
. The reduced atomic positions <a href="../input_variables/varbas.html#xred" target="kwimg">
xred</a>
are a starting approximation, and will
be replaced by our converged results in the remaining input files, as will <a href=“../input_variables/varbas.html#acell” target=“kwimg”>
acell</a>
.</p>
<p>We will work with a fixed plane wave cutoff <a href=“../input_variables/varbas.html#ecut” target=“kwimg”>
ecut</a>
(=6 Ha), but introduce <a href="../input_variables/varrlx.html#ecutsm" target="kwimg">
ecutsm</a>
(=0.5 Ha)as in <a href="lesson_base3.html">
lesson 3</a>
to smear the cutoff, which produces smoothly varying
stresses as the lattice parameters are optimized. We will keep
the same value of <a href="../input_variables/varrlx.html#ecutsm" target="kwimg">
ecutsm</a>
for the response-function calculations
as well, since changing it from the optimization run value could reintroduce non-zero forces and stresses. For the k-point grid, we must explicitly specify <a href=“../input_variables/varbas.html#shiftk” target=“kwimg”>
shiftk</a>
since the default value results in a grid
shifted so as to break hexagonal symmetry. The RF strain
calculations check this, and will exit with an error message if the grid does not have the proper symmetry. The self-consistency
procedures follow <a href="lesson_rf1.html">lesson RF1</a>
.</p>
<p>Dataset 1 optimizes the atomic positions keeping the lattice parameters
fixed, setting <a href="../input_variables/varrlx.html#ionmov" target="kwimg">
ionmov</a>
=2 as in <a href="lesson_base1.html">lesson
1</a>
. The optimization steps proceed until the maximum force
component on any atom is less than <a href="../input_variables/varrlx.html#tolmxf" target="kwimg">
tolmxf</a>
. It is always advised to relax the forces
before beginning the lattice parameter optimization. Dataset
2 optimizes the lattice parameters with <a href="../input_variables/varrlx.html#optcell" target="kwimg">
optcell</a>
=2 as in <a href="lesson_base3.html">lesson
3</a>
. However, lesson 3 treated cubic Si, and the atom positions in reduced coordinates remained fixed. In the present, more general case, the reduced atomic coordinates must be reoptimized as the lattice parameters are optimized. Note that it is necessary to include
<a href=“../input_variables/varrlx.html#getxred”> getxred</a>
= -1 so that the second dataset is initialized with the relaxed coordinates . Coordinate and lattice parameter optimizations actually take place simultaneously, with the
computed stresses at each step acting as forces on the lattice parameters. We have introduced <a href=“../input_variables/varrlx.html#strfact” target=“kwimg”>
strfact</a>
which scales the stresses so that they
may be compared with the same <a href=“../input_variables/varrlx.html#tolmxf” target=“kwimg”>
tolmxf</a>
convergence test that is applied to the
forces. The default value of 100 is probably a good choice for many systems, but you should be aware of what is happening.</p>
<p>From the hexagonal symmetry, we know that the positions of the atoms in
the a-b basal plane are fixed. However, a uniform translation
along the c axis of all the atoms leaves the structure invariant.
Only the relative displacement of the Al and As planes along
the c axis is physically relevant. We will fix the Al positions to be at reduced c-axis coordinates 0 and 1/2 (these are related by symmetry) by introducing <a href=“../input_variables/varrlx.html#natfix” target=“kwimg”>
natfix</a>
and <a href="../input_variables/varrlx.html#iatfix" target="kwimg">
iatfix</a>
to constrain the structural optimization.
This is really just for cosmetic purposes, since letting them
all slide an arbitrary amount (as they otherwise would) won't
change any results. However, you probably wouldn't want to publish the results that way, so we may as well develop good habits.<br>
</p>
<p>Now we shall examine the results of the structural optimization run. As always, we should first examine the log file to make sure the run has terminated cleanly. There are a number of warnings, but none of them are apparently serious. Next, let us edit the output file, telast_1.out. The first thing to look for is to see whether Abinit recognized the symmetry of the system. In setting up a new data file, it's easy to make mistakes, so this is a valuable check. We see</p>
<pre> DATASET 1 : space group P6_3 m c (#186); Bravais hP (primitive hexag.)<br></pre>
<p>which is correct. Next, we confirm that the structural optimization converged.
The following lines from dataset 1 and dataset2 tell us
that things are OK:</p>
<pre>At Broyd/MD step 4, gradients are converged : <br> max grad (force/stress) = 1.0674E-08 < tolmxf= 1.0000E-06 ha/bohr (free atoms) <br><br>At Broyd/MD step 11, gradients are converged : <br> max grad (force/stress) = 7.8147E-08 < tolmxf= 1.0000E-06 ha/bohr (free atoms) <br></pre>
<p>We can also confirm that the stresses are relaxed: </p>
<pre>Cartesian components of stress tensor (hartree/bohr^3)<br> sigma(1 1)= -3.76644862E-10 sigma(3 2)= 0.00000000E+00<br> sigma(2 2)= -3.76644714E-10 sigma(3 1)= 0.00000000E+00<br> sigma(3 3)= 7.81298436E-10 sigma(2 1)= 0.00000000E+00<br></pre>
<p>Now would be a good time to copy telast_2.in and telast_2.files into your working directory, since we will use the present output to start the next run. Locate the optimized lattice parameters and reduced atomic coordinates near the end of telast_1.out:</p>
<pre> acell2 7.5389648144E+00 7.5389648144E+00 1.2277795374E+01 Bohr<br><br> xred2 3.3333333333E-01 6.6666666667E-01 0.0000000000E+00<br> 6.6666666667E-01 3.3333333333E-01 5.0000000000E-01<br> 3.3333333333E-01 6.6666666667E-01 3.7608588373E-01<br> 6.6666666667E-01 3.3333333333E-01 8.7608588373E-01<br></pre>
<p>With your editor, copy and paste these into telast_2.in at the indicated
places in the "Common input data" area. Be sure to change
acell2 and xred2 to acell and xred since these common values will apply to all datasets in the next set of calculations.<br>
</p>
<hr> <p><a name=“2”></a>
<br>
</p>
<h3><b>2. Response-function calculations of several second derivatives of
the total energy.</b></h3>
<p> We will now compute second derivatives of the total energy (2DTE's) with
respect to all the perturbations we need to compute elastic
and piezoelectric properties. You may want to review <a href=“../users/respfn_help.html#0” target=“kwimg”>
sections 0 and the first paragraph of section
1</a>
of the respfn_help file which you studied
in lesson RF1. We will introduce only one new input variable
for the strain perturbation, </p>
<ul>
<li><a href="../input_variables/varrf.html#rfstrs" target="kwimg">rfstrs</a>
<br>
</li>
</ul>
<p>The treatment of strain as a perturbation has some subtle aspects. It
would be a good idea to read <cite> Metric tensor formulation
of strain in density-functional perturbation theory, by D. R. Hamann,
Xifan Wu, Karin M. Rabe, and David Vanderbilt, <a href="http://prb.aps.org/abstract/PRB/v71/i3/e035117" target="kwimg">Phys. Rev. B 71, 035117 (2005)</a>
</cite> , especially Sec. II and Sec. IV. We will
do all the RF calculations you learned in lesson RF1 together with strain, so you should review the variables</p>
<ul>
<li><a href="../input_variables/varrf.html#rfphon" target="kwimg"> rfphon</a>
</li>
<li><a href="../input_variables/varrf.html#rfatpol" target="kwimg"> rfatpol</a>
</li>
<li><a href="../input_variables/varrf.html#rfdir" target="kwimg"> rfdir</a>
</li>
<li><a href="../input_variables/varrf.html#rfelfd" target="kwimg">rfelfd</a>
<br>
</li>
</ul>
<p>It would be a good idea to copy telast_2.files into Work_elast and start
the calculation while you read (less than 2 minutes on a standard 3GHz machine).
Look at telast_2.in in your editor to follow the discussion, and double check that you have copied acell and xred as discussed in the last section.</p>
<p>This has been set up as a self-contained calculation with three datasets.
The first is simply a GS run to obtain the GS wave functions we will need for the response function (RF) calculations. We
have removed the convergence test from the common input data to remind ourselves that different tests are needed for different datasets.
We set a tight limit on the convergence of the self-consistent potential
with <a href="../input_variables/varbas.html#tolvrs" target="kwimg">
tolvrs</a>
. Since we have specified <a href="../input_variables/varbas.html#nband" target="kwimg">
nband</a>
=8, all the bands are occupied and the potential
test also assures us that all the wave functions are well converged.
This issue will come up again in section <a href="#6">6</a>
. We could have used the output wave functions telast_1o_DS2_WFK
as input for our RF calculations and skipped dataset 1, but redoing
the GS calculation takes relatively little time for this simple
system .</p>
<p>Dataset 2 involves the calculation of the derivatives of the wave functions
with respect to the Brillouin-zone wave vector, the so-called
ddk wave functions. Recall that these are auxiliary quantities
needed to compute the response to the<a href="lesson_rf1.html#5">
electric field perturbation</a>
and introduced in lesson RF1 . It would be a
good idea to review the relevant parts of <a href=“../users/respfn_help.html#1” target=“helpsimg”>
section 1</a>
of the respfn_help file. Examining this
section of telast_2.in, note that electric field as well as strain are uniform perturbations, only are defined for <a href=“../input_variables/vargs.html#qpt” target=“kwimg”>
qpt</a>
= 0 0 0. <a href="../input_variables/varrf.html#rfelfd" target="kwimg">
rfelfd</a>
= 2 specifies that we want the ddk calculation
to be performed, which requires <a href="../input_variables/varbas.html#iscf" target="kwimg">
iscf</a>
= -3. The ddk wave functions will be used
to calculate both the piezoelectric tensor and the Born effective
charges, and in general we need them for <b> k</b> derivatives
in all three (reduced) directions, <a href=“../input_variables/varrf.html#rfdir” target=“kwimg”>
rfdir</a>
= 1 1 1. Since there is no potential self-consistency
in the ddk calculations, we must specify convergence in terms
of the wave function residuals using <a href="../input_variables/varbas.html#tolwfr" target="kwimg">
tolwfr</a>
.<br>
</p>
<p>Finally, dataset 3 performs the actual calculations of the needed 2DTE's
for the elastic and piezoelectric tensors. Setting <a href="../input_variables/varrf.html#rfphon" target="kwimg">
rfphon</a>
= 1 turns on the atomic displacement perturbation,
which we need for all atoms (<a href="../input_variables/varrf.html#rfatpol" target="kwimg">
rfatpol</a>
= 1 4) and all directions (<a href="../input_variables/varrf.html#rfdir" target="kwimg">
rfdir</a>
= 1 1 1). Abinit will calculate first-order
wave functions for each atom and direction in turn, and use
those to calculate 2DTE's with respect to all pairs of atomic displacements and with respect to one atomic displacement and one component of electric field. These quantities, the interatomic force constants (at gamma) and the Born effective charges will be used later to compute the atomic relaxation contribution to the elastic
and piezoelectric tensor.</p>
<p>First-order wave functions for the strain perturbation are computed next.
Setting <a href="../input_variables/varrf.html#rfstrs" target="kwimg">rfstrs</a>
= 3 specifies that we want both uniaxial
and shear strains to be treated, and <a href=“../input_variables/varrf.html#rfdir” target=“kwimg”>
rfdir</a>
= 1 1 1 cycles through strains xx, yy, and
zz for uniaxial and yz, xz, and xy for shear. We note that while
other perturbations in Abinit are treated in reduced coordinates, strain is better dealt with in Cartesian coordinates for reasons discussed in the reference cited above. These wave functions
are used to compute three types of 2DTE's. Derivatives with respect
to two strain components give us the so-called rigid-ion elastic
tensor. Derivatives with respect to one strain and one electric field
component give us the rigid-ion piezoelectric tensor. Finally, derivatives with respect to one strain and one atomic displacement yield the
internal-strain force-response tensor, an intermediate quantity that will be necessary to compute the atomic relaxation corrections to the rigid-ion quantities. As in lesson RF1, we specify convergence in terms of the residual of the potential (here the first-order potential)
using <a href="../input_variables/varbas.html#tolvrs" target="kwimg">
tolvrs</a>
.</p>
<p>Your run should have completed by now. Abinit should have created quite a few files.</p>
<ul>
<li>telast_2.log (log file)</li>
<li>telast_2.out (main output file)</li>
<li>telast_2o_DS1_DDB (first derivatives of the energy from GS calculation)</li>
<li>telast_2o_DS3_DDB (second derivatives from the RF calculation)</li>
<li>telast_2o_DS1_WFK (GS wave functions)</li>
<li>telast_2o_DS2_1WF* (ddk wave functions)</li>
<li>telast_2o_DS3_1WF* (RF first-order wave functions from various perturbations)</li>
</ul>
The log and out files are diagnostics and
readable output information for a wide variety of properties. The derivative database DDB files are ascii and readable, but primarily for subsequent analysis by anaddb which we will undertake in the next section. Finally, the various wave function binary files are primarily of use for subsequent calculations, where they could cut the number of needed iterations in, for example, convergence testing.
We take note of a few conventions in the file names. The root output
file name telast_2o is from the 4th line of the “files” file. The dataset producing the file is next. Finally, the first-order wave function 1WF files have a final “pertcase” number described in <a href=“../users/respfn_help.html#1” target=“kwimg”>
section 1</a>
of the respfn_help file. While telast_2.in
specifies all atomic displacements, only the symmetry-inequivalent perturbations are treated, so the “pertcase” list is incomplete.
All cases specified in the input data are treated for the strain
perturbation. <br>
<p>First, take a look at the end of the telast_2.log file to make sure the
run has completed without error. You might wish to take a
look at the WARNING's, but they all appear to be harmless. Next, edit your telast_2.out file. Searching backwards for ETOT you will find</p>
<pre> iter 2DEtotal(Ha) deltaE(Ha) residm vres2<br>-ETOT 1 2.3955210936959 -6.519E+00 6.313E-01 4.126E+02<br> ETOT 2 1.3034866746082 -1.092E+00 4.874E-04 4.710E+00<br> ETOT 3 1.2898922627828 -1.359E-02 1.856E-05 3.514E-01<br> ETOT 4 1.2891989643464 -6.933E-04 2.648E-07 1.382E-02<br> ETOT 5 1.2891785442295 -2.042E-05 8.156E-09 1.945E-04<br> ETOT 6 1.2891783810507 -1.632E-07 6.814E-11 4.395E-05<br> ETOT 7 1.2891783087573 -7.229E-08 2.750E-11 3.704E-06<br> ETOT 8 1.2891783033031 -5.454E-09 2.224E-12 1.001E-07<br> ETOT 9 1.2891783031248 -1.783E-10 7.592E-14 6.584E-10<br> ETOT 10 1.2891783031235 -1.276E-12 1.112E-15 4.697E-11<br><br> At SCF step 10 vres2 = 4.70E-11 < tolvrs= 1.00E-10 =>converged.<br></pre>
<p>Abinit is solving a set of Schroedinger-like equations for the first-order
wave functions, and these functions minimize a variational expression
for the 2DTE. (Technically, they are called self-consistent Sternheimer
equations.) The energy convergence looks similar to that
of GS calculations. The fact that vres2, the residual of the self-consistent
first-order potential, has reached <a href="../input_variables/varbas.html#tolvrs" target="kwimg">
tolvrs</a>
well within <a href="../input_variables/varbas.html#nstep" target="kwimg">nstep</a>
(40) iterations indicates that the 2DTE calculation for this perturbation
(xy strain) has converged . It would pay to examine a few
more cases for different perturbations (unless you have looked
through all the warnings in the log). </p>
<p>Another convergence item to examine in your .out file is</p>
<pre> Seventeen components of 2nd-order total energy (hartree) are<br> 1,2,3: 0th-order hamiltonian combined with 1st-order wavefunctions<br> kin0= 9.10477366E+00 eigvalue= 3.11026172E-01 local= -3.66858410E+00<br> 4,5,6,7: 1st-order hamiltonian combined with 1st and 0th-order wfs<br> loc psp = -8.91644866E+00 Hartree= 4.33575581E+00 xc= -6.58530125E-01<br> kin1= -8.62111357E+00<br> 8,9,10: eventually, occupation + non-local contributions<br> edocc= 0.00000000E+00 enl0= 6.43290213E-01 enl1= -1.55388913E-01<br> 1-10 gives the relaxation energy (to be shifted if some occ is /=2.0)<br> erelax= -7.62521951E+00<br> 11,12,13 Non-relaxation contributions : frozen-wavefunctions and Ewald<br> fr.hart= -1.18530360E-01 fr.kin= 5.20015318E+00 fr.loc= 4.18792396E-01<br> 14,15,16 Non-relaxation contributions : frozen-wavefunctions and Ewald<br> fr.nonl= 2.94970653E-01 fr.xc= 9.41457939E-02 Ewald= 3.02486615E+00<br> 17 Non-relaxation contributions : pseudopotential core energy<br> pspcore= 0.00000000E+00<br> Resulting in :<br> 2DEtotal= 0.1289178303E+01 Ha. Also 2DEtotal= 0.350803264954E+02 eV<br> (2DErelax= -7.6252195079E+00 Ha. 2DEnonrelax= 8.9143978110E+00 Ha)<br> ( non-var. 2DEtotal : 1.2891783532E+00 Ha)<br></pre>
<p>This detailed breakdown of the contributions to 2DTE is probably
of limited interest, but you should compare "2DEtotal" and "non-var.
2DEtotal“ from the last three lines. While the first-order wave function for the present perturbation minimizes a variational expression for the second derivative with respect to this perturbation as we just saw, the various 2DTE given as elastic tensors, etc. in the output and in the DDB file are all computed using non-variational expressions.
Using the non-variational expressions, mixed second derivatives
with respect to the present perturbation and all other perturbations of interest can be computed directly from the present first-order wave functions. The disadvantage is that the non-variational result has errors which are linearly proportional to convergence errors in the GS and first-order wave functions. Since errors in the variational 2DEtotal are second-order
in wave-function convergence errors, comparing this to the non-variational
result for the diagonal second derivative will give an idea of the
accuracy of the latter and perhaps indicate the need for tighter convergence
tolerances for both the GS and RF wave functions. This is discussed
in <cite> X. Gonze and C. Lee, Phys. Rev. B 55, 10355 (1997)</cite>
, Sec. II. For an atomic-displacement perturbation, the
corresponding breakdown of the 2DTE is headed “Thirteen components.”</p>
<p>Now let us take a look at the results we want, the various 2DTE's.
They begin</p>
<pre> ==> Compute Derivative Database <==<br> <br> 2nd-order matrix (non-cartesian coordinates, masses not included,<br> asr not included )<br> cartesian coordinates for strain terms (1/ucvol factor <br> for elastic tensor components not included) <br> j1 j2 matrix element<br> dir pert dir pert real part imaginary part<br> <br> 1 1 1 1 5.4508667670 0.0000000000<br> 1 1 2 1 -2.7254333834 0.0000000000<br> 1 1 3 1 0.0000000000 0.0000000000<br> …..<br></pre>
<p>These are the “raw” 2DTE's, in reduced coordinates for atom-displacement
and electric-field perturbations, but Cartesian coordinates
for strain perturbations. This same results with the same organization
appear in the file telast_2_DS3_DDB which will be used later
as input for automated analysis and converted to more useful notation
and units by anaddb. A breakout of various types of 2DTE's follows
(all converted to Cartesian coordinates and in atomic units): <br>
</p>
<pre> Dynamical matrix, in cartesian coordinates,<br> if specified in the inputs, asr has been imposed<br> j1 j2 matrix element<br> dir pert dir pert real part imaginary part<br> <br> 1 1 1 1 0.0959051953 0.0000000000<br> 1 1 2 1 0.0000000000 0.0000000000<br> 1 1 3 1 0.0000000000 0.0000000000<br> …..<br></pre>
<p>This contains the interatomic force constant data that will be used later
to include atomic relaxation effects. "asr" refers to
the acoustic sum rule, which basically is a way of making sure that forces
sum to zero when an atom is displaced.</p>
<pre> Effective charges, in cartesian coordinates,<br> (from phonon response) <br> if specified in the inputs, asr has been imposed<br> j1 j2 matrix element<br> dir pert dir pert real part imaginary part<br> <br> 1 6 1 1 1.8290468197 0.0000000000<br> 2 6 1 1 0.0000000000 0.0000000000<br> 3 6 1 1 0.0000000000 0.0000000000<br> …..<br></pre>
<p>The Born effective charges will be used to find the atomic relaxation contributions of the piezoelectric tensor.</p>
<pre> Rigid-atom elastic tensor , in cartesian coordinates,<br> j1 j2 matrix element<br> dir pert dir pert real part imaginary part<br> <br> 1 7 1 7 0.0056418398 0.0000000000<br> 1 7 2 7 0.0013753713 0.0000000000<br> 1 7 3 7 0.0007168444 0.0000000000<br> …..<br></pre>
<p><font face=“Times New Roman, Times, serif”>The rigid-atom elastic tensor
is the 2DTE with respect to a pair of strains. We recall that "pert" = natom+3 and natom+4 for unaxial and shear strains, respectively.</font></p>
<pre> Internal strain coupling parameters, in cartesian coordinates,<br> zero average net force deriv. has been imposed <br> j1 j2 matrix element<br> dir pert dir pert real part imaginary part<br> <br> 1 1 1 7 0.1249319229 0.0000000000<br> 1 1 2 7 -0.1249319273 0.0000000000<br> 1 1 3 7 0.0000000000 0.0000000000<br> …..<br></pre>
<p>These 2DTE's with respect to one strain and one atomic displacement are
needed for atomic relaxation corrections to both the elastic tensor and piezoelectric tensor. While this set of parameters is of limited direct interest, it should be examined in cases when you think that
high symmetry may eliminate the need for these corrections. You are
probably wrong, and any non-zero term indicates a correction.</p>
<pre> Rigid-atom proper piezoelectric tensor, in cartesian coordinates,<br> j1 j2 matrix element<br> dir pert dir pert real part imaginary part<br> <br> 1 6 1 7 0.0000000000 0.0000000000<br> 1 6 2 7 0.0000000000 0.0000000000<br> 1 6 3 7 0.0000000000 0.0000000000<br></pre>
<p>Finally, we have the piezoelectric tensor, the 2DTE with respect to one
strain and one uniform electric field component. (Yes,
there are non-zero elements.)<br>
</p>
<hr> <p><a name=“3”></a>
<br>
</p>
<h3><b>3. anaddb calculation of atom-relaxation effects.</b></h3>
<p>In this section, we will run the program anaddb, which analyzes DDB files
generated in prior RF calculations. You should copy telast_3.in and telast_3.files in your Work_elast directory. You should
now go to the <a href=”../users/anaddb_help.html“ target=“kwimg”>
anaddb help file</a>
, and read the short introduction. The bulk
of the material in this help file is contained in the description
of the variables. You should read the descriptions of<br>
</p>
<ul>
<li><a href="../users/anaddb_help.html#elaflag" target="kwimg">elaflag</a>
</li>
<li><a href="../users/anaddb_help.html#piezoflag" target="kwimg">piezoflag</a>
</li>
<li><a href="../users/anaddb_help.html#instrflag" target="kwimg">instrflag</a>
<br>
</li>
<li><a href="../users/anaddb_help.html#chneut" target="kwimg">chneut</a>
</li>
</ul>
For the theory underlying the incorporation
of atom-relaxation corrections, it is recommended you see X. Wu, D. Vanderbilt, and D. R. Hamann, <a href=“http://prb.aps.org/abstract/PRB/v72/i3/e035105” target=“kwimg”>Phys. Rev, B 72, 035105 (2005)</a> .<br>
<br>
Anaddb can do lots of other things, such as
calculate the frequency-dependent dielectric tensor, interpolate the phonon spectrum to make nice phonon dispersion plots, calculate Raman spectra, etc., but we are focusing on the minimum needed for the elastic and piezoelectric constants at zero electric field.
<br>
<br>
We also mention that <a href="../users/mrgddb_help.html" target="helpsimg">
mrgddb</a>
is another utility program that can be used
to combine DDB files generated in several different datasets or in different runs into a single DDB file that can be analyzed by anaddb. One particular usage would be to combine the DDB file produced by the GS run, which contains first-derivative information
such as stresses and forces with the RF DDB. It is anticipated that
anaddb in a future release will implement the finite-stress corrections
to the elastic tensor discussed in <a href="../theory/elasticity-oganov.pdf" target="kwimg">
notes by A. R. Oganov</a>
.<br>
<br>
Now would be a good time to edit telast_3.in
and observe that it is very simple, consisting of nothing more than the four variables listed above set to appropriate values.
The telast_3.files file is used with anaddb in the same manner as
the abinit .files you are by now used to. The first two lines specify
the .in and .out files, the third line specifies the DDB file, and
the last two lines are dummy names which would be used in connection with other capabilities of anaddb. Now you should run the calculation, which is done in the same way as you are now used to for abinit:<br>
<br>
<small><font face="Courier New, Courier, monospace">
../../anaddb <telast_3.files >&telast_3.log</font></small><br>
<br>
This calculation should only take a few seconds.
You should edit the log file, go to the end, and make sure the
calculation terminated without error. Next, examine telast_3.out.
After some header information, we come to tables giving the "force-response"
and "displacement-response" internal strain tensors. These represent,
respectively, the force on each atom and the displacement of each
atom in response to a unit strain of the specified type. These numbers are of limited interest to us, but represent important intermediate quantities in the treatment of atomic relaxation (see the X. Wu paper cited above). <br>
<br>
Next, we come to the elastic tensor output:<br>
<br>
<pre> Elastic Tensor(clamped ion)(unit:10^2GP):<br><br> 1.6598864 0.4046482 0.2109029 0.0000000 0.0000000 0.0000002<br> 0.4046481 1.6598863 0.2109029 0.0000000 0.0000000 0.0000002<br> 0.2109030 0.2109030 1.8258574 0.0000000 0.0000000 0.0000002<br> 0.0000000 0.0000000 0.0000000 0.4081819 0.0000000 0.0000000<br> 0.0000000 0.0000000 0.0000000 0.0000000 0.4081822 0.0000000<br> 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.6276191<br><br> Elastic Tensor(relaxed ion)(unit:10^2GP):<br>(at fixed electric field boundary condition)<br><br> 1.3526230 0.5445033 0.3805291 0.0000000 0.0000000 0.0000002<br> 0.5445032 1.3526228 0.3805291 0.0000000 0.0000000 0.0000002<br> 0.3805292 0.3805293 1.4821105 0.0000000 0.0000000 0.0000002<br> 0.0000000 0.0000000 0.0000000 0.3055073 0.0000000 0.0000000<br> 0.0000000 0.0000000 0.0000000 0.0000000 0.3055072 0.0000000<br> 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.4040599<br><br></pre>
<br>
While not labeled, the rows and columns 1-6 here
represent xx, yy, zz, yz, xz, xy strains and stresses in the conventional Voigt notation. <br> The clamped-ion results were calculated in the telast_2 RF run, and are simply converted to standard GPa units by anaddb (the terms “clamped ion,” “clamped atom,” and “rigid atom” used in various places are interchangeable, similarly for “relaxed.”) <br>
The relaxed-ion result was calculated by anaddb by combining 2DTE's
for internal strain and interatomic force constants which are stored in the input DDB file. Comparing the clamped and relaxed results, we see that all the diagonal elastic constants have decreased in value. <br> This is plausible, since allowing the internal degrees of freedom to relax should make a material less stiff. These tensors should be symmetric, and certain tensor elements should be zero or identical by symmetry. <br> It's a good idea to check these properties against a standard text such as J. F. Nye, Physical Properties of Crystals (Oxford U. P., Oxford 1985). Departures from expected symmetries (there are a few in the last decimal place here) are due to either convergence
errors or, if large, incorrectly specified geometry.<br>
<br>
Next in telast_3.out we find the piezoelectric tensor
results:<br>
<br>
<pre> Proper piezoelectric constants(clamped ion)(unit:c/m^2)<br><br> 0.00000000 0.00000000 0.38490082<br> 0.00000000 0.00000000 0.38490078<br> 0.00000000 0.00000000 -0.73943037<br> 0.00000000 0.43548809 0.00000000<br> 0.43548801 0.00000000 0.00000000<br> 0.00000000 0.00000000 0.00000000<br><br> Proper piezoelectric constants(relaxed ion)(unit:c/m^2)<br><br> 0.00000000 -0.00000002 -0.01187139<br> 0.00000000 0.00000002 -0.01187157<br> 0.00000000 0.00000000 0.06462740<br> 0.00000000 -0.04828794 0.00000000<br> -0.04828881 0.00000000 0.00000000<br> 0.00000001 -0.00000001 0.00000000<br></pre>
<br>
The 3 columns here represent x, y, and z electric
polarization, and the 6 rows the Voigt strains. The clamped-ion result was calculated in the telast_2 RF run, and is simply scaled to conventional units by anaddb. The ion relaxation contributions are based on 2DTE's for internal strain, interatomic force constants, and Born effective charges, and typically constitute much larger corrections to the piezoelectric tensor than to the elastic tensor.
Once again, symmetries should be checked. (The slight discrepancies
seen here can be removed by setting tolvrs3=1.0d-18 in telast_2.in.)
One should be aware that the piezoelectric tensor is identically zero
in any material which has a center of symmetry.<br>
<br>
Since we are dealing with a hypothetical material,
there is no experimental data with which to compare our results.
In the next section, we will calculate a few of these numbers by
a finite-difference method to gain confidence in the RF approach.<br>
<br>
<hr> <p><a name=“4”></a>
<br>
</p>
<h3><b>4. Finite-difference calculation of elastic and piezoelectric constants.</b></h3>
You should copy telast_4.in and telast_4.files into your Work_elast directory. Editing telast_4.in, you will see
that it has four datasets, the first two with the c-axis contracted 0.01% and the second two with it expanded 0.01%, which we specified
by changing the third row of <a href="../input_variables/varbas.html#rprim" target="kwimg">
rprim</a>
. The common data is essentially the same as telast_2.in,
and the relaxed <a href="../input_variables/varbas.html#acell" target="kwimg">
acell</a>
values and <a href="../input_variables/varbas.html#xred" target="kwimg">
xred</a>
from telast_1.out have already been included.
Datasets 1 and 3 do the self-consistent convergence of the GS wave functions for the strained lattices and compute the stress. Datasets 2 and 4 introduce a new variable.<br>
<ul>
<li><a href="../input_variables/varff.html#berryopt" target="kwimg">berryopt</a>
</li>
</ul>
Electric polarization in solids is a subtle topic
which has only recently been rigorously resolved. It is now understood to be a bulk property, and to be quantitatively described by a Berry phase formulation introduced by R. D. King-Smith and D. Vanderbilt, Phys. Ref. B 47, 1651(1993) . It can be calculated in a GS calculation by integrating the gradient with respect to <b>k</b> of the GS wave functions over the Brillouin zone. In GS calculations, the gradients are approximated by finite-difference
expressions constructed from neighboring points in the <b>k</b> mesh.
These are closely related to the ddk wave functions used in RF calculations in <a href=”#2“> 2</a>
and introduced in <a href="lesson_rf1.html#5" target="kwimg">
lesson RF1, section 5</a>
. We will use <a href="../input_variables/varff.html#berryopt" target="kwimg">
berryopt</a>
= -1, which utilizes an improved coding of the calculation,
and must specify <a href="../input_variables/varrf.html#rfdir" target="kwimg">rfdir</a>
= 1 1 1 so that the Cartesian components of the polarization
are computed.<br>
<br>
Now, run the telast_4 calculation, which should only
take a minute or two, and edit telast_4.out. To calculate the
elastic constants, we need to find the stresses <small><font face=“Courier New, Courier, monospace”>
sigma(1 1)</font></small> and<small><font face="Courier New, Courier, monospace">
sigma(3 3)</font></small> . We see that each of
the four datasets have stress results, but that there are slight differences between those from, for example dataset 1 and dataset 2, which should be identical. Despite our tight limit, this is still a convergence issue. Look at the following convergence results, <br>
<br>
<pre>Dataset 1:<br> At SCF step 21 vres2 = 4.14E-19 < tolvrs= 1.00E-18 =>converged.<br><br>Dataset 2:<br> At SCF step 1 vres2 = 8.54E-20 < tolvrs= 1.00E-18 =>converged.<br></pre>
<p>Since dataset 2 has better convergence, we will use this and the dataset
4 results, choosing those in GPa units,</p>
<pre>- sigma(1 1)= -2.11921106E-03 sigma(3 2)= 0.00000000E+00<br>- sigma(3 3)= -1.82392096E-02 sigma(2 1)= 0.00000000E+00<br><br>- sigma(1 1)= 2.09884071E-03 sigma(3 2)= 0.00000000E+00<br>- sigma(3 3)= 1.82778626E-02 sigma(2 1)= 0.00000000E+00<br></pre>
<p>Let us now compute the numerical derivative of <small><font face=“Courier New, Courier, monospace”>
sigma(3 3)</font></small>and compare to our RF result. Recalling that our dimensionless strains were ±0.0001,
we find 182.5853 GPa. This compares very well with the value <font face=“Times New Roman, Times, serif”>
182.58574</font> GPa, the 3,3 element of the Rigid-ion
elastic tensor we found from our anaddb calculation in <a href="#3">
3</a>
. (Recall that our strain and stress were both 3
3 or z z or Voigt 3.) Similarly, the numerical derivative of <small><font face=“Courier New, Courier, monospace”>
sigma(1 1)</font></small>is 21.09025 GPa, compared
to <font face=“Times New Roman, Times, serif”> 21.09030</font>
GPa, the 3,1 elastic-tensor element.</p>
<p>The good agreement we found from this simple numerical differentiation
required that we had accurately relaxed the lattice so that the
stress of the unstrained structure was very small. Similar numerical-derivative
comparisons for systems with finite stress are more complicated,
as discussed in <a href=”../theory/elasticity-oganov.pdf“ target=“kwimg”>
notes by A. R. Oganov</a>
. Numerical-derivative comparisons for the relaxed-ion
results are extremely challenging since they require relaxing atomic
forces to exceedingly small limits.</p>
<p>Now let us examine the electric polarizations found in datasets 2 and 4, focusing on the C/m^2 results,</p>
<pre> Polarization -1.578184222E-11 C/m^2<br> Polarization 1.578180951E-11 C/m^2<br> Polarization -2.979936117E-01 C/m^2<br><br> Polarization -1.577713293E-11 C/m^2<br> Polarization 1.577662674E-11 C/m^2<br> Polarization -2.981427295E-01 C/m^2<br></pre>
<p>While not labeled as such, these are the Cartesian x, y, and z
components, respectively, and the x and y components are zero
within numerical accuracy as they must be from symmetry. Numerical differentiation of the z component yields -0.745589 C/m^2. This is to be compared with the z,3 element of our rigid-ion piezoelectric tensor from <a href=”#3“> 3</a>
, <font face="Times New Roman, Times, serif">-0.73943037</font>
C/m^2, and the two results do not compare as well as we
might wish.</p>
<p>What is wrong? There are two possibilities. The first is that the RF
calculation produces the proper piezoelectric tensor, while numerical
differentiation of the polarization produces the improper piezoelectric
tensor. This is a subtle point, for which you are referred to
D. Vanderbilt, J. Phys. Chem. Solids 61, 147 (2000)
. The improper-to-proper transformation only effects certain
tensor elements, however, and for our particular combination of crystal symmetry and choice of strain there is no correction. The second possibility is the subject of the next section.</p>
<p><br>
</p>
<hr> <p><a name=“5”></a>
<br>
</p>
<h3><b>5. Alternative response-function calculation of some piezoelectric
constants.</b></h3>
<p>Our GS calculation of the polarization in <a href=”#4“>4</a>
used, in effect, a finite-difference approximation
to ddk wave functions, while our RF calculations in <a href=”#2“>
2</a>
used analytic results based on the RF approach.
Since the <b> k</b> grid determined by <a href=”../input_variables/varbas.html#ngkpt“ target=“kwimg”>
ngkpt</a>
= 4 4 4 and <a href="../input_variables/varbas.html#nshiftk" target="kwimg">
nshiftk</a>
= 1 is rather coarse, this is a probable source of
discrepancy. Since this issue was noted previously in connection with the calculation of Born effective charges by Na Sai, K. M. Rabe, and D. Vanderbilt, Phys. Rev. B 66, 104108 (2002)
, Abinit has incorporated the ability to use finite-difference ddk
wave functions from GS calculations in RF calculations of electric-field-related
2DTE's. Copy telast_5.in and telast_5.files into Work_elast, and edit telast_5.in.</p>
<p>You should compare this with our previous RF data, telast_2.in, and note
that dataset1 and the Common data (after entering relaxed structural results) are essentially identical. Dataset 2 has been replaced
by a non-self-consistent GS calculation with <a href=”../input_variables/varff.html#berryopt“ target=“kwimg”>
berryopt</a>
= -2 specified to perform the finite-difference ddk wave function
calculation. (The finite-difference first-order wave functions
are implicit but not actually calculated in the GS polarization calculation.)
We have restricted <a href="../input_variables/varrf.html#rfdir" target="kwimg">
rfdir</a>
to 0 0 1 since we are only interested in the 3,3
piezoelectric constant. Now compare dataset 3 with that in telast_2.in.
Can you figure out what we have dropped and why? Run the telast_5
calculation, which will only take about a minute with our simplifications.</p>
<p>Now edit telast_5.out, looking for the piezoelectric tensor,</p>
<pre> Rigid-atom proper piezoelectric tensor, in cartesian coordinates,<br> j1 j2 matrix element<br> dir pert dir pert real part imaginary part<br> <br> 3 6 3 7 -0.0130314055 0.0000000000<br></pre>
<br>
We immediately see a problem -- this output, like
most of the .out file, is in atomic units, while we computed our numerical derivative in conventional C/m^2 units. While you might think to simply run anaddb to do the conversion as before, its present version is not happy with such an incomplete DDB file as telast_5 has generated and will not produce the desired result. While it should be left as an exercise to the student to dig the conversion factor out of the literature, or better yet out of the source code, we will cheat
and tell you that 1a.u.=57.2147606 C/m^2 Thus the new RF result for the
3,3 rigid-ion piezoelectric constant is -0.7455887 C/m^2 compared to the result found in <a href=”#4“> 4</a>
by a completely-GS finite difference calculation,
-0.745589 C/m^2. The agreement is now excellent!
<p></p>
<p>The fully RF calculation in <a href=”#2“>2</a>
in fact will converge much more rapidly with <b>
k</b> sample than the partial-finite-difference method introduced here. Is it worthwhile to have learned how to do this? We believe that is always pays to have alternative ways to test results, and besides, this didn't take much time. (Have you found the conversion factor on your own yet?)<br>
</p>
<p><br>
</p>
<hr> <p><a name=“6”></a>
<br>
</p>
<h3><b>6. Response-function calculation of the elastic constants of Al metal.</b></h3>
For metals, the existence of partially occupied bands is
a complicating feature for RF as well as GS calculations. Now would be a good time to review <a href=“lesson_base4.html”>lesson 4</a> which dealt in detail with the interplay between <b> k</b>-sample convergence and Fermi-surface broadening, especially section <a href=“lesson_base4.html#43”> 4.3</a>
. You should copy telast_6.in and telast_6.files into Work_elast, and begin your run while you read on, since it involves
a convergence study with multiple datasets and may take about two minutes.<br>
<br>
While the run is in progress, edit telast_6.in. As
in t43.in, we will set <a href=”../input_variables/varbas.html#udtset“ target=“kwimg”>
udtset</a>
to specify a double loop. In the present case, however,
the outer loop will be over 3 successively larger meshes of <b>k</b>
points, while the inner loop will be successively<br>
<ol>
<li>GS self-consistent runs with optimization of acell.</li>
<li>GS density-generating run for the next step.</li>
<li>Non-self-consistent GS run to converge unoccupied or slightly-occupied
bands.</li>
<li>RF run for symmetry-inequivalent elastic constants.</li>
</ol>
In Section<a href="#1"> 1</a>
, we did a separate GS structural optimization run and transferred
the results by hand to RF run <a href="#2">2</a>
. Because we are doing a convergence test here, we
have combined these steps, and use <a href=”../input_variables/varrlx.html#getcell“ target=“kwimg”>
getcell</a>
to transfer the optimized coordinates from the first dataset
of the inner loop forward to the rest. If we were doing a more
complicated structure with internal coordinates that were also optimized,
we would need to use both this and <a href="../input_variables/varrlx.html#getxred" target="kwimg">
getxred</a>
to transfer these, as in telast_1.in.<br>
<br>
The specific data for inner-loop dataset 1 is very similar
to that for telast_1.in. Inner-loop dataset 2 is a bit of a hack.
We need the density for inner-loop dataset 3, and while we could
set <a href="../input_variables/varfil.html#prtden" target="kwimg"> prtden</a>
= 1 in dataset 1, this would produce a separate density file
for every step in the structural optimization, and it isn't clear how to automatically pick out the last one. So, dataset 2 picks up the wave functions from dataset 1 (only one file of these is produced, for the optimized structure), does one more iteration with fixed geometry,
and writes a density file. <br>
<br>
Inner-loop dataset 3 is a non-self-consistent run
whose purpose is to ensure that all the wave functions specified by <a href=”../input_variables/varbas.html#nband“ target=“kwimg”> nband</a>
are well converged. For metals, we have to specify enough
bands to make sure that the Fermi surface is properly calculated. Bands
above the Fermi level which have small occupancy or near-zero
occupancy if their energies exceed the Fermi energy by more than a few times<a href=”../input_variables/vargs.html#tsmear“ target=“kwimg”> tsmear</a>
, will have very little effect on the self-consistent potential,
so the <a href="../input_variables/varbas.html#tolvrs" target="kwimg"> tolvrs</a>
test in dataset 1 doesn't ensure their convergence. Using
<a href="../input_variables/varbas.html#tolwfr" target="kwimg"> tolwfr</a>
in inner-loop dataset 3 does. Partially-occupied
or unoccupied bands up to <a href=”../input_variables/varbas.html#nband“ target=“kwimg”>
nband</a>
play a different role in constructing the first-order
wave functions than do the many unoccupied bands beyond <a href="../input_variables/varbas.html#nband" target="kwimg">
nband</a>
which aren't explicitly treated in Abinit, as discussed
in S. de Gironcoli, Phys. Rev. B 51, 6773 (1995). By
setting <a href="../input_variables/varbas.html#nband" target="kwimg"> nband</a>
exactly equal to the number of occupied bands for RF calculations
for semiconductors and insulators, we avoid having to deal with the
issue of converging unoccupied bands. Could we avoid the extra
steps by simply using <a href="../input_variables/varbas.html#tolwfr" target="kwimg">
tolwfr</a>
instead of <a href="../input_variables/varbas.html#tolvrs" target="kwimg">
tolvrs</a>
in dataset 1? Perhaps, but experience has shown that
this does not necessarily lead to as well-converged a potential, and
it is not recommended. These same considerations apply to phonon
calculations for metals, or in particular to <a href="../input_variables/vargs.html#qpt">
qpt</a>
= 0 0 0 phonon calculations for the interatomic force constants
needed to find atom-relaxation contributions to the elastic constants
for non-trivial structures as in <a href="#2">2</a>
and <a href="#3">3</a>
.<br>
<br>
The data specific to the elastic-tensor RF calculation in
inner-loop dataset 4 should by now be familiar. We take advantage of the fact that for cubic symmetry the only symmetry-inequivalent elastic constants are C<small><small> 11</small></small>, C<small><small> 12</small></small> , and C<small><small> 44</small></small> . Abinit, unfortunately, does not do this analysis automatically, so we specify <a href=”../input_variables/varrf.html#rfdir“ target=“kwimg”> rfdir</a>
=1 0 0 to avoid duplicate calculations. (Note that
if atom relaxation is to be taken into account for a more complex
structure, the full set of directions must be used.)<br>
<br>
When the telast_6 calculations finish, first look at telast_6.log
as usual to make sure they have run to completion without error. Next,
it would be a good idea to look at the band occupancies occ?? (where
?? is a dual-loop dataset index) reported at the end (following <small><font face=“Courier New, Courier, monospace”>
==END DATASET(S)==</font></small>). The highest band,
the fourth in this case, should have zero or very small occupation,
or you need to increase <a href="../input_variables/varbas.html#nband" target="kwimg">
nband</a>
or decrease <a href="../input_variables/vargs.html#tsmear" target="kwimg">
tsmear</a>
. Now, use your newly perfected knowledge of the Abinit
perturbation indexing conventions to scan through telast_6.out and
find C<small><small> 11</small></small> , C<small><small>12</small></small>
, and C<small><small> 44</small></small> for each of the
three <b>k</b>-sample choices, which will be under the ”<small><font face=“Courier New, Courier, monospace”>
Rigid-atom elastic tensor</font></small>" heading. Also
find the lattice constants for each case, whose convergence you studied in lesson 4. You should be able to cut-and-paste these into a table like the following,<br>
<pre> C_11 C_12 C_44 acell<br><br>ngkpt=3*6 0.0037773594 0.0022583552 0.0013453703 7.5710952267<br>ngkpt=3*8 0.0042004471 0.0020423400 0.0013076775 7.5693986688<br>ngkpt=3*10 0.0042034439 0.0020343450 0.0012956781 7.5694820863<br></pre>
<p>We can immediately see that the lattice constant converges considerably
more rapidly with <b>k</b> sample than the elastic constants. For
<a href="../input_variables/varbas.html#ngkpt" target="kwimg"> ngkpt</a>
=3*6, acell is converged to 0.02%, while the C's have 5-10%
errors. For <a href=“../input_variables/varbas.html#ngkpt” target=“kwimg”>ngkpt</a>
=3*8, the C's are converged to better than 1%, much better
for the largest, C<small><small>11</small></small>, which should be acceptable.</p>
<p>As in lesson 4, the <a href=“../input_variables/varbas.html#ngkpt” target=“kwimg”>ngkpt</a>
convergence is controlled by <a href="../input_variables/vargs.html#tsmear" target="kwimg">
tsmear</a>
. The smaller the broadening, the denser the <b>k</b>
sample that is needed to get a smooth variation of occupancy, and
presumably stress, with strain. While we will not explore <a href=“../input_variables/vargs.html#tsmear” target=“kwimg”>
tsmear</a>
convergence in this lesson, you may wish to do so on your
own. We believe that the value <a href=“../input_variables/vargs.html#tsmear” target=“kwimg”>
tsmear</a>
= 0.02 in telast_6.in gives results within 1% of the
fully-converged small-broadening limit.</p>
<p><b>We find that</b> <b><a href=“../input_variables/varbas.html#occopt” target=“kwimg”>
occopt</a>
</b><b>=3, standard Fermi-Dirac broadening</b><b>, gives
<u> much better</u> convergence of the C's than “cold smearing.”</b>
Changing <a href="../input_variables/varbas.html#occopt" target="kwimg">
occopt</a>
to 4 in telast_6.in, the option used in lesson 4, the C's
show no sign of convergence. At ngkpt=3*16, errors are still
~5%. The reasons that this supposedly superior smoothing function performs so poorly in this context is a future research topic. The main thing to be learned is that checking convergence with respect to all relevant parameters is <b> always</b> the user's responsibility. Simple systems that include the main physical features of a complex system of interest will usually suffice for this testing. Don't get caught publishing a result that another researcher refutes on convergence grounds, and don't blame such a mistake on Abinit!</p>
<p>Finally, we conclude the lesson with a comparison with experiment. Converting
the C's to standard units (Ha/Bohr^3 = 2.94210119E+04 GPa) and using zero-temperature extrapolated experimental results from <cite>P. M.
Sutton, Phys. Rev. 91, 816 (1953)</cite>, we find<br>
</p>
<pre> C_11(GPa) C_12(GPa) C_44(GPa)<br><br>Calculated 123.7 59.9 38.1<br>Experiment (T=0) 123.0 70.8 30.9<br></pre>
Is this good agreement? There isn't much literature
on DFT calculations of full sets of elastic constants. Many calculations
of the bulk modulus (K=(C<small><small>11</small></small>+2C<small><small>
12</small></small> )/3 in the cubic case) typically are within
10% of experiment for the LDA. Running telast_6 with ixc=11, the Perdew-Burke-Enzerhof GGA, increases the calculated C's by 1-2%, and wouldn't be expected to make a large difference for a nearly-free-electron metal.<br>
<p></p>
<p> </p>
<p> </p>
<hr> <p> <br>
This ABINIT tutorial is now finished... </p>
<script type=“text/javascript” src=“list_internal_links.js”> </script>
</body> </html>