# 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:

- Density-of-States at the Fermi energy and the Fermi velocity. Results are compared against free-electron values, which Al resembles.
- 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.
- The three independent elastic constants
*c*_{11},*c*_{12},*c*_{44}.

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

Self-consistency:

```
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 bnds.al
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 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*^{3}/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*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.

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 n

_{in}and n_{out}^{*}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*(*E _{F}*), the density-of-states at the Fermi level, is 4.238 / Ry / cell = 0.0408 / Ry / a

_{0}

^{3}, using 110 a

_{0}

^{3}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 Δ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 ^{2}*,

*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*) = 5.5 / Ry / cell = 0.050 / Ry / a

_{F}_{0}

^{3}.

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 *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

Since scales as , we can readily find = 0.59 a.u. = 0.64×10^{6} 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*^{3}**k** *δ*(*E*(**k**)-*E*), the quantity 1/2 ∫ *d*^{3}**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*near

_{F}**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^{6} 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^{6} m/s, about 80% of the free electron result (2.03×10^{6} 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*(*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** ^{1}**, 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. *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 . 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 cv.al -e1 x4/x1
```

You should find that /atom. To convert to (mJ/mole/K^{2}) 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 *c*_{11}−*c*_{12}, 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*_{44} 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: *c*_{11}=1.08, *c*_{12}=0.62, *c*_{44}=0.28, so that *c*_{11}−*c*_{12} = 0.46. The calculated *c*_{11}−*c*_{12} 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*_{11}=1.14, *c*_{12}=0.62, *c*_{44}=0.32, so that *c*_{11}−*c*_{12} = 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*_{44}. 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*_{11}−*c*_{12}. 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 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*_{11}+2*c*_{12})/3 = 0.77 erg/cm^{3} 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^{3} at the experimental lattice constant, and 0.65 erg/cm^{3} 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} *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.