Equation of State A. Mahmoud and A. Erba

 Introduction

The CRYSTAL program can perform an automated scan over the volume in order to compute energy $$E$$ $$versus$$ volume $$V$$ curves that are then fitted to various universal Equations of State (EOS) as Murnaghan's, third-order isothermal Birch-Murnaghan's, logarithmic Poirier-Tarantola's and exponential Vinet's, in order to compute equilibrium properties such as bulk modulus $$B_0$$, its first derivative with respect to the pressure $$B_0^\prime$$ and volume/pressure dependence of the energy, enthalpy and bulk modulus. For each volume, a full $$V$$-constrained geometry optimization is performed. This option is activated by inserting the keyword $${\bf{EOS}}$$ at the end of the Geometry input block, which opens a sub-block that must be closed with a keyword $${\bf{END}}$$:

 MgO BULK CRYSTAL 0 0 0 225 4.26068985 2 12 0. 0. 0. 8 0.5 0.5 0.5 Title Dimensionality of the system Space Group Cell parameter Number of non equivalent atoms Atomic number and Cartesian coordinates of the symmetry-irreducible atoms EOS [Optional sub-keywords] END Keyword to perform an EOS calculation End of the EOS input block END End of the geometry input block

A volume range and a number of volumes $$N_V$$ can be defined in input (default values are used otherwise) with the optional $${\bf{RANGE}}$$ sub-keyword. The initial geometry is assumed to be fully optimized; if not, the optional $${\bf{PREOPTGEOM}}$$ sub-keyword must be used to perform a preliminary optimization of both lattice parameters and atomic positions. For each considered volume $$V_i$$, a $$V_i$$-constrained optimization is performed and the corresponding minimum energy $$E_i$$ determined. The set of $$N_V$$ data points $$\left\lbrace V_i, E_i \right\rbrace$$ is fitted to the various EOSs implemented. The fitted energy, enthalpy and bulk modulus are printed at the end of the calculation, for each EOS, as a function of volume and pressure at various points. These volume/pressure ranges can be defined with optional sub-keywords $${\bf{VRANGE}}$$ and $${\bf{PRANGE}}$$, respectively. Typically, this kind of final information analysis can be performed with a complete restart of the calculation, using the $${\bf{RESTART2}}$$ optional sub-keyword. A partial restart from a previous incomplete run can be activated with the optional sub-keyword $${\bf{RESTART}}$$. These two restart options require the external restart file EOSINFO.DAT to be provided.

 A Few Theoretical Remarks

