Questaal Home

Elastic Constants and Electronic Contribution to Specific Heat in Al

This tutorial applies the lmf code to elemental Al to calculate, within the local-density approximation, several common materials properties:

  1. Density-of-States at the Fermi energy and the Fermi velocity. Results are compared against free-electron values, which Al resembles.
  2. Electronic contribution to the specific heat. The LDA predicts a value close to the free-electron result, but it deviates significantly from the experimental value.
  3. The three independent elastic constants c11, c12, c44.

In addition, this tutorial illustrates several features of the Questaal package. Also it does not use the automatic input file maker or basis set generator; instead a complete input file is supplied.

Calculations execute very quickly and can easily be performed on a laptop.

Command summary


lmfa  al
lmf al > out.lmfsc

Make energy bands

lmchk al --syml~n=41~lbl=LGX~q=1/2,1/2,1/2,0,0,0,1,0,0,1,1/2,0,0,0,0,3/4,3/4,0
lmf al --band~fn=syml
echo -12,5,5,10 | plbnds -fplot~scl=.7,.6 -ef=0 -scl=13.6 -lbl=L,G,X,W,G,K
fplot -f plot.plbnds

For Fermi velocity:

lmf al -vnk=72 --quit=rho --pr20
lmdos al -vnk=72 --dos~bands=1:3~rdm~window=-.05,.05~npts=101~mode=1~vec=1,0,0  --kres=0.001834

For electronic contribution to the specific heat:

lmf al -vfermi=t -vnk=72 --quit=rho --rs=1,1,0,1 --cvK:100:900:100

For shear modulus:

rm -f out
for x in -0.04 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 0.04
    rm -f
    lmf al -vdist=$x -vnk=72 --pr20,20 >>out


The input file uses some features of the preprocessor (mainly it uses symbolic variables), so you may wish to go through the standard lmf tutorial, first. This tutorial explains in more detail the input file, and the workings of the lmf basis set. See also the PbTe tutorial, which explains the input and output in more detail.

This tutorial uses a number of Questaal executables and scripts, e.g. lmfa, lmf, lmdos, pfit, dval, and vextract. They are assumed to be in your path, and that the source directory is ~/lmf.

1. Building The Input File

The input file for Al is provided here. (It was built by hand, not autogenerated from, e.g. an init file).

Its structure and categories are explained in some detail here.

2. Getting Started: Some Preliminary Checks (Optional)

Use lmchk to verify that:

  • The preprocessor works as advertised:

    lmchk al --showp

    For example, the line containing RSMH= should have been turned into

    RSMH=1.8,1.8,1.8    EH=-.1,-.1,-.1
  • Invoke:

    lmchk al

    and verify from the output that the sphere overlap is about 0.6% (a safe number) and that sum-of-sphere volumes equals about 75% of the total cell volume.

    The lattice constant a=7.606 and unit cell volume are printed near the beginning of the output. Confirm that its value is a3/4 = 110.004, as appropriate for an fcc structure. A typical atom takes up roughly 100 a.u.3; this number is slightly larger than typical because Al is a column III element.

    Use lmfa to generate densities for the free atom. Atomic densities will be overlapped to make a first trial density for the solid (Mattheis construction).

    lmfa al
  • You can efficiently optimize shape of the envelope function controlled by RSMH with:

    lmf al --optbas

    You can take the result of the optimized calculation (saved in basp2.ext) and copy them to your input file. Or, you can have lmf automatically read these parameters from file; tokens in HAM_AUTOBAS control what lmf reads from this file, and what lmfa writes to it.

    In this case the optimizer didn’t improve the basis, so there is nothing to do.

  • We might consider the effect of enlarging the basis, which as we noted we can do from the command line with -vbigbas=1. Do this:

    lmf al -vhf=1
    lmf al -vbigbas=1 -vhf=1

    (-vhf=1 causes the preprocessor to set token OPTIONS_HF to 1, which tells lmf to run in non-self-consistent “Harris Foulkes” mode.)

    In the first line, lmf uses a “minimum basis” case with 9 orbitals; the Harris-Foulkes total energy is ehf=-.2923995. The latter case adds additional envelope functions of s and d types, increasing the basis to 15 orbitals. Thus lowers ehf by about 2 mRy to -.2942652. The extra basis will have a small, but non-negligible effect on the elastic constants and specific heat.

    Read more about basis sets here.

