Periodic systems III. Equation of state

An alternative to geometry optimization for periodic systems is to fit energy versus cell volume to an equation of state (EOS). The fit yields the equilibrium volume $V_0$, the equilibrium lattice energy $E_0$, and the bulk modulus $B_0$.

Birch–Murnaghan equation of state

One standard form is the Birch–Murnaghan EOS, written as

$$\begin{gathered} E(V) = E_0 + \frac{9V_0 B_0}{16}\Big[\big[(\tfrac{V_0}{V})^{2/3}-1\big]^3 B_0' \\ + \big[(\tfrac{V_0}{V})^{2/3}-1\big]^2\big[6-4(\tfrac{V_0}{V})^{2/3}\big]\Big] \end{gathered}$$

Here $B_0 = -V(\partial P/\partial V)_T$ is the stiffness at equilibrium, and $B_0' = (\partial B_0/\partial P)_T$ is not redundant with it. The bulk modulus itself changes as the crystal is compressed, and $B_0'$ (dimensionless, around 4 for many solids) is the fit parameter that captures that stiffening away from equilibrium.

Calculations scan $\pm 7.5\%$ of diamond’s experimental cell volume in steps of $2.5\%$, with the GTH-cc-pVDZ basis, GTH-HF-rev pseudopotential, and RSDF. The fits at each mesh size are summarized below.

	$a_0$ (Å)	$E_0$ (ha)	$B_0$ (GPa)
$N_k=1^3$	3.923	−10.307	556.383
$N_k=2^3$	3.607	−10.958	689.53
$N_k=3^3$	3.573	−11.022	687.106
$N_k=4^3$	3.568	−11.029	678.893
$N_k=5^3$	3.567	−11.028	676.799
$N_k=6^3$	3.567	−11.027	675.894
$N_k=\infty$	3.567	−11.026	674.652

Fitted EOS with HF energies.

	$a_0$ (Å)	$E_0$ (ha)	$B_0$ (GPa)
$N_k=1^3$	4.004	−10.499	438.281
$N_k=2^3$	3.625	−11.193	590.371
$N_k=3^3$	3.582	−11.277	591.995
$N_k=4^3$	3.575	−11.292	579.976
$N_k=5^3$	3.573	−11.295	578.070
$N_k=6^3$	3.572	−11.296	576.769
$N_k=\infty$	3.572	−11.297	574.986

Fitted EOS with MP2 energies.

For both HF and MP2, $a_0$ and $E_0$ converge sharply after $N_k = 2^3$ and stabilize at the second decimal after $5^3$, showing that the k-mesh is not the most important error source relative to experiment.

The basis set, and the quality of the HF reference

We now compare against a larger basis, GTH-cc-pVQZ, at a $5^3$ mesh.

	$a_0$ (Å)	$E_0$ (ha)	$B_0$ (GPa)
HF	3.557	−11.035761160443705	490.7
MP2	3.555	−11.374095434242212	453.0

Neither $a_0$ nor $B_0$ from the small basis agrees well with these, and the $B_0$ error is large. Curiously, the small-basis MP2 $B_0$ is closest to the reference at the smallest mesh, which smells like cancellation between the calculation’s errors and the fit’s. More telling, MP2 in the small basis is seemingly less accurate than HF for $a_0$, despite being the higher-level method. MP2 inherits the HF orbitals, so a poor HF reference gets amplified rather than corrected, the same pattern I saw in my undergraduate benchmark, where systems that HF fails tend to suffer poor MP2 results too.

Readjusting the MP2 correlation energy onto the larger-basis HF numbers makes the point directly.

	$a_0$ (Å)	$E_0$ (ha)	$B_0$ (GPa)
MP2 @ QZ-HF	3.564	−11.303	393.591

With the better HF reference, $a_0$ moves toward experiment and $E_0$ toward the large-basis value. The $B_0$ from this mixed treatment overshoots downward, consistent with combining a QZ-quality HF curvature with a DZ-quality correlation correction. The two basis errors no longer cancel the way they do in a consistent small-basis fit. The lesson stands either way. The quality of the HF reference matters as much as the correlation treatment on top of it.