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 c₁₁, c₁₂, c₄₄.

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

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 bnds.al
fplot -f plot.plbnds

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

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

rm -f out save.al rst.al
for x in -0.04 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 0.04
do
    rm -f mixm.al wkp.al
    lmf al -vdist=$x -vnk=72 --pr20,20 >>out
done

Preliminaries

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 a³/4 = 110.004, as appropriate for an fcc structure. A typical atom takes up roughly 100 a.u.³; 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 basp.al; 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.

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(E_F), the density-of-states at the Fermi level, is 4.238 / Ry / cell = 0.0408 / Ry / a₀³, using 110 a₀³ 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 h², 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(E_F) = 5.5 / Ry / cell = 0.050 / Ry / a₀³.

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 $$D^\mathrm{free}(E_F) = \frac{2m}{\hbar^2} \frac{k_F}{2\pi^2} \quad\mathrm{with}\quad k_F = (3\pi^2 n)^{1/3}$$

In Atomic Rydberg units ($\hbar = 2m = e^2 = a_0 = 1$), use $n = 3/110 a_0^{-3}$ to get $\quad k_F a_0 = 0.931 \quad \mathrm{and} \quad D^\mathrm{free}(E_F) = 0.0472 \mathrm{Ry/a_0^3}$.

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 $W = (k_F a_0)^2 = 0.866$ 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 $D(E_F)/D^\mathrm{free}(E_F) = 0.050/0.0472 = 1.06$.

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 $v_F$ 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 $$v_F = \left({\partial E}/{\hbar \partial k}\right)_{E_F} = \frac{2\hbar^2 k_F}{2m} = 2 \frac{E_F}{\hbar} \frac{a_0}{k_F a_0}$$

which, for Al, evaluates to 1.86 in atomic units. Expressing v_F 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 $$v_F^\mathrm{free} = 2 \frac{0.866\,\mathrm{Ry}\times 13.6\,\mathrm{eV/Ry}}{6.582\times 10^{-16}\,\mathrm{eV\,sec}} \frac{0.5292\times 10^{-10}\,\mathrm{m}}{0.931} = 2.03 \times 10^6\,\mathrm{m/sec}$$

Since $v^\mathrm{free}(E)$ scales as $\sqrt{E}$, we can readily find $v^\mathrm{free}(E_F/10)$ = 0.59 a.u. = 0.64×10⁶ 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.   ∫ d³k δ(E(k)-E), the quantity   1/2 ∫ d³k δ(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 bnds.al
fplot -f plot.plbnds  ← this makes postscript file fplot.ps

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 E_F, trace the second band from L (which is the first band folded back on itself) and note that it crosses E_F 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 dos.al, 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:

 <|v|> (spin 1) at Ef = 0.4647 a.u. = 0.6153 x 10^6 m/s.  Resolve by band:
  ib   <|v|>       ib   <|v|>       ib   <|v|>
   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 dosq.al. 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×10⁶ m/sec).

Next consider the vector velocity at the Fermi level, E_F=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 dosq.al dosq.al.1
$ lmdos al -vnk=72 --dos~bands=1:3~rdm~window=-.05,.05~npts=101~mode=1~vec=0,1,0  --kres=0.001834
$ cp dosq.al dosq.al.2
$ lmdos al -vnk=72 --dos~bands=1:3~rdm~window=-.05,.05~npts=101~mode=1~vec=0,0,1  --kres=0.001834
$ cp dosq.al dosq.al.3

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

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

Only bands 2 and 3 cross E_F, 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×10⁶ m/s, about 80% of the free electron result (2.03×10⁶ m/s).

4. Electronic contribution to the specific heat

The electronic contribution to the specific heat at constant volume is given by $$C_v^\mathrm{el} = T \left(\frac{\partial S}{\partial T}\right)$$

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 ¹ $$c_v^\mathrm{el} = \frac{\pi^3}{3} k_B^2 T D({E_F}) \rightarrow \frac{\pi^2}{2}\left(\frac{k_BT}{E_F}\right) n k_B$$

which is proportional to D(E_F). 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 Mermin¹, Eq. 2.84) $$C^\mathrm{el,free} = \frac{\pi^2}{2} ZR \frac{k_B}{E_F} T = \gamma T$$

where Z is the number of valence electrons, and $R = 8.314\,\mathrm{joule/mole} = 1.99\,\mathrm{cal/mole}$ the gas constant, and $E_F/k_B$ is the the “Fermi temperature” $T_F$. Expressing $k_B = 6.337{\times}10^{-6}$ Ry/K in atomic Ry units, we obtain $\gamma^\mathrm{free} = 0.90$ mJ-mole$^{-1}$ K$^{-2}$.

The experimental value, $\gamma^\mathrm{exp} \approx 1.35$ mJ-mole$^{-1}$ K$^{-2}$, is larger by a factor 1.4. The LDA predicts $\gamma^\mathrm{LDA}/\gamma^\mathrm{free} = D(E_F)/D^\mathrm{free}(E_F) = 1.06$, and thus it significantly underestimates $\gamma$.

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

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

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 cv.al. You must use Fermi-Dirac sampling.

cv.al 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 $k_B$. It is essentially linear in T, showing that the next order term (proportional to T³) is negligible.

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

$ mcx cv.al  -e1 x4/x1

You should find that $\gamma = 0.000113 k_B$/atom. To convert to (mJ/mole/K²) units substitute $k_B \rightarrow R = N_A k_B$ to obtain $\gamma = 0.000113 \times 8.314 = 0.94$ mJ-mole$^{-1}$ K$^{-2}$.

5. Elastic constants

This section shows how to calculate two of the three independent shear constants in Al. We first calculate c₁₁−c₁₂, 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 save.al rst.al
foreach x ( -0.04 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 0.04)
    rm -f mixm.al wkp.al
    lmf al -vdist=$x -vnk=72 --pr20,20 >>out
end

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
do
...
done

Note that the mixing file eigenvalue weights file (mixm.al and wkp.al) are deleted for each new shear calculation. save.al 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 save.al | ~/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 c₄₄ 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⁻⁷. 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: c₁₁=1.08, c₁₂=0.62, c₄₄=0.28, so that c₁₁−c₁₂ = 0.46. The calculated c₁₁−c₁₂ 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 c₁₁=1.14, c₁₂=0.62, c₄₄=0.32, so that c₁₁−c₁₂ = 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 c₄₄. 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 c₁₁−c₁₂. 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 ².

rm -f out save.al rst.al
foreach x ( -.4 -.3 -0.2 -0.1 0 0.1 0.2 0.3 .4 )
    rm -f mixm.al wkp.al
    lmf al -vda=$x -vnk=72 --rs=0 --pr20,20 -vbigbas=1 >> out
end

$ cat save.al | ~/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 (c₁₁+2c₁₂)/3 = 0.77 erg/cm³ 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/cm³ at the experimental lattice constant, and 0.65 erg/cm³ 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

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

² ctrl.al 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.