CRYSTAL home page CRYSTAL home page

MgO_eos

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: \begin{equation} \label{eq:Bmod} B_0 = -V \left( \frac{\partial P}{\partial V} \right)_T\; . \end{equation} A dimensionless parameter \(B_0^\prime\) can then be defined as its first derivative with respect to the pressure, at constant temperature \(T\): \begin{equation} \label{eq:Bprim} B_0^\prime = \left( \frac{\partial B_0}{\partial P}\right)_T\; . \end{equation} Let us recall that the pressure \(P\) may be written as a function of the volume \(V\) as: \begin{equation} \label{eq:PV} P(V) = - \left( \frac{\partial E}{\partial V}\right)_S\; . \end{equation} 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: \begin{equation} \label{eq:Bmod2} B(V) = V \left( \frac{\partial^2 E}{\partial V^2} \right)_{T,S}\; . \end{equation} 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: \begin{equation} \label{eq:H} H(V) = E(V) + P(V)\times V\; . \end{equation}

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: \begin{equation} \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} \; , \end{equation} 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: \begin{equation} \label{eq:MUR2} P(V) = \frac{B_0}{B_0^\prime} \left[\left(\frac{V_0}{V}\right)^{B_0^\prime}-1 \right] \; . \end{equation}

Birch-Murnaghan EOS

The third-order Birch-Murnaghan isothermal equation of state (based on the Eulerian strain), published in 1947, reads like: \begin{equation} \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 \; . \end{equation} Again, according to equation (\ref{eq:PV}), we can get \(P(V)\) third-order Birch-Murnaghan's EOS: \begin{equation} \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 \; . \end{equation}

Poirier-Tarantola logarithmic EOS

The third-order Poirier-Tarantola logarithmic equation of state (derived from the natural strain), proposed in 1998, is: \begin{equation} \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) \; , \end{equation} while \(P(V)\) Poirier-Tarantola's EOS is: \begin{equation} \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] \; . \end{equation}

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: \begin{equation} \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] \; . \end{equation} 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
 
CRYSTAL home page CRYSTAL home page