3. Self-consistent LDA calculation

A self-consistent density and potential can generated as follows:

lmfa al
lmf al > out.lmfsc

For a detailed description the steps lmf takes in an iteration to self-consistency, see the PbTe tutorial.

In the first iteration, an output density is generated. With the output density, the HKS energy functional can be evaluated. Also the logarithmic derivative parameters P are floated to the band centers-of-gravity. How the P’s are floated is prescribed by tokens IDMOD.

At self-consistency, the input density that generates the potential, is the same as the output density generated by the potential.

A new density is constructed as follows:

  • The output density is screened using the Lindhard function as an estimate for the static dielectric function , provided the Lindhard parameter ELIND is nonzero:

  • An estimate for the self-consistent density is made by mixing nin and nout* using a mixing scheme described here.

    Because the Lindhard function provides a reasonable description for , lmf should converge rapidly with iteration. You can see this by monitoring the RMS change in charge density as a function of iteration (grep DQ out.lmfsc) You should see that it converges to a very tight tolerance in 6 iterations.

    For another measure of approach to convergence, look at parameter tj printed out at each mixing step. tj indicates what portion of prior iterations are included in the mixing to make a trial density for the next iteration. (This quantity is made only with Andersen mixing.) The number of tj printed is the number of prior iterations included in mixing for that iteration. Ideally, tj should be as close to zero as possible. tj→1 means that lmf much preferred the prior iteration to the current one, and if this happens regularly you should turn down the mixing parameter beta (see Eq (3) on this page). If tj is consistently significantly less than 0, try increasing beta. In this case tj is consistently close to 0, which is a reflection of the fact that Al is a nearly free-electron metal and the Lindhard function is a fairly good description of the screening.

  • The self-consistent density is saved in rst.ext, unless you specify otherwise using –rs=. To see how to use –rs=, see this page.

  • At the end of the iteration, the total energies are printed and checks for convergence against tolerances specified by ITER_CONV and ITER_CONVC. The RMS DQ generated at the mixing step is the measure compared against tolerance CONVC; the change in energy from one iteration to the next is tested against tolerance CONV. Both tolerances must be satisfied, unless you set CONV=0 or CONVC=0 (0 tolerance tells lmf not to test that parameter).

Fermi energy and Density of States at the Fermi level

The physical quantities computed here (specific heat, Fermi velocity, shear modulus) require a fine k mesh to converge. (The density just calculated can be converged with a coarser k mesh).

In the last iteration the following data can be seen in out.lmfsc:

 BZWTS : --- Tetrahedron Integration ---
 BZINTS: Fermi energy:      0.003696;   3.000000 electrons;  D(Ef):    4.478
         Sum occ. bands:   -0.9777570  incl. Bloechl correction:  -0.0071308

D(EF), the density-of-states at the Fermi level, is 4.238 / Ry / cell = 0.0408 / Ry / a03, using 110 a03 for the cell volume.

This value is not well converged in the k integration. The Bloechl correction, -0.007365 Ry here, is a measure of the change in the band sum, including effectively knocking out the second-order error in the standard tetrahedron integration. To see how these quantities converge with k mesh, try

lmf al -vnk=24 --quit=band
lmf al -vnk=48 --quit=band
lmf al -vnk=72 --quit=band
lmf al -vnk=96 --quit=band

You should find the following

nk    D(Ef)      E<sub>F</sub>        band sum    Bloechl      &Delta;E
16    4.906    0.003582  -0.9779456  -0.0043892  0.0000131
24    5.245    0.002829  -0.9779523  -0.0020664  0.0000064
48    5.498    0.002030  -0.9779603  -0.0005447 -0.0000016
72    5.503    0.001834  -0.9779591  -0.0002444 -0.0000004
96    5.504    0.001767  -0.9779586  -0.0001377  0.0000001

ΔE is the difference in either the band sum or the Harris-Foulkes energy relative to the nk=128 result (not shown). The Bloechl correction (an estimate of the error in the band sum that scales in the standard tetrahedron method h2, h being the spacing between k points, i.e. h∝1/nk), greatly improves the convergence in nk. It vanishes in the limit nk→∞, as the Table shows. More precisely it knocks out most of the error second-order in h, making the final error scale approximately third-order. In any case, at convergence D(EF) = 5.5 / Ry / cell = 0.050 / Ry / a03.

