Adsorption energy calculation using DFT

6.2. Adsorption energy calculation using DFT#

6.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}\)

6.2.2. Calculating Adsorption Energy in QE#

Consider O₂ 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 O₂ 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 O₂ 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()

O₂ 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()

O₂ 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))

Adsorption energy calculation using DFT

Contents

6.2. Adsorption energy calculation using DFT#

6.2.1. Adsorption Energy#

6.2.2. Calculating Adsorption Energy in QE#