CRYSTAL Tutorial Project - Phase transitions tutorial

CaO phase transition

Phase Transitions in Crystalline Systems
(Using the MULFAS program)
M. Llunell

Introduction
Theoretical background
Using MULFAS
Volume and energy data
Exercise 1: CaO phase transition
Exercise 2: Multi-phase transition
Bibliography

Introduction

This tutorial is devoted to the study of pressure dependent phase transitions. The use of CRYSTAL is combined with the MULFAS program, which permits the study of phase transitions that depend on pressure, by analysing the enthalpy difference between products and reactants.

It is assumed that you know how to run CRYSTAL, and how to change the system geometry from a model input deck.

Theoretical background

Pressure is an important variable for condensed matter: under pressure, intermolecular and interatomic distances show variations usually larger than those induced by temperature; for instance, structural modifications can be of the order of ten percent at 50 GPa for a molecular solid. Various possible geometries of the system and the associated relative energies can be determined experimentally as a function of pressure. This evolution as a function of pressure can involve one (polymorphism) or many (solid-state reaction driven by pressure) types of systems.

This tutorial presents a general scheme and the related utilities to deal with phase transitions, by means of two examples:

The structural B1 $\rightleftharpoons$ B2 phase transition of CaO which involves two phases with simple structure [1,2].
The decomposition reactions of MgAl₂O₄ spinel into its oxide components [3,4]:
MgAl₂O₄ $\rightleftharpoons$ MgO + Al₂O₃.

We proceed as follows: for each crystal phase, the total energy E is computed at a number of unit-cell volumes; at each volume, lattice (unit-cell edges) and internal (atomic coordinates) parameters that minimize E are determined. An analytical representation of E vs V is obtained by using a polynomial expression or the Murnaghan equation of state (or any other fitting function). The Murnaghan function, by far the most universally adopted [5], is as follows:

