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
  2. Electronic contribution to the specific heat
  3. The three elastic constants

Results are compared against free-electron values, which Al resembles.

This tutorial 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


lmchk al
lmfa  al
lmf al > out

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  ← this makes postscript file

For Fermi velocity:

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

For 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
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=24 --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. Its structure and categories are explained in some detail here.

2. Getting Started: Analyzing The Results Of A Band Pass

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

Note: (the name of the input file) can be in place of al when it appears on the command line. You can also use mpirun -n # lmf in place of lmf.

3. Optimizing The Basis Set

  • One can optimize RSMH efficiently 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.

  • 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

(Switch -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=-.2929211. 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 -.2947882. The extra basis will have a small, but non-negligible effect on the elastic constants and specific heat.

Read more about basis sets here.

4. Self-consistent LDA calculation

A fully self-consistent full-potential calculation can be performed by invoking:

$ lmf al > out

Now the 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.

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

    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).

For a more detailed description of the self-consistency cycle, see the PbTe tutorial.

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 output the following data can be seen:

 BZWTS : --- Tetrahedron Integration ---
 BZINTS: Fermi energy:      0.003302;   3.000000 electrons;  D(Ef):    4.238
         Sum occ. bands:   -0.9779714  incl. Bloechl correction:   -0.007365

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     ehf         &Delta;E        Bloechl
16    4.765    0.004289  -0.9776583  -0.290651  -0.000026    -0.004324
24    5.250    0.002993  -0.9775536  -0.290546   0.000079    -0.002026
48    5.487    0.002135  -0.9776274  -0.290620   0.000005    -0.000539
72    5.503    0.001952  -0.9776318  -0.290624   0.000001    -0.000244
96    5.493    0.001877  -0.9776323  -0.290625   0.000000    -0.000138

ΔE is the difference in either the band sum or the Harris-Foulkes energy relative to the nk=96 result ( columns 4 and 5 exactly track each other because the density is kept fixed.) The Bloechl correction (a measure the second-order term in k spacing), is smaller than the actual error (approximately third-order), and it vanishes in the limit nk→∞. Thus these values change little with a still finer k mesh. 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 20225, k=  0.00000  0.00000  0.00000
 -0.8186  0.9378  0.9378  0.9378  1.0823  1.0823  1.0823  1.4662  1.4662

Relative to the band bottom, the Fermi level is then 0.0019–0.8186 = 0.821 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.821/0.866 = 0.948, 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.8186 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.8186+0.821/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.8186+0.821/10 is a switch telling lmdos to evaluate the object (|v| in this case) for each tetrahedron that encloses the energy (-0.8186+0.821/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.4637 a.u. = 0.614 x 10^6 m/s.  Resolve by band:
  ib   &lt;|v|&gt;       ib   &lt;|v|&gt;       ib   &lt;|v|&gt;
   1   0.4637  |    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 dosq.ext. 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.001952. 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.001952
$ cp
$ lmdos al -vnk=72 --dos~bands=1:3~rdm~window=-.05,.05~npts=101~mode=1~vec=0,1,0  --kres=0.001952
$ cp
$ lmdos al -vnk=72 --dos~bands=1:3~rdm~window=-.05,.05~npts=101~mode=1~vec=0,0,1  --kres=0.001952
$ cp

You should find for all three cases a table like this one:

 &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.2143  |

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).

5. 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,1,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.8946e-6  0.0112692
    200.0  0.0012666  2.9519e-5  0.0216945
    300.0  0.0018999  6.5285e-5  0.0337297
    400.0  0.0025332  1.1582e-4  0.0457702
    500.0  0.0031666  1.8113e-4  0.0575095
    600.0  0.0037999  2.6107e-4  0.0690252
    700.0  0.0044332  3.5556e-4  0.0804263
    800.0  0.0050665  4.6452e-4  0.0917801
    900.0  0.0056998  5.8793e-4  0.1031180

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.

6. 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=48 --pr20,20 >>out

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 48 divisions of k-points.

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 -.2901507   -.2900076
  -.03 -.2903599   -.2902814
  -.02 -.2905064   -.2904724
  -.01 -.2905906   -.2905803
     0 -.2906163   -.2906163
   .01 -.2905912   -.2905807
   .02 -.2905093   -.2904743
   .03 -.2903652   -.2902898
   .04 -.2901593   -.2900276

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=-.29061679  f'=  -3.69767e-5  .518998  -.475787  843.889  2653.85  -7633330

The first derivative is very small (it should be exactly 0; the minimum shear from this polynomial fit falls at 7×10−7, which is zero within the numerical resolution of the data). The second derivative, 0.52 Ry, can be converted to standard 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.46 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. Agreement is fortuitously good! This is in part because the lattice constant used here (7.606) is smaller than the experimental one at room temperature (7.652 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 .69 Ry.

Conversion to standard units in this case is:

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

The final shear constant, the bulk modulus, you can calculate by varying the lattice constant. We do not do it here, but 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=48 --rs=0 --pr20,20 -vbigbas=1
$ cat | ~/lmf/utils/vextract c da ehf | mcx . -e2 x1+7.606 x2 > dat

You should get:

    7.206000   -0.284829
    7.306000   -0.288145
    7.406000   -0.290102
    7.506000   -0.290873
    7.606000   -0.290616
    7.706000   -0.289473
    7.806000   -0.287569
    7.906000   -0.285020
    8.006000   -0.281926

Use pfit to determine the minimum-energy lattice constant. You should find that the minimum is at 7.529 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.

$ set a = 7.529
$ 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.78 erg/cm3 at the LDA lattice constant.

7. Other things to read and do

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

The generation of the core-level spectroscopy, Mulliken analysis or density-of-states is done first by invoking lmf with the appropriate command-line switch, followed by lmdos. The lmdos step is illustrated (including a way to plot results) in the ASA tutorial.

To compute the density of states, see this tutorial. To compute core-level EELS spectra or Mulliken analysis in Fe, try running

    ~/lmf/fp/test/test.fp fe 2

To compute total or partial DOS in hcp Co, try running

    ~/lmf/fp/test/test.fp co 2

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.