Using a shifted k-mesh, or another integration schemes, are alternate ways to control k convergence; see the Brillouin zone integration page.

Al is similar to a free-electron metal, and it is useful for later sections to compare this value and the density-of-states to the free electron gas which is

In Atomic Rydberg units (), use to get .

You can also extract the bottom of the band from the lowest eigenvalue at the Γ point:

 bndfp:  kpt 1 of 413, k=  0.00000  0.00000  0.00000
 -0.8187  0.9376  0.9376  0.9376  1.0821  1.0821  1.0821  1.4660  1.4660

Relative to the band bottom, the Fermi level is then 0.003696–0.8187 = 0.822 Ry, which can be compared to the free-electron gas, whose bandwidth is in atomic Ry units.

In summary, the bandwidth as a fraction of the free electron gas result is 0.822/0.866 = 0.949, and .

Fermi velocity

Another useful quantity is the Fermi velocity. lmdos can generate the BZ-integrated average velocity; it can also resolve band velocities by individual band and also by k. Before showing to how to evaluate for a real material, consider the the free electron gas where the Fermi surface is simple, with spherical symmetry and a constant Fermi velocity. It is

which, for Al, evaluates to 1.86 in atomic units. Expressing vF in the second form makes it convenient to convert to other units since the first and second ratio have units of inverse time and length, respectively. Evaluate as

Since scales as , we can readily find = 0.59 a.u. = 0.64×106 m/sec.

lmdos has a special “ballistic conductivity mode” which can calculate the absolute value of the velocity, or well one the x-, y-, and z- components of it. This mode is activated by the mode modifier in the --dos switch mode=1 tells lmdos to calculate, in place of the DOS, i.e.   ∫ d3k δ(E(k)-E), the quantity   1/2 ∫ d3k δ(E(k)-E) ∇E(k). This is the Landauer formula for the ballistic conductivity. mode=3 tells lmdos to calculate the expectation value of |∇E(k)|, which is what we will use here.

Note: the conductivity modes require that you use tetrahedron integration. The BZ is divided into tetrahedra; any tetrahedron that encompasses a given energy will contribute to the number of states at that energy.

When additionally combined with the --kres switch, lmdos will resolve this whatever quantity is calculated by k and by band for a fixed energy you select. Typically you are interested in the average or k- and band- resolved Fermi velocity.

To see what we should expect, first draw the energy bands for Al along some high-symmetry lines. Make the symmetry lines file, generate and draw the bands:

lmchk al --syml~n=41~lbl=LGX~q=1/2,1/2,1/2,0,0,0,1,0,0,1,1/2,0,0,0,0,3/4,3/4,0
lmf al --band~fn=syml
echo -12,5,5,10 | plbnds -fplot~scl=.7,.6 -ef=0 -scl=13.6 -lbl=L,G,X,W,G,K
fplot -f plot.plbnds  &larr; this makes postscript file

Energy bands of Al

Look along the L-Γ line. L is point (1/2,1/2,1/2) in units of 2π/a, and the quadratic band reaches L well below the Fermi level. To reach EF, trace the second band from L (which is the first band folded back on itself) and note that it crosses EF near k=(1/3,1/3,1/3). In the unfolded zone this would correspond approximately to k=(1/2,1/2,1/2)+(1/6,1/6,1/6)2π/a, or |k|=0.95 a.u., close to the free-electron result (0.93 a.u.)

Consider first an energy just a little above the band bottom. Recalling that the band bottom is -0.8187 Ry, use a “Fermi level” equal to [ bottom + bandwidth/10 ]. Do the following:

lmf al -vnk=72 --quit=rho --pr20
lmdos al -vnk=72 --dos~bands=1:3~rdm~window=-.75,-.65~npts=101~mode=3 --kres=-0.8187+0.822/10

lmdos generates DOS for bands 1,2,3. Tags rdm~window=-.75,-.65~npts=101 affect the energy window and formatting of the output, but that is not of concern here. mode=3 tells lmdos to replace the DOS with the absolute value of the velocity.