$\begin{displaymath} E(V)=B_0V_0\left\lbrack \frac{1}{B'(B'-1)} \left( \frac{V_0}... ...^{B'-1} + \frac{V}{B'V_0} - \frac{1}{B'-1} \right\rbrack + E_0 \end{displaymath}$

(1)

The fitting parameters are V₀ (volume of minimum energy), B₀ (zero-pressure bulk modulus) and B' (pressure derivative of the bulk modulus B at P=0), and E₀ (minimum energy).

From the $P(V) = -\partial E/\partial V$ relationship, we get

$\begin{displaymath} V(P) \; = \; V_0 \; \left( \; \frac{B'}{B_0} P + 1 \; \right)^{-\frac{1}{B'}} \end{displaymath}$

(2)

Inserting equation (2) in (1), one obtains the analytic E vs P dependence; by adding the PV term, the enthalpy as a function of pressure is obtained:

$\begin{displaymath} H(P) = E + PV = E_0 + \frac{B_0V_0}{B'-1} \left\lbrack \left( \frac{B'}{B_0} P + 1 \right)^{1-1/B'} - 1 \right\rbrack. \end{displaymath}$

(3)

At T=0 K, the transition pressure, P_t, corresponds to the point where all the systems have the same enthalpy:

$\begin{displaymath} \Delta H(P_t) = 0 \end{displaymath}$

(4)

This equation is solved numerically yielding the transition or the decomposition pressure. Knowing this pressure, the equation of state followed by the system during the process can be deduced.

A rough estimate of the transition pressure, P_t', can be obtained just from the knowledge of the equilibrium values (E₀ and V₀) for each phase. At T=0K, the enthalpy as a function of pressure can be evaluated as follows:

$\begin{displaymath} H(P) = H_0 + \int_0^P dH = E_0 + \int_0^P VdP \end{displaymath}$

(5)

Using this relation, equation (4) becomes:

$\begin{displaymath} \Delta H(P_t) \; = \; \Delta E_0 \; + \; \int_0^{P_t} \Delta V dP \; = \; 0 \end{displaymath}$

(6)

If the pressure dependence of $\Delta V$ is negligible, one obtains

$\begin{displaymath} P_t \simeq P_t^{'} = -\Delta E_0 / \Delta V_0 \end{displaymath}$

(7)

Using MULFAS

To analyse the phase transition we will use the MULFAS program. The information concerning this program can be found in the MULFAS User's Guide. The only information that must be provided in the input files (one name.dat file for each phase) is the volume and energy data. During the interactive execution the user is asked for other details.

To run MULFAS you should go to the directory where the input files are located and type mulfas. The input files for the proposed exercises have already been generated and you can find them in the subdirectory of your phase_transition directory corresponding to each exercise.

In this tutorial, only the third execution type available in MULFAS (transition analysis between two sets of phases) will be used (see dialog 1 in the User's Guide). In the next dialogs, the program will ask you for further details: the name of the input files, the fitting function to be used, the pressure interval to be explored (see dialogs 2 to 4 in the User's Guide). Additionally, a name for the output files is asked (see dialog 5). It is also possible to take into account the zero point vibrational energy difference associated to the reaction under study, but in this tutorial this option will not be used. So, just type 0 when you will be asked (see dialog 6).

Volume and energy data

When searching for a transition pressure, the final result depends on many variables such as the Hamiltonian used in the energy calculation (HF, LDA, ...), the number of initial volume values considered, their distribution around the equilibrium volume or the function used to fit the E vs V data.

First of all, we need the E(V) function for each phase involved in the transition. A set of energy points (obtained for instance from CRYSTAL calculations) is fitted by using an adequate equation of state (or any other fitting functions, such as a polynomial). For a given volume, the energy should correspond to its minimal value, in other words, a restricted geometry optimization should be performed at each volume.

For sake of simplicity, the volume and energy data for the exercises of this tutorial have already been obtained by using CRYSTAL, however it could be interesting to discuss a bit how to obtain them. We will analyse two different phase transitions:

$CaO (B1) \rightleftharpoons CaO (B2)$
$MgO+Al_{2}O_{3} \rightleftharpoons MgAl_{2}O_{4}$

The first step is to decide which volumes are to be considered for constructing the E(V) curve. The explored volume interval must be such to include the transition pressure, P_t (if it exists). An approximate estimate of P_t can be obtained from the E₀ and V₀ data. Depending on the computational cost of any restricted geometry optimization (V=constant), it is possible to minimize the total energy at a large or small number of points. What is really important, especially when using a polynomial function, is to choose these values around the equilibrium volume, obtained from experimental data, or calculated with a non restricted geometry optimization.

The next step is the geometry optimization under the constant volume restriction for each volume. The global cost of these optimizations depends on both the cost of any energy calculation and the number of steps necessary to find the minimal energy. The first one is directly related to the system size and the second one to the degrees of freedom of the structure.

It is also important to notice that both E and V should refer to the same unit formula, and in the multi-phase transitions it should be chosen according to the reaction stoichiometry. We will return to this point when studying the transition examples included in this tutorial.

Once the volume-energy data for each phase are obtained, they are fitted with MULFAS. There are a few points you must take into account when fitting the E(V) curve:

It is obvious that the number of points must be larger than the number of parameters in the fitting function.
The bulk modulus, B (second derivative of E) is much more dependent on the number of points than E₀ and V₀ (related to the first derivative of E). It is important to remember that pressure, which is the quantity we are interested in, is related to the first derivative of E.
The Murnaghan equation of state is usually more stable than a polynomial fit under modifications of the initial data. However, if a minimal number of well distributed initial volume values are available, a medium order polynomial (4-6th order) can also be used.

As an example, the E(V) function has been determined for MgO, using two different sets of initial volumes (see "Fitting one phase" in the MULFAS User's Guide). The Murnaghan and polynomial results are reported in Table 1.

**Table 1:** Fitting results for different MgO initial data (HF from **CRYSTAL**). The equilibrium values (E₀, V₀, B₀) and the accuracy (sigma or RMS) are tabulated for each case.
Fitting function	Sigma	Volume	Energy	Bulk Modulus
	RMS	(Å³)	(u.a.)	(GPa)
9 points; $\pm$ 2% steps with respect to the equilibrium volume
Murnaghan	0.182^.10^-5	18.3541	-274.68205	183.969
Polyn. 2nd	0.852^.10^-4	18.3562	-274.68215	224.816
Polyn. 3rd	0.476^.10^-5	18.3465	-274.68205	183.000
Polyn. 4th	0.516^.10^-6	18.3534	-274.68205	183.503
Polyn. 5th	0.563^.10^-6	18.3533	-274.68205	183.555
Polyn. 6th	0.505^.10^-6	18.3534	-274.68205	184.223
Polyn. 7th	0.477^.10^-6	18.3541	-274.68205	184.550
5 points; $\pm$ 4% steps with respect to the equilibrium volume
Murnaghan	0.872^.10^-6	18.3545	-274.68205	183.898
Polyn. 2nd	0.967^.10^-4	18.3667	-274.68217	225.159
Polyn. 3rd	0.476^.10^-5	18.3458	-274.68205	183.326
Polyn. 4th	0.000	18.3544	-274.68205	183.568

Exercise 1: CaO phase transition

In this exercise we will study the phase transition of CaO, which presents two polymorphic phases:

$CaO (B1) \rightleftharpoons CaO (B2)$

In both cases, the energy is just depending on the cell parameter a of the cubic cells. No internal coordinates must be optimized.

CaO B1
CRYSTAL
0 0 0
225
"a"
2
20 0. 0. 0.
8 0.5 0.5 0.5
END Title
Dimensionality of the system

Space Group
Cell parameters
Number of non equivalent atoms
Atomic number and cartesian coordinates

CaO B2
CRYSTAL
0 0 0
221
"a"
2
20 0. 0. 0.
8 0.5 0.5 0.5
END Title
Dimensionality of the system

Space Group
Cell parameters
Number of non equivalent atoms
Atomic number and cartesian coordinates

In the EXERCISE_1 directory there are two subdirectories (PHASE_B1 and PHASE_B2) which contain the input CRYSTAL files (Hartree-Fock calculations) for a set of a parameter. Compare the cao01.d12 and cao02.d12 input files. Why has the FIXINDEX option been used in the second one?

After running CRYSTAL calculations, use the information contained in the output files, to complete the MULFAS.dat input files for both phases which have been partially constructed (caoHF_b1.dat and caoHF_b2.dat). Use them to run MULFAS and compute the transition pressure.

There are also other MULFAS input files ready for different hamiltonians:

caoLDA_b1.dat caoBLYP_b1.dat caoB3LYP_b1.dat caoPWGGA_b1.dat caoPBE_b1.dat
caoLDA_b2.dat caoBLYP_b2.dat caoB3LYP_b2.dat caoPWGGA_b2.dat caoPBE_b2.dat

Use them to find the transition pressure in each case. Some results of a complete study of this CaO phase transition are reported in Tables 2-4 and Figures 1-3.


Figure 1: CaO Hartree-Fock and LDA energy as a function of the cell volume. Circle and continuous line and square and dashed line refer to B1 and B2 phases, respectively. Symbols indicate calculated energy points.	Figure 2: Hartree-Fock and LDA enthalpy of the CaO phases. Symbols as in Figure 1.

**Table 2:** CaO equilibrium parameters obtained with various equations of state at the Hartree-Fock level. V₀, B and B' are the equilibrium volume (Å³), the bulk modulus (GPa) and the derivative *dB/dP*, respectively. (Ref. a)
Equation of State	B1			B2
	V₀	B	B'	V₀	B	B'
Murnaghan	28.8	121	3.87	25.4	132	3.85
Birch-Murnaghan 2^nd order	28.8	123	0.	25.4	136	0.
Birch-Murnaghan 3^rd order	28.8	123	3.86	25.4	136	3.79
Lagrangian 3^rd order	28.9	130	3.42	25.6	146	3.05
Davis and Gordon	28.8	124	3.82	25.4	137	3.73
Bardeen	28.8	123	3.85	25.4	135	3.79
Slater	28.8	123	3.85	25.4	137	3.77
Brennan and Stacey	28.8	124	3.84	25.4	137	3.76

a: The program "STATEQ" for obtaining the data of Table 2 has been kindly provided by Prof. Ph. D'Arco.

**Table 3:** CaO equilibrium lattice parameter a₀ (Å), bulk modulus B (GPa) (obtained with the Murnaghan equation of state) and binding energies B.E. (hartree) (evaluated with respect to the atomic references) calculated at the HF and DFT level. The DFT notation corresponds to the exchange-correlation potentials used in the local and non-local approximations: LV: LDA/Vosko-Wilk-Nusair [6, 7]; BL: Becke/Lee-Yang-Parr [8, 9]; B3L: B3LYP [10]; PP: Perdew-Wang/Perdew-Wang [11, 12,13]; PBE: Perdew-Burke-Ernzerhof/Perdew-Burke-Ernzerhof [14].
Hamiltonian	B1			B2
	a₀	B	B.E.	a₀	B	B.E.
HF	4.87	123	0.286	2.94	130	0.243

LV	4.72	139	0.466	2.85	150	0.435

BL	4.86	117	0.388	2.95	123	0.348
B3L	4.84	122	0.383	2.93	130	0.343

PP	4.82	121	0.410	2.92	129	0.373
PB	4.82	120	0.404	2.92	128	0.367
Exp.	4.81^a	111^a	0.404^b	2.91^c	115^c	-
		112^d			130^a
		115^c

^a: Ref. [15]. ^b: deduced from the Janaf Thermodynamical Tables Ref. [16]. ^c: Ref. [17]. ^d: Ref. [18].

**Table 4:** Phase transition parameters obtained at the Hartree-Fock and DFT level. P_t, V_B1 and V_B2 are the transition pressure (GPa) and the transition volumes (Å³ of the B1 and B2 phases, respectively. Symbols as in Table 3.
Hamiltonian	P_t	V_B1	V_B2
HF	69.4	21.2	19.0

LV	57.4	20.4	18.3

BL	75.2	20.5	18.6
B3L	72.5	20.5	18.5

PP	66.1	20.6	18.6
PB	65.6	20.7	18.7
Exp.	60.0^a	-	-
	65.0^b	-	-
	63.0^c	20.7^c	18.7^c

^a: Ref. [15]. ^b: Ref. [18]. ^c: Ref. [17].

Figure 3: Equation of state (volume, Å³, versus pressure, GPa), for the two CaO phases, as obtained at the Hartree-Fock and LDA level. Vertical lines represent the transition pressure.

Exercise 2: Multi-phase transition

In this exercise we will study a solid state reaction which implies more than two compounds:

$MgO+Al_{2}O_{3} \rightleftharpoons MgAl_{2}O_{4}$

In this case to construct the volume-energy list for each phase, more than one energy calculation for each volume is required for MgAl₂O₄ and Al₂O₃, because their geometries depend on more than one parameter (2 and 3, respectively).

MgO
CRYSTAL
0 0 0
225
"a"
2
12 0. 0. 0.
8 0.5 0.5 0.5
END Title
Dimensionality of the system

Space Group
Cell parameters
Number of non equivalent atoms
Atomic number and cartesian coordinates

Al2O3
CRYSTAL
0 0 0
167
"a" "c"
2
13 0. 0. "z_Al"
8 "x_O" 0. 0.25
END Title
Dimensionality of the system

Space Group
Cell parameters
Number of non equivalent atoms
Atomic number and cartesian coordinates

MgAl2O4
CRYSTAL
0 0 0
227
"a"
3
13  0.500   0.500   0.500
12  0.125   0.125   0.125
8  "x_O"   "x_O"   "x_O"
END Title
Dimensionality of the system

Space Group
Cell parameters
Number of non equivalent atoms
Atomic number and cartesian coordinates

The cost of energy minimization for a fixed volume is higher in Al₂O₃, where three geometry parameter must be optimized (a/c, z_Al and x_O). We remind that E and V should refer to the same unit formula and should be chosen according to the reaction stoichiometry. In this example:

for MgO, with a fcc cubic cell whose primitive cell contains just one [MgO] formula (Z=1), volumes and energies must refer to the primitive cell; that is to say, the energy obtained in the CRYSTAL output and a volume equal to a³/4, where a is the lattice parameter of the conventional cubic cell.
for both Al₂O₃ and MgAl₂O₄, with a primitive cell containing two formula units (Z=2), the volumes and energies listed in .dat file must be 1/2 of those extracted from the CRYSTAL output file.

In the EXERCISE_2 directory there are the necessary input data files for running MULFAS using data from both HF and LDA calculations:

mgo_HF.dat al2o3_HF.dat mgal2o4_HF.dat
mgo_LDA.dat al2o3_LDA.dat mgal2o4_LDA.dat

Use these files to compute the transition pressure. Some results for this study are reported in Table 5 and Figures 4-6.

**Table 5:** Volume (Å³) and energy (hartree) differences per formula unit at T = 0, P = 0 between MgAl₂O₄ and the sum of MgO + a-Al₂O₃. Pressure (GPa) of the spinel decomposition reaction MgAl₂O₄ $\rightarrow$ MgO + a-Al₂O₃ at T = 0 (approximate P_t', see equation (7), and exact P_t values). Theoretical (calculated at the Hartree-Fock and LDA levels) and experimental values are given.
	$\Delta V_0$	$\Delta E_0$	P_t'	P_t
HF	5.03	-0.020	17.1	18.7
LDA	5.03	-0.008	6.8	7.0
Exp.^a	4.80	-0.008	7.7	8.0
Exp.^b				13.3

^a: Data from compression experiment (15 GPa at 1800 K, see Ref. [19]).
^b: Previous data extrapolated to 0 K.

Hartree-Fock	LDA



V (3)	V (3)
Figure 4: HF (left) and LDA (right) total energy of MgO (top), Al₂O₃ (middle) and MgAl₂O₄ (bottom) as a function of the cell volume. Circles indicate calculated energy points.


Figure 5: Calculated enthalpies (E+PV) per formula unit (hartree) of spinel (full line) and of the assembly MgO + a-Al₂O₃ (dashed line) vs pressure (GPa) at the Hartree-Fock and LDA levels.	Figure 6: Equation of state, volume (Å³) versus pressure (GPa), of spinel (full line) and of the MgO + a-Al₂O₃ assembly (dashed line) followed during the spinel decomposition into its oxide components determine at the Hartree-Fock and LDA levels.

Bibliography

1: M. Caus�, R. Dovesi, C. Pisani and C. Roetti, Phys. Rev. B 33, 1308 (1986).
2: M.P. Habas, R. Dovesi and A. Lichanot, J. Phys.: Condens. Matter 10, 6897 (1998).
3: M. Catti, G. Valerio, R. Dovesi and M. Caus�, Phys. Rev. B 49, 14179 (1994).
4: M. Catti, F. Freyria Fava, C. Zicovich and R. Dovesi, Phys. Chem. Minerals 26, 389 (1999).
5: F. D. Murnaghan, Proc. Natl. Acad. Sci. USA 30, 244 (1944).
6: P.A.M. Dirac, Proc. Cambridge Phil. Soc. 26, 376 (1930).
7: S.H. Vosko, L. Wilk and M. Nusair, Can. J. Phys. 58, 1200 (1980).
8: A.D. Becke, Phys. Rev. A 38, 3098 (1988).
9: C. Lee, W. Yang and R.G. Parr, Phys. Rev. B 37, 785 (1988).
10: A.D. Becke, J. Chem. Phys. 98, 5648 (1993).
11: J.P. Perdew and Y. Wang, Phys. Rev. B 33, 8800 (1986).
12: J.P. Perdew and Y. Wang, Phys. Rev. B 40, 3399 (1989).
13: J.P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992).
14: J.P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
15: P. Richet, H.K. Mao and P.M. Bell, J. Geophys. Res. 93, 15279 (1986).
16: Janaf Thermodynamical Tables, J. Phys. Chem. Ref. Data Suppl. 1 14 (1985).
17: J.F. Mammone, H.K. Mao and P.M. Bell, Geophys. Rev. Lett. 8, 140 (1981).
18: R. Jeanloz, T.J. Ahrens, H.K. Mao and P.M. Bell, Science 206, 829 (1979).
19: T. Irifune, K. Fujino and E. Ohtani, Nature 349, 409 (1991).