The Equation of State is a pressure-volume or energy-volume relation describing the behavior of a solid under compression or expansion. The simplest isothermal EOS for a solid is the bulk modulus which is a measure of the ability of a material in resisting to changes in volume under uniform compression or expansion. The equilibrium bulk modulus $$B_0$$ of a crystal can be defined as follows: $$\label{eq:Bmod} B_0 = -V \left( \frac{\partial P}{\partial V} \right)_T\; .$$ A dimensionless parameter $$B_0^\prime$$ can then be defined as its first derivative with respect to the pressure, at constant temperature $$T$$: $$\label{eq:Bprim} B_0^\prime = \left( \frac{\partial B_0}{\partial P}\right)_T\; .$$ Let us recall that the pressure $$P$$ may be written as a function of the volume $$V$$ as: $$\label{eq:PV} P(V) = - \left( \frac{\partial E}{\partial V}\right)_S\; .$$ According to equation (\ref{eq:PV}) we can redefine the bulk modulus in equation (\ref{eq:Bmod}) as the second energy derivative with respect to the volume: $$\label{eq:Bmod2} B(V) = V \left( \frac{\partial^2 E}{\partial V^2} \right)_{T,S}\; .$$ We can now define the enthalpy $$H$$ (coinciding with Gibbs' free energy $$G$$ at $$T = 0$$ $$K$$, the standard temperature of ab initio calculations) as a function of the volume $$V$$ simply as: $$\label{eq:H} H(V) = E(V) + P(V)\times V\; .$$

## Murnaghan EOS

In 1944, Murnaghan proposed his famous equation of state which is based on the assumption that the bulk modulus varies linearly with pressure: $$\label{eq:MUR} E(V) = E_0 + \frac{B_0V}{B_0^\prime} \left[\left(\frac{V_0}{V}\right)^{B_0^\prime}\frac{1}{B_0^\prime -1} +1 \right] - \frac{B_0 V_0}{B_0^\prime -1} \; ,$$ where $$V_0$$ and $$E_0$$ are the equilibrium volume and energy, at zero pressure. Application of equation (\ref{eq:PV}) to equation (\ref{eq:MUR}), gives $$P(V)$$ Murnaghan's EOS: $$\label{eq:MUR2} P(V) = \frac{B_0}{B_0^\prime} \left[\left(\frac{V_0}{V}\right)^{B_0^\prime}-1 \right] \; .$$

## Birch-Murnaghan EOS

The third-order Birch-Murnaghan isothermal equation of state (based on the Eulerian strain), published in 1947, reads like: $$\label{eq:BIRMUR} E(V) = E_0 + \frac{9V_0B_0}{16}\left\lbrace \left[ \left( \frac{V_0}{V} \right)^{\frac{2}{3}} -1 \right]^3 B_0^\prime + \left[ \left( \frac{V_0}{V} \right)^{\frac{2}{3}} -1 \right]^2 \left[ 6 -4 \left( \frac{V_0}{V} \right)^{\frac{2}{3}} \right] \right\rbrace \; .$$ Again, according to equation (\ref{eq:PV}), we can get $$P(V)$$ third-order Birch-Murnaghan's EOS: $$\label{eq:BIRMUR2} P(V) = \frac{3B_0}{2} \left[ \left( \frac{V_0}{V} \right)^{\frac{7}{3}} - \left( \frac{V_0}{V} \right)^{\frac{5}{3}} \right]\left\lbrace 1+ \frac{3}{4} (B_0^\prime -4) \left[ \left( \frac{V_0}{V} \right)^{\frac{2}{3}} -1 \right] \right\rbrace \; .$$

## Poirier-Tarantola logarithmic EOS

The third-order Poirier-Tarantola logarithmic equation of state (derived from the natural strain), proposed in 1998, is: $$\label{eq:LOGE} E(V) = E_0 + \frac {B_0 V_0}{2}\left[ \ln \left( \frac{V_0}{V}\right) \right]^{2} + \frac{B_0 V_0}{6} \left[ \ln \left(\frac{V_0}V{} \right) \right]^3 \left(B_0^\prime - 2 \right) \; ,$$ while $$P(V)$$ Poirier-Tarantola's EOS is: $$\label{eq:LOGP} P(V) = B_0 \frac{V_0}{V} \left[\ln \left( \frac{V_0}{V} \right) + \frac{\left( B_0^\prime -2 \right) }{2} \left[ \ln \left( \frac{V_0}{V} \right) \right]^2\right] \; .$$

## Vinet exponential EOS

The exponential Vinet's equation of state, published in 1987, reads: \begin{eqnarray} \label{eq:VinE} E(V)&=&E_0 + \frac{2B_0V_0}{\left(B_0^\prime-1\right)^2}\left\lbrace 2 -\left[ 5 + 3\left( \frac{V}{V_0}\right)^\frac{1}{3} (B_0^\prime -1) -3B_0^\prime \right] \times \right. \nonumber \\ &\times &\left. \exp \left[ -\frac{3}{2} \left( B_0^\prime-1\right)\left[\left( \frac{V}{V_0}\right)^\frac{1}{3} -1\right]\right]\right\rbrace \end{eqnarray} According to equation (\ref{eq:PV}), we get $$P(V)$$ Vinet's EOS: $$\label{eq:VinPV} P(V)=3B_0\left( \frac{V}{V_0}\right)^{-\frac{2}{3}}\left[ 1-\left(\frac{V}{V_0} \right)^\frac{1}{3}\right]\exp \left[- \frac{3}{2}\left( B_0^\prime -1 \right)\left[\left(\frac{V}{V_0}\right)^\frac{1}{3}-1\right]\right] \; .$$ This EOS is found to be quite accurate especially at large compressions and to be exact for harmonic crystals.

 A Guided Tour of Input and Output files

Exercise: Compute the EOS of the MgO crystal. To do so, copy the input file MgO_EOS.d12 in your working directory. The initial part of the input file reads as follows:

 MgO BULK CRYSTAL 0 0 0 225 4.26068985 2 12 0. 0. 0. 8 0.5 0.5 0.5 Title Dimensionality of the system Space Group Cell parameter Number of non equivalent atoms Atomic number and Cartesian coordinates of the symmetry-irreducible atoms EOS RANGE 0.92 1.08 6 PRANGE 15 20 15 END Keyword to perform an EOS calculation Keyword to specify the range of volumes and number of points in the E(V) curve Minimum and maximum variation of the initial volume; number of points in the selected range Keyword to define the pressure range and the number of points, for fitted values Minimum and maximum pressure (in GPa); number of points in the selected range End of the EOS input block END End of the geometry input block

The explored volume interval is specified by providing with the minimum (compression) and maximum (expansion) variation of the equilibrium volume. The set of volumes is then defined according to the number of points in the selected range. For instance, in this case, 6 points have been chosen between 0.92$$\times V_{eq}$$ and 1.08$$\times V_{eq}$$ , where $$V_{eq}$$ is the volume of the equilibrium geometry given as input (assumed to be the fully optimized structure) or as obtained after a preliminary geometry optimization ($${\bf PREOPTGEOM}$$).

At the end of the calculation, 6+1 (where +1 refers to the equilibrium values $$V_{eq}$$ and $$E_{eq}$$) volume/energy data points, as obtained after the geometry optimizations are sorted and printed as follows:

  VOLUME (A^3) ENERGY (a.u.) 17.789658 -2.753023029137E+02 18.382125 -2.753039525749E+02 18.987603 -2.753047829172E+02 19.336585 -2.753049167772E+02 19.606232 -2.753048617406E+02 20.238155 -2.753042539881E+02 20.883512 -2.753030141001E+02 

The following table is then reported with the fitted values of the minimum volume $$V_0$$, energy $$E_0$$, bulk modulus $$B_0$$ and its first derivative $$B_0^\prime$$ with different Equations of State:

  +++++++ FITTING USING ALL POINTS +++++++ EQUATION OF STATE VOL(A^3) E(AU) BM(GPa) BM PRIME ------------------------------------------------------------------------------- MURNAGHAN 1944 19.3620 -275.30491639 156.02 3.70 BIRCH-MURNAGHAN 1947 19.3618 -275.30491678 156.34 3.72 POIRIER-TARANTOLA 1998 19.3617 -275.30491699 156.51 3.73 VINET 1987 19.3617 -275.30491690 156.43 3.73 POLYNOMIAL FITTING VOL(A^3) E(AU) BM(GPa) ------------------------------------------------------------------------------- THIRD ORDER POLYNOMIAL 19.3610 -275.30491935 158.40 FOURTH ORDER POLYNOMIAL 19.3616 -275.30491664 156.23 FIFTH ORDER POLYNOMIAL 19.3612 -275.30491662 156.22 

Additionally, for each EOS (not for the polynomial ones), the following fitted data are reported in the pressure range (from 0 GPa to 50 GPa, as defined in the input with the $${\bf PRANGE}$$ sub-keyword). We report here only the data obtained with Vinet's EOS:

  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ THERMODYNAMIC FUNCTIONS OBTAINED WITH EOS: VINET 1987 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ V = VOLUME, P = PRESSURE, E = ENERGY, G = GIBBS FREE ENERGY, B = BULK MODULUS V (A^3) P (GPa) E (a.u.) G (a.u.) B (GPa) 19.33 0.00 -275.30507192 -275.30506973 156.32 18.91 3.57 -275.30490304 -275.28940772 169.47 18.53 7.14 -275.30443814 -275.27407294 182.28 18.19 10.71 -275.30372779 -275.25903546 194.80 17.87 14.29 -275.30281074 -275.24427100 207.07 17.57 17.86 -275.30171739 -275.22975855 219.11 17.29 21.43 -275.30047224 -275.21548243 230.95 17.03 25.00 -275.29909483 -275.20142431 242.61 16.79 28.57 -275.29760166 -275.18757278 254.10 16.56 32.14 -275.29600594 -275.17391229 265.44 16.34 35.71 -275.29431956 -275.16043582 276.65 16.14 39.29 -275.29255227 -275.14713352 287.72 15.94 42.86 -275.29071215 -275.13399447 298.67 15.76 46.43 -275.28880664 -275.12101213 309.51 15.58 50.00 -275.28684204 -275.10817940 320.24 

Exercise: Calculate the thermodynamic properties of MgO in the volume range 15-20 Å$$^3$$using the keywords $${\bf RESTART2}$$ and $${\bf VRANGE}$$. You can find the corresponding input here: MgO_EOS_vrange.d12. The restart option requires the external restart file MgO_EOS_vrange.eosinfo to be provided. The initial part of the input file reads as follows:

 MgO BULK CRYSTAL 0 0 0 225 4.26068985 2 12 0. 0. 0. 8 0.5 0.5 0.5 Title Dimensionality of the system Space Group Cell parameter Number of non equivalent atoms Atomic number and Cartesian coordinates of the symmetry-irreducible atoms EOS RESTART2 RANGE 0.92 1.08 6 VRANGE 15 20 15 END Keyword to perform an EOS calculation Keyword for complete restart of the calculation Keyword to specify the range of volumes and number of points in the E(V) curve Minimum and maximum variation of the initial volume; number of points in the selected range Keyword to define the volume range and the number of points, for fitted values Minimum and maximum volume; number of points in the selected range End of the EOS input block END End of the geometry input block

The fitted energy, Gibbs free energy (coinciding with the enthalpy at $$T = 0$$ $$K$$ ) and bulk modulus are printed at the end of the output in the requested volume range (the corresponding pressure is also reported):

  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ THERMODYNAMIC FUNCTIONS OBTAINED WITH EOS: VINET 1987 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ V = VOLUME, P = PRESSURE, E = ENERGY, G = GIBBS FREE ENERGY, B = BULK MODULUS V (A^3) P (GPa) E (a.u.) G (a.u.) B (GPa) 15.00 63.37 -275.27885241 -275.06081288 359.56 15.36 55.20 -275.28370375 -275.08927669 335.66 15.71 47.74 -275.28791513 -275.11585247 313.44 16.07 40.93 -275.29154250 -275.14067243 292.77 16.43 34.71 -275.29463652 -275.16385781 273.52 16.79 29.02 -275.29724311 -275.18551991 255.56 17.14 23.82 -275.29940394 -275.20576109 238.81 17.50 19.05 -275.30115694 -275.22467560 223.16 17.86 14.69 -275.30253660 -275.24235034 208.53 18.21 10.70 -275.30357440 -275.25886558 194.84 18.57 7.04 -275.30429903 -275.27429548 182.02 18.93 3.69 -275.30473673 -275.28870871 170.00 19.29 0.62 -275.30491151 -275.30216889 158.74 19.64 -2.20 -275.30484532 -275.31473503 148.17 20.00 -4.77 -275.30455831 -275.32646194 138.24 

 Directory Input Output Description MgO MgO_EOS.d12 MgO_EOS.out EOS of MgO with a defined pressure range MgO_EOS_vrange.d12 MgO_EOS_vrange.out Restart from the previous EOS calculation with a defined volume range