--kres=-0.8187+0.822/10 is a switch telling lmdos to evaluate the object (|v| in this case) for each tetrahedron that encloses the energy (-0.8187+0.822/10 Ry), and write information for that tetrahedron to file dosq.ext. On exit lmdos should write to stdout the following:

 &lt;|v|&gt; (spin 1) at Ef = 0.4647 a.u. = 0.6153 x 10^6 m/s.  Resolve by band:
  ib   &lt;|v|&gt;       ib   &lt;|v|&gt;       ib   &lt;|v|&gt;
   1   0.4647  |    2   0.0000  |    3   0.0000  |

This shows the BZ average of the velocity, followed by a table resolving this average by band. Note that the only band crossing this level is the first one, as the band plot shows.

Inspect file Its contents are documented here; essentially it gives the velocity for each band in every tetrahedron where it crossed the the specified energy. In this case, the velocity is almost independent of k, reflecting the fact that the band is very nearly parabolic near the band bottom. Moreover, the velocity is near the free-electron value for that energy (0.64×106 m/sec).

Next consider the vector velocity at the Fermi level, EF=0.001834. To collect all the components of the velocity lmdos must be run three times

$ lmdos al -vnk=72 --dos~bands=1:3~rdm~window=-.05,.05~npts=101~mode=1~vec=1,0,0  --kres=0.001834
$ cp
$ lmdos al -vnk=72 --dos~bands=1:3~rdm~window=-.05,.05~npts=101~mode=1~vec=0,1,0  --kres=0.001834
$ cp
$ lmdos al -vnk=72 --dos~bands=1:3~rdm~window=-.05,.05~npts=101~mode=1~vec=0,0,1  --kres=0.001834
$ cp

You should find for all three cases a table like this one (they should be the same since Al has cubic symmetry)

 &lt;v/2.c&gt; (spin 1) at Ef = 0.3137 a.u. = 0.4154 x 10^6 m/s.  Resolve by band:
  ib   &lt;v/2.c&gt;     ib   &lt;v/2.c&gt;     ib   &lt;v/2.c&gt;
   1   0.0000  |    2   0.3340  |    3   0.2130  |

Only bands 2 and 3 cross EF, and band 3 has a lower average value than band 2, as might be guessed from the band plot.

In this case the x-, y-, or z- component of v/2 is shown, which is the Landauer expression for the ballistic conductivity. The BZ average of this quantity should be independent of component because Al has cubic symmetry. However this is not true at any particular k point. A comparison to the free-electron value can be done in different ways. One way (perhaps not the most physical one) is to take the average value of the scalar v, which you can get by running lmdos in mode 3 at the Fermi level. The result is 1.606×106 m/s, about 80% of the free electron result (2.03×106 m/s).

4. Electronic contribution to the specific heat

The electronic contribution to the specific heat at constant volume is given by

In a one-particle description such as the LDA, there is a contribution to S from the Fermi distribution of occupation numbers. To leading order in temperature the specific heat per atom is 1

which is proportional to D(EF). The second expression applies to the free-electron gas.

Heat capacities are usually quoted in joules or calories per mole per degree. Making the conversion the following results (see Ashcroft Mermin1, Eq. 2.84)

where Z is the number of valence electrons, and the gas constant, and is the the “Fermi temperature” . Expressing Ry/K in atomic Ry units, we obtain mJ-mole K.

The experimental value, mJ-mole K, is larger by a factor 1.4. The LDA predicts , and thus it significantly underestimates .

You can tell lmf to generate the entropy and heat capacity by switching from tetrahedron integration to sampling integration with a Fermi-Dirac weighting. is already set up to do this; see the lines

%ifdef fermi
        TETRA=0 N=-1 W={w}

The middle line is included in the input stream if variable fermi is nonzero. Then the tetrahedron method is turned off, and the code uses sampling integration. To see what N and W signify, see the documentation or simply invoke

$ lmf --input

and you should find the following:

 BZ_N                   opt    i4       1,  1     default = 0
   N>0: Polynomial order for Methfessel-Paxton sampling
   N=0: Conventional Gaussian sampling
   N<0: Broadening by Fermi-Dirac distribution
   To be used in conjunction with W, next
 BZ_W                   opt    r8       1,  1     default = 5e-3
   N>=0: Line broadening for sampling integration
   N<0 : Temperature for Fermi distribution (Ry)

