9. Adsorption Energy
The adsorption energy denotes the energy variation associated with the attachment of a molecule to a solid surface. This parameter is fundamental for the quantitative analysis of interfacial molecular interactions and plays a pivotal role in governing key physicochemical properties, including particle adhesion, wettability, and dispersibility [GCK]. In this tutorial we’ll show how to use the RESCU code to calculate the adsorption energy of an water molecule on the anatase \(TiO_{2} (001)\) surface, as shown in Fig. 9.1.
Fig. 9.1 The \(TiO_{2}\) anatase along (001) slab with attached \(H_{2}O\) molecule.
The calculation steps are given as follows:
Build the structures:
Generate the \(TiO_{2}\) conventional cell;
Create the surface slab of \(TiO_{2}\) along (001);
Place the \(H_{2}O\) molecule on the \(TiO_{2}\) (001) surface.
Structural Relaxation:
Relax the surface layer together with \(H_{2}O\) molecule constrained the bottom layers.
Self-consistent (SCF) calculations.
Simulation results:
Bond length and Adsorption energy analysis.
This simple model of \(TiO_{2}/H_{2}O\) interaction will be useful to explore all methodology steps.
9.1. Build the structures
Generating the conventional cell
If the user do not have the conventional cell of \(TiO_{2}\), or no database with the avaible structure, it conventional cell could be automatically generated from scratch using the provided Python script gen_TiO2.py.
To execute this script, the user must install the required Python packages, which include the following modules:
pymatgen
ase (Atomic Simulation Environment)
numpy
These packages can be installed using the command below:
pip install ase numpy pymatgen
Note
It is highly recommended to use an isolated environment, such as a Python virtual environment or a Conda environment, when installing additional Python modules on your system, in order to prevent potential conflicts with dependencies from other projects.
Generate the \(TiO_{2}\) conventional cell in CIF format (TiO2_anatase_I41amd.cif) by running:
python gen_TiO2.py
The structure is constructed using lattice parameters of a = b = 3.784 Å; c = 9.514 Å and space group \(I4_{1}/amd\) (Fig. 9.1.2). Nonetheless, this structure may serve as an initial configuration for a full relaxation (see the Equation of State tutorial).
Fig. 9.1.2 The \(TiO_{2}\) conventional cell.
The surface TiO2 (001)
Once the \(TiO_{2}\) crystal (bulk structure) has been set up, the user can generate the slab surface along the [0 0 1] direction using the provided Python script gen_TiO2_slab.py.
In this script the parameters were seted as follows:
The Miller indices were set to h = 0, k = 0, l = 1, in order to cleave the \(TiO_{2}\) crystal along the [0 0 1] direction.
Slab thickness was set to 12 Å, corresponding to 5 atomic layers;
A vacuum region of 15 Å was added along the z-axis to prevent interactions between periodic images;
To create a \(2\times 2\times 1\) supercell, the following lattice vector multiplications were applied:
u = (2a, 0b, 0c),
v = (0a, 2b, 0c),
w = (0a, 0b, 1c)
Generate the \(TiO_{2}\) (001) slab in CIF format (TiO2_anatase_001_slab_2x2x1.cif) by running:
python gen_TiO2_slab.py
The surface TiO2 (001)/H2O
Next, we’ll place an \(H_{2}O\) molecule on top of the previously generated \(TiO_{2}\) (001) slab, as illustrated in Fig. 9.1.3.
Launch the ASE GUI with the following command:
ase gui TiO2_anatase_001_slab_2x2x1.cif
This will open the ASE GUI and load the file TiO2_anatase_001_slab_2x2x1.cif.
In the top menu, go to Edit ➟ add atoms:
A new window will pop up. In the Add field, type H2O, and set the coordinates to 2, 2, 20.6;
Click in the Add buttom.
Adjust the view to see the molecule better:
Go to the menu bar and click View ➟ Change View ➟ zx-plane to change the visualization of the structure.
Use the left mouse button to click and drag, and select only the \(H_{2}O\) molecule.
Rotate the H₂O molecule:
Back to the menur bar and click Tools ➟ Rotate selected atoms;
Use the up or down arrow keys from keyboard to rotate the selected molecule by \(180^{o}\).
Go to File ➟ Save to export the structure in the CIF format and name the file as
tio2_water.rx.cif.
Fig. 9.1.3 Graphical interface to visualize and manipulate molecular atoms objects (ASE GUI).
9.2. Structural Relaxation
To begin we’ll convert the tio2_water.rx.cif file into a RESCU input file using the provided Python script cif2RESCU.py, as shown bellow:
python cif2RESCU.py -i tio2_water.rx.cif
After running this command, two files should be generated:
atoms.xyz(sorted by z-axis)
These files will be used for the structural relaxation step.
Warning
To ensure a clean and organized workflow, we strongly recommend creating a separate folder for each step of the calculation process.
Go to the folder where all the input files are located. Before running the calculation, we’ll review the contents of the input files.
Atoms.xyz file
63
O 0.000000000000000 1.891999999999999 8.142191009999998
O 0.000000000000000 5.676000000000000 8.142191009999998
O 3.783999999999998 1.891999999999999 8.142191009999998
O 3.783999999999998 5.676000000000000 8.142191009999998
Ti 5.676000000000000 5.676000000000000 8.542730340000000
Ti 5.676000000000000 1.891999999999999 8.542730340000000
Ti 1.891999999999999 5.676000000000000 8.542730340000000
Ti 1.891999999999999 1.891999999999999 8.542730340000000
O 1.891999999999999 0.000000000000000 8.943269670000001
O 1.891999999999999 3.783999999999998 8.943269670000001
O 5.676000000000000 0.000000000000000 8.943269670000001
O 5.676000000000000 3.783999999999998 8.943269670000001
O 5.676000000000000 5.676000000000000 10.520690940000001
O 5.676000000000000 1.891999999999999 10.520690940000001
O 1.891999999999999 5.676000000000000 10.520690940000001
O 1.891999999999999 1.891999999999999 10.520690940000001
Ti 1.891999999999999 0.000000000000000 10.921230270000001
Ti 1.891999999999999 3.783999999999998 10.921230270000001
Ti 5.676000000000000 0.000000000000000 10.921230270000001
Ti 5.676000000000000 3.783999999999998 10.921230270000001
O 0.000000000000000 0.000000000000000 11.321769599999996
O 0.000000000000000 3.783999999999998 11.321769599999996
O 3.783999999999998 0.000000000000000 11.321769599999996
O 3.783999999999998 3.783999999999998 11.321769599999996
O 5.676000000000000 3.783999999999998 12.899190870000000
O 5.676000000000000 0.000000000000000 12.899190870000000
O 1.891999999999999 3.783999999999998 12.899190870000000
O 1.891999999999999 0.000000000000000 12.899190870000000
Ti 3.783999999999998 3.783999999999998 13.299730199999997
Ti 3.783999999999998 0.000000000000000 13.299730199999997
Ti 0.000000000000000 3.783999999999998 13.299730199999997
Ti 0.000000000000000 0.000000000000000 13.299730199999997
O 3.783999999999998 5.676000000000000 13.700269800000001
O 3.783999999999998 1.891999999999999 13.700269800000001
O 0.000000000000000 5.676000000000000 13.700269800000001
O 0.000000000000000 1.891999999999999 13.700269800000001
O 3.783999999999998 3.783999999999998 15.277690800000000
O 3.783999999999998 0.000000000000000 15.277690800000000
O 0.000000000000000 3.783999999999998 15.277690800000000
O 0.000000000000000 0.000000000000000 15.277690800000000
Ti 0.000000000000000 1.891999999999999 15.678230399999999
Ti 0.000000000000000 5.676000000000000 15.678230399999999
Ti 3.783999999999998 1.891999999999999 15.678230399999999
Ti 3.783999999999998 5.676000000000000 15.678230399999999
O 5.676000000000000 5.676000000000000 16.078769730000001
O 5.676000000000000 1.891999999999999 16.078769730000001
O 1.891999999999999 5.676000000000000 16.078769730000001
O 1.891999999999999 1.891999999999999 16.078769730000001
O 3.783999999999998 5.676000000000000 17.656191000000000
O 3.783999999999998 1.891999999999999 17.656191000000000
O 0.000000000000000 5.676000000000000 17.656191000000000
O 0.000000000000000 1.891999999999999 17.656191000000000
Ti 1.891999999999999 1.891999999999999 18.056730330000001
Ti 1.891999999999999 5.676000000000000 18.056730330000001
Ti 5.676000000000000 1.891999999999999 18.056730330000001
Ti 5.676000000000000 5.676000000000000 18.056730330000001
O 1.891999999999999 0.000000000000000 18.457269659999994
O 1.891999999999999 3.783999999999998 18.457269659999994
O 5.676000000000000 0.000000000000000 18.457269659999994
O 5.676000000000000 3.783999999999998 18.457269659999994
O 1.999999999999999 1.999999999999999 20.202460666666660
H 1.999999999999999 1.236761000000001 20.798769666666658
H 1.999999999999999 2.763238999999997 20.798769666666658
Important
In this example we’ll relax only two layers of \(TiO_{2}\) (001) slab along with the \(H_{2}O\) molecule, while keeping the innermost layers (those buried beneath the surface) fixed. This is a common practice to ensure that the structural characteristics of the bulk are preserved during surface relaxation.
For the purpose of this tutorial, we are using a relatively thin slab (5 layers). However, in a realistic simulation that aims to improve the behavior of a true surface-bulk system, a greater number of layers could be included.
The atoms.xyz file contains a total of 63 atoms, divided as follows:
From 1 to 36: Represent the innermost layers of the slab, which will remain fixed during relaxation.
From 37 to 63: Represent the surface region plus molecule, i.e., 24 atoms from \(TiO_{2}\) slab and 3 atoms from the \(H_{2}O\) molecule. These atoms will be relaxed during the simulation.
Now, we are ready to start the calculation. However, before doing that we’ll briefly present the simulation parameters specified in the generated file.
The relax input file
info.calculationType = 'relaxation'
info.savepath = './results/tio2_water.rx.mat'
info.outfile = 'tio2_water.rx.out'
units.length = 'Angstrom'
atom.xyz = 'atoms.xyz'
sr.tol = 1e-4
sr.moveableAtomList = [37:63]
sr.relaxedstructure = 'relaxedstructure.xyz'
domain.latvec = [
[7.568000000000 0.000000000000 0.000000000000];
[0.000000000000 7.568000000000 0.000000000000];
[0.000000000000 0.000000000000 27.000000000000]
]
domain.boundary = [1 1 2]
LCAO.status = 1
kpoint.sigma = 0.05
eigensolver.emptyBand = 32
units.energy = 'eV'
kpoint.gridn = [2,2,1]
domain.lowres = 0.35
functional.libxc = true % required for vdw
functional.list = {'XC_GGA_X_OPTPBE_VDW','XC_GGA_C_OP_PBE'}
option.maxSCFiteration = 200
mixing.type = 'density'
mixing.method = 'pulay'
mixing.tol = [1e-05,1e-05]
mixing.maxhistory = 20
spin.type = 'degenerate'
element(1).species = 'H'
element(1).path = './pseudos/H_AtomicData.mat'
element(2).species = 'O'
element(2).path = './pseudos/O_AtomicData.mat'
element(3).species = 'Ti'
element(3).path = './pseudos/Ti_AtomicData.mat'
The system information and simulation parameters are defined by keywords. These keywords generally have self-explanatory names, but you may refer to the input reference for more information.
The keywords specify the following:
info.calculationTypeThis input parameter specifies the type of the calculation that RESCU should perform.info.savepathThe path and name prefix of the file where the results will be writer. RESCU appends the rank number to the file name prefix, which for serial calculations is always 0.info.outfileFile in which writes the simulation information.units.lengthUnit of length used for all structural inputs (atomic positions, lattice vectors, cell size). Accepted values are bohr (default) and angstrom.atom.xyzIs an array of row-vectors containing the position of each atom.sr.tolForce tolerance. The relaxation process stops once the residual forces on the movable atoms are less than the specified tolerance. The default unit is Hartree/Bohr.sr.moveableAtomListArray containing the indices of the relaxed atoms. The default is to relax all atoms except the first atom.sr.relaxedstructureFile in which the final relaxed atomic coordinates can be found.domain.latvecLattice vectors. The matrix must contain three row-vectors.domain.boundaryDetermine the boundary conditions along each dimension.LCAO.statusIf LCAO.status = true, RESCU performs a calculation using the atomic orbital basis. This is generally faster than performing a purely real space calculation, but the atomic orbital basis must be validated for correctness.kpoint.sigmaUse Fermi-Dirac smearing if option.sigma \(>\) 0; sigma in Hartree by default.eigensolver.emptyBandNumber of bands included in the calculation in addition to the number of valence bandsunits.energyDefines the default units of length in the input file.kpoint.gridnsets the number of k-points used to sample the first Brillouin zone.domain.lowrescontrols the real space resolution.functional.libxcUse the LibXC wrapper whenfunctional.libxc = trueto include functional from the LibXC.functional.listList of exchange-correlation functional used in the calculation.option.maxSCFiterationMaximal number of self-consistent steps; the Kohn-Sham solver will return an answer even if it has not converged.mixing.typeDetermine whether to mix the density, density-matrix or the potential.mixing.methodIf equal to linear and pulay the mixing method of the same name will be used. If equal to broyden method of Srivastava will be condidered.mixing.tolFor self-consistent (ground state) calculations, the self-consistent procedure stops if the charge variation per valence electron is less than mixing.tol(1) and the total energy variation per valence electron is less than mixing.tol(2). Type [1,2].mixing.maxhistoryMaximal number of trial functions retained to compute the new trial function.spin.typeDetermines how to treat spin: degenerate, collinear, non-collinear.element(1).speciesIt is a tag for the element structure. It is usually the chemical symbol of the corresponding element.element(1).pathis a string containing the path to the pseudo-potential file associated with the element.
Warning
It is highly recommended to compile the RESCU code with LibXC support enabled. The LibXC library provides access to a wide range of exchange-correlation functionals, including both commonly used and advanced functionals that improve the default RESCU implementation.
By enabling LibXC during compilation, users can:
Take advantage of a broader set of XC functionals, such as non-local van der Waals (vdW-DF), hybrid and meta-GGA density functionals.
Improve the accuracy and flexibility of simulations, especially for surface calculations, molecular adsorption, and complex materials.
To compile RESCU with LibXC please look in our documentation Compiling the LibXC.
Now, we are ready to start the calculations. Then copy the H, O, and Ti Pseudopotential file from the database to ./pseudos directory. Execute RESCU and use the -i option to point to the input file.
rescu -i tio2_water.rx.input
After the calculation is finished, the following output files are generated:
relaxationprocess.txt;
tio2_water.rx.out;
tio2_water.rx.h5;
relaxedstructure.xyz;
tio2_water.rx.mat.
The file relaxationprocess.txt contains the geometry optimization steps and can be useful for visualizing how the \(H_{2}O\) molecule relaxes on the \(TiO_{2}\) (001) anatase surface, as illustrated in Fig. 9.2.2.
Fig. 9.2.2 \(H_{2}O\) molecule relaxation steps on the \(TiO_{2}\) (001) anatase surface.
The total energy of the optimized system can be loaded from the tio2_water.rx.mat file. In the MATLAB environment, run the following command:
load('tio2_water.rx.mat')
energy.Etot
This will return:
ans =
-883.8213
As shown above, the total energy of the system is -883.8213 Hartree. The next step is to perform two additional single-point calculations: one for the isolated \(H_{2}O\) molecule and another for the isolated \(TiO_{2}\) (001) surface.
9.3. Single-point calculations
To compute the interaction energy, you need the total energy of each isolated component (\(TiO_{2}\) surface and \(H_{2}O\) molecule). Follow the steps below:
Start with the relaxedstructure.xyz file obtained from the previous calculation.
Create two new structure files by copying and editing:
tio2.xyz: keep only the first 60 atoms (TiO₂ surface)
water.xyz: keep only the last 3 atoms (H₂O molecule)
Update the atom count at the top of each file accordingly.
Prepare the input files for each component by duplicating tio2_water.rx.input:
Rename the copies as:
tio2.scf.input for the \(TiO_{2}\) surface
water.scf.input for the \(H_{2}O\) molecule
Edit both input files:
Set the calculation type:
info.calculationType = 'self-consistent'
Remove lines related to structural relaxation:
sr.tol = 1e-4 sr.moveableAtomList = [37:63] sr.relaxedstructure = 'relaxedstructure.xyz'
Update input file settings:
info.outfile: set to a descriptive name (e.g., tio2.out, water.out)
info.savepath: set to the desired result file (e.g., tio2.mat, water.mat)
atoms.xyz: change to the correct structure file (tio2.xyz or water.xyz)
At the end of the file, ensure that only the correct atomic elements are listed and the structure file matches the component being calculated
Alternatively, you can download all the prepared files from the following link:
Go to the directory \(TiO_{2}\) where the input files are located, and execute RESCU as follows:
rescu -i tio2.scf.input
Repeat the same procedure for the water calculation. As with the previous energy check, the total energy for each isolated system can be retrieved from the .mat output.
9.4. Simulation results
Bond lengthds analysis
We’ll begin by analyzing the optimized bond lengths of the composed system, as shown in the Optimized bond lengths table below:
Bond (Å) |
RESCU / GGA (PBE-vdW, NC) |
Literature / GGA (Ultrasoft) [VSRK] |
\(\Delta\) (%) |
|---|---|---|---|
Ti–O |
2.247 |
2.33 |
3.56 |
O–H |
1.082 |
1.05 |
3.05 |
H···O (H-bond) |
1.598 |
1.55 |
3.10 |
As one can see, the bond lengths obtained using RESCU with the GGA OptPBE-vdW functional and norm-conserving pseudopotentials show good agreement with literature values calculated using GGA and ultrasoft pseudopotentials, with relative deviations below 4%. The largest deviation was found for the Ti-O bond (3.56%), followed by the O-H (3.05%) and H···O hydrogen bond (3.10%). These small differences can be attributed to the distinct treatment of electron-core interactions, especially the use of norm-conserving pseudopotentials in RESCU, which are typically more rigid, and the inclusion of van der Waals corrections in the OptPBE-vdW functional, which were not present in the reference calculation. Overall, the results indicate that the RESCU approach is reliable for describing the structural features of the TiO₂/H₂O system.
Adsorption energy analysis
When modeling surfaces with adsorbed molecules (e.g., \(H_{2}O\) on \(TiO_{2}\)), van der Waals forces non-covalent dispersion interactions play a crucial role in determining both structural and energetic properties. The adsorption energy is calculated as:
Where:
\(E_{\text{adsorbate}}\) is the energy of the isolated \(H_{2}O\) molecule in vacuum.
\(E_{\text{surface}}\) is the energy of the clean anatase (001):math:TiO_{2} surface.
\(E_{\text{adsorbate/surface}}\) is the total energy of the adsorbed system.
Reference |
Method |
\(E_{ads}\) (kcal/mol) |
\(\Delta\) (%) |
|---|---|---|---|
Vittadini et al., 1998 [VSRK] |
PW: DFT-GGA (Car–Parrinello) |
-19 |
72.73 |
Erdogan et al., 2010 [EO] |
VASP / GGA (PW91) |
-15 |
36.36 |
Srnak et al., 1992 [SDC] |
Experimental |
-11 |
– |
This tutorial |
RESCU / GGA (OptPBE-vdW, NC) |
-11 |
0.00 |
As one can see in Table 2: Adsorption energy with RESCU using the OptPBE-vdW functional and norm-conserving pseudopotentials is -11 kcal/mol, in perfect agreement with the experimental value. In contrast, other DFT-based studies reported significantly stronger adsorption, with deviations of 36.36% [EO] and 72.73% [EO]. The goodd agreement obtained here highlights the importance of including van der Waals interactions for accurate modeling of physisorption on oxide surfaces.
S. Gim, K. J. Cho, H.-K. Lim, and H. Kim. Structure, Dynamics, and Wettability of Water at Metal Interfaces. Sci. Rep. 9, (2019), p. 14805.
A. Vittadini, A. Selloni, F. P. Rotzinger, and M. Grätzel. Structure and Energetics of Water Adsorbed at TiO2 Anatase (101) and (001) Surfaces. Phys. Rev. Lett. 81, (1998), p. 2954.
R. Erdogan, I. Onal. An ONIOM and DFT study of water and ammonia adsorption on anatase TiO2 (001) Cluster nternational Journal of Quantum Chemistry 111, (2011), p. 2149.
T.Z. Srnak, J.A. Dumesic, B.S. Clausen, E. Tornqvist, N.Y. Topsoe. Temperature-programmed desorption/reaction and in situ spectroscopic studies of vanadia/titania for catalytic reduction of nitric oxide J. Catal. 135, (1992), p. 246.