5.2. Adsorption energy calculation using DFT#
5.2.1. Adsorption Energy#
Adsorption energy refers to the change of energy when an adsorbate is attached to the surface of an adsorbent.
Adsorption energy, \({\Delta}E_{a}\), can be obtained by calculating the difference between the energy of the adsorbed surface and the sum of the energy for each structure composing the adsorbed surface:
\({\Delta}E_{a} = E_{ads,slab} - E_{slab} - E_{ads}\)
5.2.2. Calculating Adsorption Energy in QE#
Consider O2 adsorption on Platinum (Pt) surface with constrained atom positions
For this situation, we will first build the Atoms object for Platinum slab and adsorb an oxygen molecule to the slab.
Building atoms
In the following script, a Pt slab is created. The bottom two layers of Pt (in z-direction) are fixed. The exact position to divide the atoms is determined by visualizing the ‘input.traj’ file with ase gui command. By clicking on the atoms in the pop-out window, the exact position within a cell will be displayed. After applying the constraint to the system, visualizing ‘input.traj’ file with ase gui command, you should be able to see the bottom two layers of Pt atoms have a cross on them, indicating that they are fixed. Another way to check or to manipulate the constraint is by ‘POSCAR’, where it lists the positions and constraints of all atoms in the system.
from ase.build import bulk , molecule , surface , add_adsorbate
from ase.constraints import FixAtoms
from ase.io import read
#building Platinum slab
bulk = bulk ('Pt')
miller_indices = (1, 1, 1)
layers = 4
slab = surface (bulk , indices = miller_indices , layers = layers , vacuum =10)
slab *= (2 ,2 ,1)
slab.write('slab.traj')
Build an O2 molecule with ase library.
from ase.build import molecule, bulk
from ase.data.pubchem import pubchem_atoms_search
import numpy as np
from shutil import copy
import os
#building oxygen molecule
atoms = molecule('O2')
max_pos = np.amax(atoms.get_positions(), axis=0)
atoms.set_cell(np.array([10, 10, 10]) + max_pos)
atoms.center()
atoms.write('O2.traj')
Attach the O2 to the Pt slab using ase.build.add_adsorbate().
from ase.build import bulk , molecule , surface , add_adsorbate
from ase.constraints import FixAtoms
from ase.io import read
slab = read('slab.traj')
adsorbate = molecule('O2')
add_adsorbate(slab, adsorbate, height = 2.5, position = (1.3 ,1) )
c = FixAtoms ( indices = [ atom . index for atom in slab if atom . position [2] <14 ]) #constrained atoms (fixing the bottom 2 layers)
slab.set_constraint(c)
slab.write('slab_ads.traj')
DFT calculation
In the following script, we define a SPARC calculator object and attach it to an ase.Atoms object so that ase knows how to get the data. After defining the calculator, a DFT calculation can be run. The potential energy can be obtained by the method: get_potential_energy(). The potential energy result can be stored in a text file e.g. ‘energy.txt’.
Details on how to submit these calculations can be seen in Lecture ‘Running_SPARC_on_QE’.
Clean slab calculation
from ase.io import read
from ase.optimize import BFGSLineSearch
from espresso import iEspresso as Espresso
import numpy as np
image = read('slab.traj')
# setup calculator
parameters = dict(xc='BEEF-vdW',
kpts=(4, 4, 1), #only need 1 kpt in z-direction
pw=400.,
dw=4000.,
spinpol=False,
convergence={'energy':1e-6,
'mixing':0.05,
'maxsteps':1000,
'diag':'david'},
startingwfc='atomic',
smearing='fd', #fermi-dirac electron smearing
sigma=0.1, #smearing width
dipole={'status':True}, #dipole corrections True turns them on
outdir ='esp.log')
calc = Espresso(atoms = image, **parameters)
image.set_calculator(calc)
relax = BFGSLineSearch(image,
trajectory='opt.traj',
logfile='opt.log',
restart='opt.pckl')
relax.run(fmax=0.05)
eng = image.get_potential_energy()
image.write('converged.traj')
with open('slab_energy.txt', 'w') as f:
f.write(str(eng))
image.calc.close()
O2 on the Pt slab
from ase.io import read
from ase.optimize import BFGSLineSearch
from espresso import iEspresso as Espresso
import numpy as np
image = read('slab_ads.traj')
# setup calculator
parameters = dict(xc='BEEF-vdW',
kpts=(4, 4, 1), #only need 1 kpt in z-direction
pw=400.,
dw=4000.,
spinpol=False,
convergence={'energy':1e-6,
'mixing':0.05,
'maxsteps':1000,
'diag':'david'},
startingwfc='atomic',
smearing='fd', #fermi-dirac electron smearing
sigma=0.1, #smearing width
dipole={'status':True}, #dipole corrections True turns them on
outdir ='esp.log')
calc = Espresso(atoms = image, **parameters)
image.set_calculator(calc)
relax = BFGSLineSearch(image,
trajectory='opt.traj',
logfile='opt.log',
restart='opt.pckl')
relax.run(fmax=0.05)
eng = image.get_potential_energy()
image.write('converged.traj')
with open('slab_ads_energy.txt', 'w') as f:
f.write(str(eng))
image.calc.close()
O2 gas reference energy calculation
from ase.io import read
from ase.optimize import BFGSLineSearch
from espresso import iEspresso as Espresso
import numpy as np
image = read('O2.traj')
# setup calculator
parameters = dict(xc='BEEF-vdW',
kpts=(4, 4, 1), #only need 1 kpt in z-direction
pw=400.,
dw=4000.,
spinpol=False,
convergence={'energy':1e-6,
'mixing':0.05,
'maxsteps':1000,
'diag':'david'},
startingwfc='atomic',
smearing='fd', #fermi-dirac electron smearing
sigma=0.1, #smearing width
dipole={'status':True}, #dipole corrections True turns them on
outdir ='esp.log')
calc = Espresso(atoms = image, **parameters)
image.set_calculator(calc)
relax = BFGSLineSearch(image,
trajectory='opt.traj',
logfile='opt.log',
restart='opt.pckl')
relax.run(fmax=0.05)
eng = image.get_potential_energy()
image.write('converged.traj')
with open('O2_energy.txt', 'w') as f:
f.write(str(eng))
image.calc.close()
Analysis of adsorption energies
import os
from ase.io import read
def readFile(path):
f = open(path, 'r')
content = f.readlines()
return content
def read_energy(path):
energy = float(readFile(path)[0])
return energy
E_slab = read_energy('slab_energy.txt')
E_slab_ads = read_energy('slab_ads_energy.txt')
E_O2 = read_energy('O2_energy.txt')
# E_adsorption is the adsorption energy
E_adsorption = E_slab_ads - E_slab - E_O2
print('The adsorption energy of O2 on Pt slab is {0:1.3f} eV'.format(E_adsorption))