Setting N=−1 uses the Fermi-Dirac distribution for sampling weights, with a temperature W. Now do:

lmf al -vfermi=t -vnk=72 --quit=rho --rs=1,1,0,1 --cvK:100:900:100

The various switches do the following:

  • -vfermi=t   modifies the input file as described above
  • -vnk=72    increases the fineness of the k mesh
  • –quit=rho   tells lmf to stop after making the density
  • –rs=1,1,1,1   tells lmf to use sampling parameters from the ctrl file. It is the fourth digit that is important here; see here for more details.
  • –cvK:100:900:100   is a switch telling lmf to calculate the electronic specific heat and write the result to file You must use Fermi-Dirac sampling. should look similar to the following:

#    T(K)    T(Ry)     S(k_B)   TdS/dT(k_B)
    100.0  6.3331e-4  7.5034e-6  0.0112793
    200.0  0.0012666  2.9350e-5  0.0226774
    300.0  0.0018999  6.5538e-5  0.0339200
    400.0  0.0025332  1.1600e-4  0.0451674
    500.0  0.0031666  1.8076e-4  0.0564593
    600.0  0.0037999  2.5982e-4  0.0677677
    700.0  0.0044332  3.5320e-4  0.0790644
    800.0  0.0050665  4.6088e-4  0.0903358
    900.0  0.0056998  5.8283e-4  0.1015817

Column 4 records the specific heat/atom, in units of . It is essentially linear in T, showing that the next order term (proportional to T 3) is negligible.

The factor is the ratio of column 4 to column 1. Use the mcx calculator or your favorite utility to calculate this ratio

$ mcx  -e1 x4/x1

You should find that /atom. To convert to (mJ/mole/K2) units substitute to obtain mJ-mole K.

5. Elastic constants

This section shows how to calculate two of the three independent shear constants in Al. We first calculate c11c12, which we will do by computing the total energy at different lattice distortions, and fitting the curvature of the total energy. The tetragonal distortion is conveniently generated using the line

    SHEAR=0 0 1 1+dist

which distorts the lattice in a way that conserves volume to all orders (this is useful because it tends to be less error-prone). The direction of distortion is set by the first three parameters; the lattice will be sheared along (001).

The first difficulty is that our specification of the FT mesh using token GMAX may cause the program to switch meshes for as parameter dist changes. This is a bad idea, since we want to resolve very small energy differences. So, the first step is to comment out the line with GMAX=gmax in the input and use instead:

    FTMESH=10 10 10

The second difficulty is that the shear constants in Al are difficult to converge, because they require many k-points, as noted above. The following steps are written in tcsh and compute the self-consistent total energy parametrically as a function of ‘dist’:

rm -f out
foreach x ( -0.04 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 0.04)
    rm -f
    lmf al -vdist=$x -vnk=72 --pr20,20 >>out

In bash, the loop is written a little differently:

for x in -0.04 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 0.04

