Effect of pressure on phase stability#
So far we have only considered relative stability at a pressure of 0 Pa. We now consider the relative stability of two phases under pressure. We will consider TiO\(_2\) in the rutile and anatase phases.
The pressure is defined by: \(P = -\left(\frac{\partial E}{\partial V}\right)_T\). So if we have an equation of state \(E(V)\) we can calculate the pressure at any volume, or alternatively, given a pressure, compute the volume. Pressure can affect the energy of two phases differently, so that one may become stable under pressure. The condition where a phase transition occurs is when the pressure in the two phases is the same, which occurs at a common tangent.
To show this, we need \(E_{rutile}(V)\) and \(E_{anatase}(V)\).
# run the rutile calculations
from vasp import Vasp
from ase import Atom, Atoms
import numpy as np
B = 'Ti'; X = 'O'; a = 4.59; c = 2.958; u = 0.305;
'''
create a rutile structure from the lattice vectors at
http://cst-www.nrl.navy.mil/lattice/struk/c4.html
spacegroup: 136 P4_2/mnm
'''
a1 = a * np.array([1.0, 0.0, 0.0])
a2 = a * np.array([0.0, 1.0, 0.0])
a3 = c * np.array([0.0, 0.0, 1.0])
atoms = Atoms([Atom(B, [0., 0., 0.]),
Atom(B, 0.5*a1 + 0.5*a2 + 0.5*a3),
Atom(X, u*a1 + u*a2),
Atom(X, -u*a1 - u*a2),
Atom(X, (0.5+u)*a1 + (0.5-u)*a2 + 0.5*a3),
Atom(X, (0.5-u)*a1 + (0.5+u)*a2 + 0.5*a3)],
cell=[a1, a2, a3])
nTiO2 = len(atoms) / 3.
v0 = atoms.get_volume()
cell0 = atoms.get_cell()
volumes = [28., 30., 32., 34., 36.] #vol of one TiO2
for v in volumes:
atoms.set_cell(cell0 * ((nTiO2 * v / v0)**(1. / 3.)),
scale_atoms=True)
calc = Vasp(label='bulk/TiO2/rutile/rutile-{0}'.format(v),
encut=350,
kpts=(6, 6, 6),
xc='PBE',
ismear=0,
sigma=0.001,
isif=2,
ibrion=2,
nsw=20,
atoms=atoms)
calc.calculate()
---------------------------------------------------------------------------
VaspNotConverged Traceback (most recent call last)
Cell In[1], line 45
32 atoms.set_cell(cell0 * ((nTiO2 * v / v0)**(1. / 3.)),
33 scale_atoms=True)
35 calc = Vasp(label='bulk/TiO2/rutile/rutile-{0}'.format(v),
36 encut=350,
37 kpts=[6, 6, 6],
(...)
43 nsw=20,
44 atoms=atoms)
---> 45 calc.calculate()
File ~/vasp/vasp/calculator.py:332, in Vasp.calculate(self, atoms, properties, system_changes)
329 raise VaspRunning(message=status.message or "", jobid=status.jobid)
331 if status.state == JobState.FAILED and not self.force:
--> 332 raise VaspNotConverged(status.message or "Calculation failed")
334 # Not started (or force=True with failed) - need to run
335 Calculator.calculate(self, atoms, properties, system_changes)
VaspNotConverged: OUTCAR incomplete and no running process
# run the anatase calculations
import numpy as np
from vasp import Vasp
from ase import Atom, Atoms
# http://cst-www.nrl.navy.mil/lattice/struk/c5.html
B = 'Ti'; X = 'O'; a = 3.7842; c = 2*4.7573; z = 0.0831;
a1 = a * np.array([1.0, 0.0, 0.0])
a2 = a * np.array([0.0, 1.0, 0.0])
a3 = np.array([0.5 * a, 0.5 * a, 0.5 * c])
atoms = Atoms([Atom(B, -0.125 * a1 + 0.625 * a2 + 0.25 * a3),
Atom(B, 0.125 * a1 + 0.375 * a2 + 0.75 * a3),
Atom(X, -z*a1 + (0.25-z)*a2 + 2.*z*a3),
Atom(X, -(0.25+z)*a1 + (0.5-z)*a2 + (0.5+2*z)*a3),
Atom(X, z*a1 - (0.25 - z)*a2 + (1-2*z)*a3),
Atom(X, (0.25 + z)*a1 + (0.5 + z)*a2 + (0.5-2*z)*a3)],
cell=[a1,a2,a3])
nTiO2 = len(atoms) / 3.
v0 = atoms.get_volume()
cell0 = atoms.get_cell()
volumes = [30., 33., 35., 37., 39.] #vol of one TiO2
for v in volumes:
atoms.set_cell(cell0 * ((nTiO2*v/v0)**(1./3.)),
scale_atoms=True)
calc = Vasp(label='bulk/TiO2/anatase/anatase-{0}'.format(v),
encut=350,
kpts=(6, 6, 6),
xc='PBE',
ismear=0,
sigma=0.001,
isif=2,
ibrion=2,
nsw=20,
atoms=atoms)
calc.calculate()
Now we will fit cubic polynomials to the data.
# fit cubic polynomials to E(V) for rutile and anatase
from vasp import Vasp
import matplotlib.pyplot as plt
import numpy as np
np.set_printoptions(precision=2)
# anatase equation of state
volumes = [30., 33., 35., 37., 39.] # vol of one TiO2 formula unit
a_volumes, a_energies = [], []
for v in volumes:
calc = Vasp(label='bulk/TiO2/anatase/anatase-{0}'.format(v))
atoms = calc.load_atoms()
nTiO2 = len(atoms) / 3.0
a_volumes.append(atoms.get_volume() / nTiO2)
a_energies.append(atoms.get_potential_energy() / nTiO2)
# rutile equation of state
volumes = [28., 30., 32., 34., 36.] # vol of one TiO2
r_volumes, r_energies = [], []
for v in volumes:
calc = Vasp(label='bulk/TiO2/rutile/rutile-{0}'.format(v))
atoms = calc.load_atoms()
nTiO2 = len(atoms) / 3.0
r_volumes.append(atoms.get_volume() / nTiO2)
r_energies.append(atoms.get_potential_energy() / nTiO2)
# cubic polynomial fit to equation of state E(V) = pars*[V^3 V^2 V^1 V^0]
apars = np.polyfit(a_volumes, a_energies, 3)
rpars = np.polyfit(r_volumes, r_energies, 3)
print('E_anatase(V) = {0:1.2f}*V^3 + {1:1.2f}*V^2 + {2:1.2f}*V + {3:1.2f}'.format(*apars))
print('E_rutile(V) = {0:1.2f}*V^3 + {1:1.2f}*V^2 + {2:1.2f}*V + {3:1.2f}'.format(*rpars))
print('anatase epars: {0!r}'.format(apars))
print('rutile epars: {0!r}'.format(rpars))
# get pressure parameters P(V) = -dE/dV
dapars = -np.polyder(apars)
drpars = -np.polyder(rpars)
print('anatase ppars: {0!r}'.format(dapars))
print('rutile ppars: {0!r}'.format(drpars))
print()
print('P_anatase(V) = {0:1.2f}*V^2 + {1:1.2f}*V + {2:1.2f}'.format(*dapars))
print('P_rutile(V) = {0:1.2f}*V^2 + {1:1.2f}*V + {2:1.2f}'.format(*drpars))
vfit = np.linspace(28, 40)
# plot the equations of state
plt.plot(a_volumes, a_energies,'bo ', label='Anatase')
plt.plot(vfit, np.polyval(apars, vfit), 'b-')
plt.plot(r_volumes, r_energies,'gs ', label='Rutile')
plt.plot(vfit, np.polyval(rpars, vfit), 'g-')
plt.xlabel('Volume ($\AA^3$/f.u.)')
plt.ylabel('Total energy (eV/f.u.)')
plt.legend()
plt.xlim([25, 40])
plt.ylim([-27, -26])
To find the conditions where a phase transition occurs, we have to find the common tangent line between the rutile and anatase phases. In other words we have to solve these two equations:
\((E_{anatase}(V1) - E_{rutile}(V2))/(V1-V2) = P_{anatase}(V1)\)
\((E_{anatase}(V1) - E_{rutile}(V2))/(V1-V2) = P_{rutile}(V2)\)
This is a nonlinear algebra problem. We use the func:scipy.optimize.fsolve to solve this problem.
from ase.units import GPa
from numpy import array, linspace, polyval
# copied from polynomial fit above
anatase_epars = array([-1.06049246e-03, 1.30279404e-01, -5.23520055e+00,
4.25202869e+01])
rutile_epars = array([-1.24680208e-03, 1.42966536e-01, -5.33239733e+00,
3.85903670e+01])
# polynomial fits for pressures
anatase_ppars = array([3.18147737e-03, -2.60558808e-01, 5.23520055e+00])
rutile_ppars = array([3.74040625e-03, -2.85933071e-01, 5.33239733e+00])
def func(V):
V1 = V[0] # rutile volume
V2 = V[1] # anatase volume
E_rutile = polyval(rutile_epars, V1)
E_anatase = polyval(anatase_epars, V2)
P_rutile = polyval(rutile_ppars, V1)
P_anatase = polyval(anatase_ppars, V2)
return [(E_anatase - E_rutile) / (V1 - V2) - P_anatase,
(E_anatase - E_rutile) / (V1 - V2) - P_rutile]
from scipy.optimize import fsolve
x0 = fsolve(func, [28, 34])
print('The solutions are at V = {0}'.format(x0))
print('Anatase pressure: {0} GPa'.format(polyval(anatase_ppars, x0[1]) / GPa))
print('Rutile pressure: {0} GPa'.format(polyval(rutile_ppars, x0[0]) / GPa))
# illustrate the common tangent
import matplotlib.pyplot as plt
vfit = linspace(28, 40)
plt.plot(vfit, polyval(anatase_epars,vfit), label='anatase')
plt.plot(vfit, polyval(rutile_epars,vfit), label='rutile')
plt.plot(x0, [polyval(rutile_epars,x0[0]),
polyval(anatase_epars,x0[1])], 'ko-', label='common tangent')
plt.legend()
plt.xlabel('Volume ($\AA^3$/f.u.)')
plt.ylabel('Total energy (eV/f.u.)')
The solutions are at V = [31.67490656 34.60893508]
Anatase pressure: 4.524949423600637 GPa
Rutile pressure: 4.524949423742939 GPa
Text(0, 0.5, 'Total energy (eV/f.u.)')
At a pressure of 4.5 GPa, we expect that anatase will start converting into rutile. Along this common tangent, a mixture of the two phases will be more stable than either pure phase.
################ add literature discussion
there is some controversy about the most stable phase. add discussion here.
################ END