Bulk reaction energies#
Alloy formation energies#
In this section we will consider how to calculate the formation energy of an fcc Cu-Pd alloy and how to use that information to discuss relative stabilities. The kinds of questions we can easily answer are:
Is the formation of an alloy at a particular composition and structure energetically favorable?
Given two alloy structures at the same composition, which one is more stable?
Given a set of alloy structures at different compositions, which ones are stable with respect to phase separation?
Each of these questions is answered by calculating the formation energy of the alloy from the parent metals. Thus, we will need the total energies of fcc Cu and fcc Pd. To get started. We get those first. Rather than compute a full equation of state for these, we will rely on the built in unit cell optimization algorithm in VASP (ISIF=3).
Basic alloy formation energy#
# get bulk Cu and Pd energies. <<pure-metal-components>>
from vasp import Vasp
from ase import Atom, Atoms
atoms = Atoms([Atom('Cu', [0.000, 0.000, 0.000])],
cell= [[ 1.818, 0.000, 1.818],
[ 1.818, 1.818, 0.000],
[ 0.000, 1.818, 1.818]])
cuc = Vasp(label='bulk/alloy/cu',
xc='PBE',
encut=350,
kpts=[13, 13, 13],
nbands=9,
ibrion=2,
isif=3,
nsw=10,
atoms=atoms)
cu = cuc.potential_energy
atoms = Atoms([Atom('Pd', [0.000, 0.000, 0.000])],
cell=[[ 1.978, 0.000, 1.978],
[ 1.978, 1.978, 0.000],
[0.000, 1.978, 1.978]])
pd = Vasp(label='bulk/alloy/pd',
xc='PBE',
encut=350,
kpts=[13, 13, 13],
nbands=9,
ibrion=2,
isif=3,
nsw=10,
atoms=atoms).potential_energy
print(f'Cu energy = {cu} eV')
print(f'Pd energy = {pd} eV')
Note that the Pd energy is more negative than the Cu energy. This does not mean anything significant. We cannot say Pd is more stable than Cu; it is not like Cu could transmutate into Pd!
Next, we will consider a question like which of two structures with composition of CuPd is more stable. These coordinates for these structures came from research of the author. The approach is pretty general, you must identify the coordinates and unit cell of the candidate structure, and then run a calculation to find the optimized geometry and unit cell. This may take some work, as previously described in the multistep process for optimizing a bulk system. Here the geometry is pretty close to optimized, so we can use the VASP optimization routines. We consider two structures with composition CuPd.
from vasp import Vasp
from ase import Atom, Atoms
atoms = Atoms([Atom('Cu', [0.000, 0.000, 0.000]),
Atom('Pd', [-1.652, 0.000, 2.039])],
cell= [[0.000, -2.039, 2.039],
[0.000, 2.039, 2.039],
[-3.303, 0.000, 0.000]])
calc = Vasp(label='bulk/alloy/cupd-1',
xc='PBE',
encut=350,
kpts=[12, 12, 8],
nbands=17,
ibrion=2,
isif=3,
nsw=10,
atoms=atoms)
cupd1 = atoms.get_potential_energy()
atoms = Atoms([Atom('Cu', [-0.049, 0.049, 0.049]),
Atom('Cu', [-11.170, 11.170, 11.170]),
Atom('Pd', [-7.415, 7.415, 7.415]),
Atom('Pd', [-3.804 , 3.804, 3.804])],
cell=[[-5.629, 3.701, 5.629 ],
[-3.701, 5.629, 5.629 ],
[-5.629, 5.629, 3.701 ]])
calc = Vasp(label='bulk/alloy/cupd-2',
xc='PBE',
encut=350,
kpts=[8, 8, 8],
nbands=34,
ibrion=2,
isif=3,
nsw=10,
atoms=atoms)
cupd2 = atoms.get_potential_energy()
print(f'cupd-1 = {cupd1} eV')
print(f'cupd-2 = {cupd2} eV')
cupd-1 = -9.17593835 eV
cupd-2 = -18.07779325 eV
Looking at these energies, you could be tempted to say cupd-2 is more stable than cupd-1 because its energy is much lower. This is wrong, however, because cupd-2 has twice as many atoms as cupd-1. We should compare the normalized total energies, that is the energy normalized per CuPd formula unit, or as an alternative the number of atoms in the unit cell. It does not matter which, as long as we normalize consistently. It is conventional in alloy calculation to normalize by the number of atoms in the unit cell.
from vasp import Vasp
calc = Vasp(label='bulk/alloy/cupd-1')
atoms = calc.load_atoms()
e1 = atoms.get_potential_energy()/len(atoms)
calc = Vasp(label='bulk/alloy/cupd-2')
atoms = calc.load_atoms()
e2 = atoms.get_potential_energy()/len(atoms)
print(f'cupd-1: {e1} eV/atom')
print(f'cupd-2: {e2} eV/atom')
After normalizing by number of atoms, we can see that cupd-1 is a more stable structure. However, we are looking at total energies, and we might ask: is cupd-1 more stable than an unreacted mixture of the parent compounds, fcc Cu and Pd? In other words, is the following reaction exothermic:
Cu + Pd \(\rightarrow\) CuPd for the two configurations we examined? Below, we show some pretty general code that computes these formation energies, and normalizes them by the number of atoms in the unit cell.
from vasp import Vasp
# bulk energy 1
calc = Vasp(label='bulk/alloy/cu')
atoms = calc.load_atoms()
cu = atoms.get_potential_energy()/len(atoms)
# bulk energy 2
calc = Vasp(label='bulk/alloy/pd')
atoms = calc.load_atoms()
pd = atoms.get_potential_energy()/len(atoms)
calc = Vasp(label='bulk/alloy/cupd-1')
atoms = calc.load_atoms()
e1 = atoms.get_potential_energy()
# subtract bulk energies off of each atom in cell
for atom in atoms:
if atom.symbol == 'Cu':
e1 -= cu
else:
e1 -= pd
e1 /= len(atoms) # normalize by number of atoms in cell
calc = Vasp(label='bulk/alloy/cupd-2')
atoms = calc.load_atoms()
e2 = atoms.get_potential_energy()
for atom in atoms:
if atom.symbol == 'Cu':
e2 -= cu
else:
e2 -= pd
e2 /= len(atoms)
print(f'Delta Hf cupd-1 = {e1:1.2f} eV/atom')
print(f'Delta Hf cupd-2 = {e2:1.2f} eV/atom')
The answer is yes. Both structures are energetically more favorable than an equal composition mixture of the parent metals. The heat of formation for both structures is exothermic, but the cupd-1 structure is more stable than the cupd-2 structure. This is shown conceptually in Figure ref:fig:alloy1.
We will now examine another structure at another composition and its stability.
from vasp import Vasp
from vasp.exceptions import VaspNotConverged
from ase import Atom, Atoms
# parent metals
atoms = Vasp(label='bulk/alloy/cu').load_atoms()
cu = atoms.get_potential_energy() / len(atoms)
atoms = Vasp(label='bulk/alloy/pd').load_atoms()
pd = atoms.get_potential_energy() / len(atoms)
atoms = Atoms([Atom('Cu', [-3.672, 3.672, 3.672]),
Atom('Cu', [0.000, 0.000, 0.000]),
Atom('Cu', [-10.821, 10.821, 10.821]),
Atom('Pd', [-7.246, 7.246, 7.246])],
cell=[[-5.464, 3.565, 5.464],
[-3.565, 5.464, 5.464],
[-5.464, 5.464, 3.565]])
calc = Vasp(label='bulk/alloy/cu3pd-1',
xc='PBE',
encut=350,
kpts=[8, 8, 8],
nbands=34,
ibrion=2,
isif=3,
nsw=10,
atoms=atoms)
try:
e3 = atoms.get_potential_energy()
for atom in atoms:
if atom.symbol == 'Cu':
e3 -= cu
else:
e3 -= pd
e3 /= len(atoms)
print(f'Delta Hf cu3pd-1 = {e3:1.2f} eV/atom')
except VaspNotConverged as e:
print(f'Calculation did not converge: {e}')
print('This structure may need different VASP parameters to converge.')
Delta Hf cu3pd-1 = -0.02 eV/atom
The formation energy is slightly exothermic, which means the structure is more stable than a mixture of the parent metals. However, let us consider whether the structure is stable with respect to phase separation into pure Cu and the cupd-1 structure. We define the following quantities:
\(H_{f,Cu}\) = 0.0 eV/atom, \(x_0\) = 0, \(H_{f,cupd-1}\) = -0.12 eV/atom, \(x_3\) = 0.5.
The composition weighted average at \(x_{Pd}=0.25\) is:
\(H_f = H_{f,Cu} + \frac{x0-x}{x0-x3}(H_{f,cupd-1} - H_{f,Cu})\)
x0 = 0.0; x3 = 0.5; x = 0.25;
Hf1 = 0.0; Hf3 = -0.12;
print(f'Composition weighted average = {Hf1 + (x0 - x) / (x0 - x3) * (Hf3 - Hf1)} eV')
Composition weighted average = -0.06 eV
We find the weighted composition formation energy of pure Cu and cupd-1 is more favorable than the formation energy of cu3pd-1. Therefore, we could expect that structure to phase separate into a mixture of pure Cu and cupd-1. Schematically what we are seeing is shown in Figure \ref:fig:alloy-phase-separation.
Finally, let us consider one more structure with the Cu3Pd stoichiometry.
from vasp import Vasp
from ase import Atom, Atoms
# parent metals
cu = Vasp(label='bulk/alloy/cu')
cu_e = cu.potential_energy / len(cu.load_atoms())
pd = Vasp(label='bulk/alloy/pd')
pd_e = pd.potential_energy / len(pd.load_atoms())
atoms = Atoms([Atom('Cu', [-1.867, 1.867, 0.000]),
Atom('Cu', [0.000, 0.000, 0.000]),
Atom('Cu', [0.000, 1.867, 1.867]),
Atom('Pd', [-1.867, 0.000, 1.86])],
cell=[[-3.735, 0.000, 0.000],
[0.000, 0.000, 3.735],
[0.000, 3.735, 0.000]])
calc = Vasp(label='bulk/alloy/cu3pd-2',
xc='PBE',
encut=350,
kpts=[8, 8, 8],
nbands=34,
ibrion=2,
isif=3,
nsw=10,
atoms=atoms)
e4 = atoms.get_potential_energy()
Vasp.wait(abort=True)
for atom in atoms:
if atom.symbol == 'Cu':
e4 -= cu_e
else:
e4 -= pd_e
e4 /= len(atoms)
print(f'Delta Hf cu3pd-2 = {e4:1.2f} eV/atom')
Delta Hf cu3pd-2 = -0.10 eV/atom
This looks promising: the formation energy is much more favorable than cu3pd-1, and it is below the composition weighted formation energy of -0.06 eV/atom. Consequently, we conclude that this structure will not phase separate into a mixture of Cu and CuPd. We cannot say, however, if there is a more stable phase not yet considered, or if it might phase separate into two other phases. We also note here that we have ignored a few other contributions to alloy stability. We have only considered the electronic energy contributions to the formation energy. At temperatures above absolute zero there are additional contributions including configurational and vibrational entropy, which may stabilize some structures more than others. Finally, our analysis is limited to comparisons of the structures computed on the fcc lattice. In fact, it is known that the CuPd alloy forms a bcc structure. We did not calculate that structure, so we can not say if it is more or less stable than the obvious fcc structure we found.
The construction of alloy phase diagrams is difficult. You are always faced with the possibility that there is a phase that you have not calculated that is more stable than the ones you did calculate. One approach is to use a tool that automates the discovery of relevant structures such as the Alloy Theoretic Automated Toolkit (ATAT) cite:vandeWalle2002539,vandeWalle2009266 which uses a cluster expansion methodology.
Metal oxide oxidation energies#
We will consider here the reaction 2 Cu2O + O2 \(\rightleftharpoons\) 4 CuO. The reaction energy is:
\(\Delta E = 4E_{CuO} - 2E_{Cu_2O} - E_{O_2}\). We need to compute the energy of each species.
Cu2O calculation#
# run Cu2O calculation
from vasp import Vasp
from ase import Atom, Atoms
# http://phycomp.technion.ac.il/~ira/types.html#Cu2O
a = 4.27
atoms = Atoms([Atom('Cu', [0, 0, 0]),
Atom('Cu', [0.5, 0.5, 0.0]),
Atom('Cu', [0.5, 0.0, 0.5]),
Atom('Cu', [0.0, 0.5, 0.5]),
Atom('O', [0.25, 0.25, 0.25]),
Atom('O', [0.75, 0.75, 0.75])])
atoms.set_cell((a, a, a), scale_atoms=True)
calc = Vasp(label='bulk/Cu2O',
encut=400,
kpts=(8, 8, 8),
ibrion=2,
isif=3,
nsw=30,
xc='PBE',
atoms=atoms)
print(atoms.get_potential_energy())
print(atoms.get_stress())
-27.27469148
[-0.01018402 -0.01018402 -0.01018402 -0. -0. -0. ]
CuO calculation#
# run CuO calculation
from vasp import Vasp
from ase import Atom, Atoms
import numpy as np
# CuO
# http://cst-www.nrl.navy.mil/lattice/struk/b26.html
# http://www.springermaterials.com/docs/info/10681727_51.html
a = 4.6837
b = 3.4226
c = 5.1288
beta = 99.54/180*np.pi
y = 0.5819
a1 = np.array([0.5*a, -0.5*b, 0.0])
a2 = np.array([0.5*a, 0.5*b, 0.0])
a3 = np.array([c*np.cos(beta), 0.0, c*np.sin(beta)])
atoms = Atoms([Atom('Cu', 0.5*a2),
Atom('Cu', 0.5*a1 + 0.5*a3),
Atom('O', -y*a1 + y*a2 + 0.25*a3),
Atom('O', y*a1 - y*a2 - 0.25*a3)],
cell=(a1, a2, a3))
calc = Vasp(label='bulk/CuO',
encut=400,
kpts=(8, 8, 8),
ibrion=2,
isif=3,
nsw=30,
xc='PBE',
atoms=atoms)
print(atoms.get_potential_energy())
-19.61568557
-19.61568557
Reaction energy calculation#
from vasp import Vasp
# don't forget to normalize your total energy to a formula unit. Cu2O
# has 3 atoms, so the number of formula units in an atoms is
# len(atoms)/3.
calc = Vasp(label='bulk/Cu2O')
atoms1 = calc.load_atoms()
cu2o_energy = atoms1.get_potential_energy()
calc = Vasp(label='bulk/CuO')
atoms2 = calc.load_atoms()
cuo_energy = atoms2.get_potential_energy()
# make sure to use the same cutoff energy for the O2 molecule!
calc = Vasp(label='molecules/O2-sp-triplet-400')
atoms3 = calc.load_atoms()
o2_energy = atoms3.get_potential_energy()
calc.stop_if(None in [cu2o_energy, cuo_energy, o2_energy])
cu2o_energy /= (len(atoms1) // 3) # use integer division
cuo_energy /= (len(atoms2) // 2)
rxn_energy = 4.0 * cuo_energy - o2_energy - 2.0 * cu2o_energy
print(f'Reaction energy = {rxn_energy} eV')
Reaction energy = -2.11600154 eV
This is the reaction energy for 2 Cu2O \(\rightarrow\) 4 CuO. In cite:PhysRevB.73.195107, the experimental reaction is estimated to be about -3.14 eV.
There are a few reasons why our number does not agree with the experimental reaction energy. One reason is related to errors in the O2 dissociation energy, and another reason is related to localization of electrons in the Cu 3\(d\) orbitals cite:PhysRevB.73.195107. The first error of incorrect O2 dissociation error is a systematic error that can be corrected empirically cite:PhysRevB.73.195107. Fixing the second error requires the application of DFT+U (see [BROKEN LINK: *DFT%2BU]).
The heat of reaction is reported to be 1000 J/g product at http://onlinelibrary.wiley.com/doi/10.1002/er.4440130107/pdffor the reaction 2CuO \(\rightleftharpoons\) Cu2O + 1/2 O2.
from ase import Atoms
atoms = Atoms('Cu2O')
MW = atoms.get_masses().sum()
H = 1. # kJ/g
print(f'rxn energy = {-2 * H * MW / 96.4:1.1f} eV') # convert to eV
rxn energy = -3.0 eV
This is pretty close to the value in cite:PhysRevB.73.195107 and might need a temperature correction to get agreement at 298K.