Note that the mixing file eigenvalue weights file ( and are deleted for each new shear calculation. contains total energies ehf for 9 values of dist using 72 divisions of k-points.

Note also the large number of k points. You could use a smaller number, but this large number turns out to get convergence to within 1% even with the tetrahedron method. The shear constants come about 2% larger with 48 divisions.

Extracting ehf and ehk parametrically as a function of dist is very easy with the vextract tool:

cat | ~/lmf/utils/vextract c dist ehf > dat

The key ‘c’ tells vextract that you want lines beginning only with ‘c’: these lines correspond to the final iteration when self-consistency was reached. You can use any regular expression for the key, and also make vextract extract any quantity associated with a variable in the save file.

dat should contain the following. The first column comes from this shear calculation; the second for c44 described next

         c11-c12       c44
  -.04 -.2901498   -.2900107
  -.03 -.2903602   -.2902803
  -.02 -.2905053   -.2904700
  -.01 -.2905905   -.2905833
    0  -.2906205   -.2906207
   .01 -.2905902   -.2905832
   .02 -.2905051   -.2904729
   .03 -.2903609   -.2902879
   .04 -.2901574   -.2900286

Fit a sixth-order polynomial through these points:

echo 0 | pfit 6 dat

You should find that the fit is very good with an RMS error of 2.9×10−7. and that the second line reads

 x=?  f=-.29061971  f'=  4.12354e-6  .574889  .21958  -49.8803  -7384.62  2026670

The first derivative is very small (it should be exactly 0, and it is zero within the numerical resolution of the data). The second derivative, E’‘=0.575 Ry, can be converted to more conventional units with the following:

         vol = a^3/4 = 110 a.u.
    1Ry/a0^3 = 147e12 erg/cm^3
    c11-c12  = (2/3*147/vol)*E'' (10^12 erg/cm^3)
             = 0.51 x 10^12 erg/cm^3

The experimental elastics are in these units: c11=1.08, c12=0.62, c44=0.28, so that c11c12 = 0.46. The calculated c11c12 is slightly larger, in part because the lattice constant used here (7.606) is smaller than the experimental one at room temperature (7.652 a.u.). (The elastic constants at 0 K are c11=1.14, c12=0.62, c44=0.32, so that c11c12 = 0.52. The excellent agreement would be a little less so if the shear were evaluated at the experimental 0 K lattice constant, 7.621 a.u.)

A trigonal shear will yield c44. To compute it, replace the SHEAR token with one corresponding to [111]:

    SHEAR=1 1 1 1+dist

and repeat the calculation in the same way as for c11c12. You should find that the second derivative is .745 Ry. Conversion to standard units in this case is:

    c44 = (1/3*147/vol)*E'' (10^12 erg/cm^3)
        = 0.33 x 10^12 erg/cm^3

The final shear constant is the bulk modulus. You can calculate by varying the lattice constant. If you do this, you are advised to use token STRUC_DALAT, rather than change ALAT 2.

rm -f out
foreach x ( -.4 -.3 -0.2 -0.1 0 0.1 0.2 0.3 .4 )
    rm -f
    lmf al -vda=$x -vnk=72 --rs=0 --pr20,20 -vbigbas=1 >> out
$ cat | ~/lmf/utils/vextract c da ehf | mcx . -e2 x1+7.606 x2 > dat

You should get:

    7.206000   -0.286449
    7.306000   -0.289838
    7.406000   -0.291876
    7.506000   -0.292739
    7.606000   -0.292583
    7.706000   -0.291548
    7.806000   -0.289761
    7.906000   -0.287334
    8.006000   -0.284370

Use pfit to determine the minimum-energy lattice constant. You should find that the minimum is at 7.539 a.u, which is about 1% smaller than the experimental lattice constant at 0K.

$ echo '7.6\nmin' | pfit 6 dat

The experimental bulk modulus B is (c11+2c12)/3 = 0.77 erg/cm3 at room temperature. There is an ambiguity whether to calculate it from the LDA using the experimental lattice constant or the LDA minimum energy one. An expression for B that is valid for fcc at any lattice constant is

B = 4 (147/a) x [ E''(a) - E'(a)/a ]  (10^12 erg/cm^3)

Use pfit to extract the first and second derivatives at the experimental lattice constant, and evaluate B at the experimental lattice constant (7.621 a.u.).

set a = 7.539
set ep = `echo $a | pfit 6 dat  | grep f= | awk '{print $4}'`
set epp = `echo $a | pfit 6 dat  | grep f= | awk '{print $5}'`
dval -a "B($a)=%1,3;3g ep=%1,3;3g epp=%1,3;3g"  "4*147/$a*($epp-$ep*2/$a)" 0+$ep 0+$epp

You should find B = 0.73 erg/cm3 at the experimental lattice constant, and 0.65 erg/cm3 at the LDA lattice constant.

6. Addtional Exercises

1) Try repeating these calculations with the larger basis set (-vbigbas=1)

2) Try repeating these calculations with the PBE functional (XCFUN=0,101,130)

3) Try computing some other properties of the Fermi surface, e.g. its area.

Footnotes and references

1 N. W. Ashcroft and D. Mermin, Solid State Physics (Brooks/Cole, 1976)

2 includes a token DALAT (which evaluates 0 if not specified). The actual lattice constant is ALAT+DALAT. The total lattice constant is separated because because some parameters scale with ALAT but are kept frozen when DALAT varies.