Slab + mol energy¶
Various processes are handled during reactions on catalysts, for example
Organic molecules are adsorbed on the catalyst surface
Reaction proceeds on the surface
After the reaction is completed, the product is desorbed from the catalyst surface
Since reactants and products adsorb and desorb on the catalyst, it is necessary to deal with adsorption structures (structures in which molecules exist on the surface structure).
In this section, we introduce the adsorption energy, the energy associated with the adsorption structure.
Adsorption energy
[1]:
from ase import Atoms
from ase.build import bulk
from ase.constraints import ExpCellFilter, StrainFilter
from ase.optimize import LBFGS, FIRE
import pfp_api_client
from pfp_api_client.pfp.calculators.ase_calculator import ASECalculator
from pfp_api_client.pfp.estimator import Estimator, EstimatorCalcMode
print(f"pfp_api_client: {pfp_api_client.__version__}")
# estimator = Estimator(calc_mode=EstimatorCalcMode.CRYSTAL, model_version="latest")
estimator = Estimator(calc_mode=EstimatorCalcMode.CRYSTAL_U0, model_version="v3.0.0")
calculator = ASECalculator(estimator)
pfp_api_client: 1.5.0
Adsorption energy¶
The adsorption energy \(E_{\rm{adsorption}}\) represents the energy difference before and after a molecule is adsorbed on a surface structure.
\[E_{\rm{adsorption}} = E_{\rm{slab+mol}} - (E_{\rm{slab}} + E_{\rm{mol}})\]
where \(E_{\rm{slab}}\) is the energy of the slab alone, \(E_{\rm{mol}}\) is the energy of the molecule alone, and \(E_{\rm{slab+mol}}\) is the energy of the adsorption structure with the molecule attached to the slab.
CO on Pt¶
[2]:
from ase.optimize import LBFGS
from ase.build import fcc111, molecule, add_adsorbate
from ase.constraints import ExpCellFilter, StrainFilter
def get_opt_energy(atoms, fmax=0.001, opt_mode: str = "normal"):
atoms.set_calculator(calculator)
if opt_mode == "scale":
opt1 = LBFGS(StrainFilter(atoms, mask=[1, 1, 1, 0, 0, 0]))
elif opt_mode == "all":
opt1 = LBFGS(ExpCellFilter(atoms))
else:
opt1 = LBFGS(atoms)
opt1.run(fmax=fmax)
return atoms.get_total_energy()
First, the Pt structure is prepared using the bulk, and the appropriate cell size is found by optimizing the cell size.
[3]:
bulk_atoms = bulk("Pt", cubic=True)
bulk_atoms.cell
[3]:
Cell([3.92, 3.92, 3.92])
[4]:
bulk_atoms = bulk("Pt", cubic=True)
bulk_atoms.calc = calculator
E_bulk = get_opt_energy(bulk_atoms, fmax=1e-4, opt_mode="scale")
E_bulk
Step Time Energy fmax
*Force-consistent energies used in optimization.
LBFGS: 0 01:40:53 -21.749339* 6.1356
LBFGS: 1 01:40:53 -21.296367* 15.5537
LBFGS: 2 01:40:53 -21.809602* 1.0217
LBFGS: 3 01:40:54 -21.811500* 0.1569
LBFGS: 4 01:40:54 -21.811548* 0.0021
LBFGS: 5 01:40:54 -21.811543* 0.0001
[4]:
-21.811543123706095
[5]:
bulk_atoms.cell
[5]:
Cell([3.9674290384673356, 3.967428901123176, 3.967429150600858])
[6]:
import numpy as np
a = np.mean(np.diag(bulk_atoms.cell))
a
[6]:
3.9674290300637893
You can see a slight change in value from the ASE default of 3.92A. We will use this cell size to create the adsorption structure below.
[7]:
slab = fcc111("Pt", a=a, size=(4, 4, 4), vacuum=40.0, periodic=True)
slab.calc = calculator
[8]:
E_slab = get_opt_energy(slab, fmax=1e-4, opt_mode="normal")
Step Time Energy fmax
*Force-consistent energies used in optimization.
LBFGS: 0 01:40:58 -329.391341* 0.3224
LBFGS: 1 01:40:58 -329.444661* 0.2434
LBFGS: 2 01:40:58 -329.511118* 0.0433
LBFGS: 3 01:40:58 -329.512264* 0.0223
LBFGS: 4 01:40:58 -329.512625* 0.0222
LBFGS: 5 01:40:59 -329.514918* 0.0108
LBFGS: 6 01:40:59 -329.515182* 0.0054
LBFGS: 7 01:40:59 -329.515256* 0.0007
LBFGS: 8 01:40:59 -329.515224* 0.0000
[9]:
mol = molecule("CO")
mol.calc = calculator
E_mol = get_opt_energy(mol, fmax=1e-4)
Step Time Energy fmax
*Force-consistent energies used in optimization.
LBFGS: 0 01:41:02 -11.497237* 0.8202
LBFGS: 1 01:41:02 -11.484529* 1.9590
LBFGS: 2 01:41:03 -11.500196* 0.0298
LBFGS: 3 01:41:03 -11.500201* 0.0011
LBFGS: 4 01:41:03 -11.500198* 0.0000
Creating adsorption structures¶
Adsorption structures can be created by using the add_adsorbate method to place molecules on the slab.
[10]:
%%time
slab_mol = fcc111("Pt", a=a, size=(4, 4, 4), vacuum=40.0)
slab_mol.pbc = True
mol2 = molecule("CO")
add_adsorbate(slab_mol, mol2, 3.0, "bridge")
E_slab_mol = get_opt_energy(slab_mol, fmax=1e-3)
Step Time Energy fmax
*Force-consistent energies used in optimization.
LBFGS: 0 01:41:04 -341.994179* 3.7071
LBFGS: 1 01:41:04 -342.131315* 4.4264
LBFGS: 2 01:41:04 -342.428967* 3.3710
LBFGS: 3 01:41:04 -342.724264* 0.9431
LBFGS: 4 01:41:05 -342.799992* 0.7936
LBFGS: 5 01:41:06 -342.828825* 0.4370
LBFGS: 6 01:41:07 -342.872854* 0.2373
LBFGS: 7 01:41:07 -342.876425* 0.1117
LBFGS: 8 01:41:07 -342.880139* 0.0750
LBFGS: 9 01:41:07 -342.883218* 0.0958
LBFGS: 10 01:41:07 -342.885115* 0.1009
LBFGS: 11 01:41:07 -342.886205* 0.0659
LBFGS: 12 01:41:07 -342.887150* 0.0422
LBFGS: 13 01:41:07 -342.887979* 0.0357
LBFGS: 14 01:41:08 -342.888542* 0.0436
LBFGS: 15 01:41:08 -342.888831* 0.0318
LBFGS: 16 01:41:08 -342.889109* 0.0229
LBFGS: 17 01:41:08 -342.889403* 0.0221
LBFGS: 18 01:41:08 -342.889636* 0.0267
LBFGS: 19 01:41:09 -342.889776* 0.0217
LBFGS: 20 01:41:09 -342.889900* 0.0178
LBFGS: 21 01:41:09 -342.890034* 0.0165
LBFGS: 22 01:41:10 -342.890133* 0.0168
LBFGS: 23 01:41:10 -342.890206* 0.0118
LBFGS: 24 01:41:10 -342.890263* 0.0100
LBFGS: 25 01:41:10 -342.890296* 0.0103
LBFGS: 26 01:41:10 -342.890354* 0.0125
LBFGS: 27 01:41:10 -342.890372* 0.0095
LBFGS: 28 01:41:10 -342.890409* 0.0074
LBFGS: 29 01:41:10 -342.890423* 0.0077
LBFGS: 30 01:41:10 -342.890432* 0.0084
LBFGS: 31 01:41:10 -342.890466* 0.0071
LBFGS: 32 01:41:11 -342.890492* 0.0061
LBFGS: 33 01:41:11 -342.890501* 0.0057
LBFGS: 34 01:41:11 -342.890509* 0.0057
LBFGS: 35 01:41:11 -342.890536* 0.0056
LBFGS: 36 01:41:12 -342.890545* 0.0053
LBFGS: 37 01:41:12 -342.890546* 0.0050
LBFGS: 38 01:41:12 -342.890562* 0.0048
LBFGS: 39 01:41:12 -342.890560* 0.0048
LBFGS: 40 01:41:12 -342.890576* 0.0044
LBFGS: 41 01:41:13 -342.890581* 0.0044
LBFGS: 42 01:41:13 -342.890600* 0.0047
LBFGS: 43 01:41:13 -342.890603* 0.0056
LBFGS: 44 01:41:13 -342.890611* 0.0047
LBFGS: 45 01:41:13 -342.890614* 0.0037
LBFGS: 46 01:41:13 -342.890630* 0.0037
LBFGS: 47 01:41:13 -342.890633* 0.0051
LBFGS: 48 01:41:13 -342.890622* 0.0047
LBFGS: 49 01:41:13 -342.890640* 0.0034
LBFGS: 50 01:41:14 -342.890634* 0.0034
LBFGS: 51 01:41:14 -342.890646* 0.0033
LBFGS: 52 01:41:14 -342.890657* 0.0030
LBFGS: 53 01:41:14 -342.890640* 0.0020
LBFGS: 54 01:41:14 -342.890636* 0.0020
LBFGS: 55 01:41:14 -342.890657* 0.0019
LBFGS: 56 01:41:14 -342.890672* 0.0017
LBFGS: 57 01:41:14 -342.890672* 0.0010
CPU times: user 502 ms, sys: 104 ms, total: 605 ms
Wall time: 11 s
[11]:
E_adsorp = E_slab_mol - (E_slab + E_mol)
print(f"E_adsorp {E_adsorp}, E_slab_mol {E_slab_mol}, E_slab {E_slab}, E_mol {E_mol}")
E_adsorp -1.8752501080074353, E_slab_mol -342.8906722017586, E_slab -329.51522411375305, E_mol -11.500197979998095
An adsorption energy of -1.87 eV was obtained.
[12]:
from pfcc_extras.visualize.view import view_ngl
view_ngl(slab_mol)
Adsorption sites¶
You can specify the adsorption site name instead of the coordinate value as the fourth positional argument to the add_adsorbate method.
[13]:
%%time
slab_mol = fcc111("Pt", a=a, size=(4, 4, 4), vacuum=40.0)
slab_mol.pbc = True
mol2 = molecule("CO")
add_adsorbate(slab_mol, mol2, 3.0, "ontop")
E_slab_mol = get_opt_energy(slab_mol, fmax=1e-4)
Step Time Energy fmax
*Force-consistent energies used in optimization.
LBFGS: 0 01:41:23 -342.400027* 0.9021
LBFGS: 1 01:41:24 -342.463759* 1.6972
LBFGS: 2 01:41:24 -342.594994* 0.4611
LBFGS: 3 01:41:24 -342.609518* 0.3743
LBFGS: 4 01:41:24 -342.655741* 0.6197
LBFGS: 5 01:41:25 -342.671914* 0.6355
LBFGS: 6 01:41:25 -342.685582* 0.4110
LBFGS: 7 01:41:25 -342.695878* 0.1851
LBFGS: 8 01:41:25 -342.707594* 0.3033
LBFGS: 9 01:41:25 -342.720403* 0.3934
LBFGS: 10 01:41:26 -342.732799* 0.2824
LBFGS: 11 01:41:26 -342.739539* 0.1139
LBFGS: 12 01:41:26 -342.744451* 0.1459
LBFGS: 13 01:41:26 -342.748675* 0.2545
LBFGS: 14 01:41:26 -342.751861* 0.2209
LBFGS: 15 01:41:26 -342.754075* 0.0854
LBFGS: 16 01:41:26 -342.755366* 0.0532
LBFGS: 17 01:41:27 -342.756432* 0.1189
LBFGS: 18 01:41:27 -342.757729* 0.1454
LBFGS: 19 01:41:27 -342.758905* 0.0954
LBFGS: 20 01:41:28 -342.759695* 0.0364
LBFGS: 21 01:41:28 -342.760252* 0.0439
LBFGS: 22 01:41:28 -342.760728* 0.0692
LBFGS: 23 01:41:28 -342.761134* 0.0537
LBFGS: 24 01:41:28 -342.761417* 0.0288
LBFGS: 25 01:41:28 -342.761598* 0.0229
LBFGS: 26 01:41:28 -342.761711* 0.0317
LBFGS: 27 01:41:28 -342.761821* 0.0234
LBFGS: 28 01:41:29 -342.761879* 0.0108
LBFGS: 29 01:41:29 -342.761931* 0.0117
LBFGS: 30 01:41:29 -342.761969* 0.0178
LBFGS: 31 01:41:29 -342.761998* 0.0124
LBFGS: 32 01:41:29 -342.762003* 0.0050
LBFGS: 33 01:41:29 -342.762012* 0.0045
LBFGS: 34 01:41:29 -342.762025* 0.0062
LBFGS: 35 01:41:29 -342.762016* 0.0037
LBFGS: 36 01:41:29 -342.762029* 0.0019
LBFGS: 37 01:41:29 -342.762035* 0.0028
LBFGS: 38 01:41:29 -342.762029* 0.0036
LBFGS: 39 01:41:30 -342.762037* 0.0022
LBFGS: 40 01:41:30 -342.762044* 0.0014
LBFGS: 41 01:41:30 -342.762035* 0.0012
LBFGS: 42 01:41:31 -342.762038* 0.0016
LBFGS: 43 01:41:31 -342.762050* 0.0010
LBFGS: 44 01:41:32 -342.762044* 0.0007
LBFGS: 45 01:41:32 -342.762032* 0.0009
LBFGS: 46 01:41:32 -342.762049* 0.0016
LBFGS: 47 01:41:33 -342.762051* 0.0015
LBFGS: 48 01:41:33 -342.762035* 0.0007
LBFGS: 49 01:41:33 -342.762033* 0.0004
LBFGS: 50 01:41:33 -342.762023* 0.0005
LBFGS: 51 01:41:33 -342.762046* 0.0006
LBFGS: 52 01:41:33 -342.762032* 0.0003
LBFGS: 53 01:41:33 -342.762058* 0.0002
LBFGS: 54 01:41:34 -342.762039* 0.0003
LBFGS: 55 01:41:34 -342.762018* 0.0003
LBFGS: 56 01:41:34 -342.762024* 0.0002
LBFGS: 57 01:41:34 -342.762020* 0.0001
CPU times: user 503 ms, sys: 99.9 ms, total: 603 ms
Wall time: 10.9 s
The name of this adsorption site can be found by checking the “sites” in the .info attribute.
[14]:
slab_mol.info
[14]:
{'adsorbate_info': {'cell': array([[2.80539597, 0. ],
[1.40269799, 2.42954418]]),
'sites': {'ontop': (0, 0),
'bridge': (0.5, 0),
'fcc': (0.3333333333333333, 0.3333333333333333),
'hcp': (0.6666666666666666, 0.6666666666666666)},
'top layer atom index': 48}}
The “fcc111”, in this example, has ontop, bridge, fcc, and hcp sites Let’s take a look at each of them.
[15]:
orthogonal = False
slab_mol_ontop = fcc111("Pt", a=a, size=(4, 4, 4), vacuum=40.0, orthogonal=orthogonal, periodic=True)
mol = molecule("CO")
add_adsorbate(slab_mol_ontop, mol, 3.0, "ontop")
slab_mol_bridge = fcc111("Pt", a=a, size=(4, 4, 4), vacuum=40.0, orthogonal=orthogonal, periodic=True)
mol = molecule("CO")
add_adsorbate(slab_mol_bridge, mol, 3.0, "bridge")
slab_mol_fcc = fcc111("Pt", a=a, size=(4, 4, 4), vacuum=40.0, orthogonal=orthogonal, periodic=True)
mol = molecule("CO")
add_adsorbate(slab_mol_fcc, mol, 3.0, "fcc")
slab_mol_hcp = fcc111("Pt", a=a, size=(4, 4, 4), vacuum=40.0, orthogonal=orthogonal, periodic=True)
mol = molecule("CO")
add_adsorbate(slab_mol_hcp, mol, 3.0, "hcp")
[16]:
view_ngl([slab_mol_ontop, slab_mol_bridge, slab_mol_fcc, slab_mol_hcp])
[17]:
from ase.io import write
from IPython.display import Image
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
write("output/pt_ontop.png", slab_mol_ontop, rotation="0x,0y,90z")
write("output/pt_bridge.png", slab_mol_bridge, rotation="0x,0y,90z")
write("output/pt_fcc.png", slab_mol_fcc, rotation="0x,0y,90z")
write("output/pt_hcp.png", slab_mol_hcp, rotation="0x,0y,90z")
fig, axes = plt.subplots(1, 4, figsize=(12, 3))
for i, adsorbate_name in enumerate(["ontop", "bridge", "fcc", "hcp"]):
ax = axes[i]
ax.imshow(mpimg.imread(f"output/pt_{adsorbate_name}.png"))
ax.set_axis_off()
ax.set_title(adsorbate_name)
fig.show()
As we can see above, each adsorption site is as follows
ontop: on one atombridge: between two atomsfcc: between three atoms with no atoms in the second layerhcp: between three atoms, directly above the atom in the second layer
The adsorption site of a molecule on a surface structure is non-trivial, depending on the slab and molecular system. It is necessary to find the location where the adsorption energy is minimium (minimum energy of the adsorbed structure).
[18]:
E_slab_mol_ontop = get_opt_energy(slab_mol_ontop, fmax=1e-3)
E_slab_mol_bridge = get_opt_energy(slab_mol_bridge, fmax=1e-3)
E_slab_mol_fcc = get_opt_energy(slab_mol_fcc, fmax=1e-3)
E_slab_mol_hcp = get_opt_energy(slab_mol_hcp, fmax=1e-3)
Step Time Energy fmax
*Force-consistent energies used in optimization.
LBFGS: 0 01:42:00 -342.400062* 0.9021
LBFGS: 1 01:42:01 -342.463767* 1.6972
LBFGS: 2 01:42:01 -342.594999* 0.4611
LBFGS: 3 01:42:01 -342.609523* 0.3743
LBFGS: 4 01:42:01 -342.655739* 0.6197
LBFGS: 5 01:42:01 -342.671902* 0.6354
LBFGS: 6 01:42:01 -342.685580* 0.4110
LBFGS: 7 01:42:02 -342.695879* 0.1851
LBFGS: 8 01:42:02 -342.707599* 0.3033
LBFGS: 9 01:42:02 -342.720400* 0.3934
LBFGS: 10 01:42:02 -342.732795* 0.2824
LBFGS: 11 01:42:02 -342.739533* 0.1139
LBFGS: 12 01:42:03 -342.744444* 0.1459
LBFGS: 13 01:42:03 -342.748672* 0.2545
LBFGS: 14 01:42:03 -342.751865* 0.2209
LBFGS: 15 01:42:03 -342.754070* 0.0853
LBFGS: 16 01:42:03 -342.755365* 0.0532
LBFGS: 17 01:42:04 -342.756441* 0.1189
LBFGS: 18 01:42:04 -342.757735* 0.1454
LBFGS: 19 01:42:04 -342.758905* 0.0954
LBFGS: 20 01:42:04 -342.759706* 0.0364
LBFGS: 21 01:42:04 -342.760252* 0.0439
LBFGS: 22 01:42:04 -342.760727* 0.0692
LBFGS: 23 01:42:05 -342.761138* 0.0537
LBFGS: 24 01:42:05 -342.761416* 0.0288
LBFGS: 25 01:42:05 -342.761603* 0.0229
LBFGS: 26 01:42:05 -342.761705* 0.0317
LBFGS: 27 01:42:05 -342.761823* 0.0234
LBFGS: 28 01:42:05 -342.761882* 0.0108
LBFGS: 29 01:42:05 -342.761928* 0.0117
LBFGS: 30 01:42:05 -342.761966* 0.0178
LBFGS: 31 01:42:06 -342.761998* 0.0124
LBFGS: 32 01:42:06 -342.762004* 0.0050
LBFGS: 33 01:42:06 -342.762015* 0.0045
LBFGS: 34 01:42:06 -342.762030* 0.0062
LBFGS: 35 01:42:06 -342.762016* 0.0037
LBFGS: 36 01:42:06 -342.762031* 0.0019
LBFGS: 37 01:42:07 -342.762037* 0.0028
LBFGS: 38 01:42:07 -342.762038* 0.0036
LBFGS: 39 01:42:08 -342.762024* 0.0022
LBFGS: 40 01:42:08 -342.762052* 0.0014
LBFGS: 41 01:42:08 -342.762044* 0.0012
LBFGS: 42 01:42:08 -342.762027* 0.0016
LBFGS: 43 01:42:09 -342.762055* 0.0010
Step Time Energy fmax
*Force-consistent energies used in optimization.
LBFGS: 0 01:42:09 -341.994158* 3.7071
LBFGS: 1 01:42:10 -342.131338* 4.4263
LBFGS: 2 01:42:10 -342.428966* 3.3710
LBFGS: 3 01:42:11 -342.724250* 0.9431
LBFGS: 4 01:42:11 -342.799982* 0.7936
LBFGS: 5 01:42:11 -342.828830* 0.4370
LBFGS: 6 01:42:11 -342.872847* 0.2373
LBFGS: 7 01:42:11 -342.876431* 0.1118
LBFGS: 8 01:42:11 -342.880149* 0.0750
LBFGS: 9 01:42:11 -342.883222* 0.0957
LBFGS: 10 01:42:12 -342.885122* 0.1009
LBFGS: 11 01:42:12 -342.886202* 0.0659
LBFGS: 12 01:42:12 -342.887140* 0.0422
LBFGS: 13 01:42:12 -342.887979* 0.0357
LBFGS: 14 01:42:12 -342.888544* 0.0436
LBFGS: 15 01:42:12 -342.888836* 0.0318
LBFGS: 16 01:42:13 -342.889117* 0.0229
LBFGS: 17 01:42:13 -342.889394* 0.0221
LBFGS: 18 01:42:13 -342.889629* 0.0267
LBFGS: 19 01:42:13 -342.889779* 0.0217
LBFGS: 20 01:42:13 -342.889895* 0.0178
LBFGS: 21 01:42:13 -342.890030* 0.0165
LBFGS: 22 01:42:14 -342.890136* 0.0169
LBFGS: 23 01:42:14 -342.890211* 0.0118
LBFGS: 24 01:42:14 -342.890256* 0.0100
LBFGS: 25 01:42:14 -342.890293* 0.0103
LBFGS: 26 01:42:14 -342.890358* 0.0125
LBFGS: 27 01:42:14 -342.890379* 0.0095
LBFGS: 28 01:42:14 -342.890403* 0.0074
LBFGS: 29 01:42:15 -342.890409* 0.0077
LBFGS: 30 01:42:15 -342.890440* 0.0084
LBFGS: 31 01:42:15 -342.890469* 0.0071
LBFGS: 32 01:42:15 -342.890489* 0.0061
LBFGS: 33 01:42:15 -342.890500* 0.0057
LBFGS: 34 01:42:15 -342.890516* 0.0057
LBFGS: 35 01:42:15 -342.890528* 0.0056
LBFGS: 36 01:42:15 -342.890544* 0.0053
LBFGS: 37 01:42:15 -342.890540* 0.0050
LBFGS: 38 01:42:15 -342.890559* 0.0048
LBFGS: 39 01:42:16 -342.890561* 0.0048
LBFGS: 40 01:42:16 -342.890571* 0.0044
LBFGS: 41 01:42:16 -342.890576* 0.0044
LBFGS: 42 01:42:16 -342.890606* 0.0047
LBFGS: 43 01:42:16 -342.890602* 0.0056
LBFGS: 44 01:42:16 -342.890612* 0.0047
LBFGS: 45 01:42:16 -342.890608* 0.0037
LBFGS: 46 01:42:16 -342.890629* 0.0037
LBFGS: 47 01:42:17 -342.890630* 0.0051
LBFGS: 48 01:42:17 -342.890624* 0.0047
LBFGS: 49 01:42:18 -342.890642* 0.0034
LBFGS: 50 01:42:18 -342.890638* 0.0034
LBFGS: 51 01:42:18 -342.890640* 0.0033
LBFGS: 52 01:42:19 -342.890654* 0.0030
LBFGS: 53 01:42:20 -342.890638* 0.0020
LBFGS: 54 01:42:20 -342.890631* 0.0020
LBFGS: 55 01:42:20 -342.890653* 0.0019
LBFGS: 56 01:42:20 -342.890665* 0.0017
LBFGS: 57 01:42:21 -342.890658* 0.0010
Step Time Energy fmax
*Force-consistent energies used in optimization.
LBFGS: 0 01:42:21 -341.816981* 3.6377
LBFGS: 1 01:42:21 -341.917168* 4.3622
LBFGS: 2 01:42:21 -342.143828* 3.0800
LBFGS: 3 01:42:21 -342.591136* 2.0656
LBFGS: 4 01:42:21 -342.676549* 1.2542
LBFGS: 5 01:42:21 -342.853229* 0.3963
LBFGS: 6 01:42:22 -342.861437* 0.2819
LBFGS: 7 01:42:22 -342.875649* 0.2979
LBFGS: 8 01:42:22 -342.883638* 0.2914
LBFGS: 9 01:42:22 -342.887090* 0.1377
LBFGS: 10 01:42:22 -342.890143* 0.1399
LBFGS: 11 01:42:22 -342.894157* 0.2395
LBFGS: 12 01:42:22 -342.897206* 0.2040
LBFGS: 13 01:42:23 -342.898653* 0.0849
LBFGS: 14 01:42:23 -342.899499* 0.0606
LBFGS: 15 01:42:23 -342.900390* 0.1245
LBFGS: 16 01:42:24 -342.901232* 0.1252
LBFGS: 17 01:42:24 -342.901736* 0.0635
LBFGS: 18 01:42:24 -342.901959* 0.0226
LBFGS: 19 01:42:24 -342.902119* 0.0516
LBFGS: 20 01:42:24 -342.902340* 0.0632
LBFGS: 21 01:42:24 -342.902449* 0.0408
LBFGS: 22 01:42:24 -342.902552* 0.0104
LBFGS: 23 01:42:24 -342.902567* 0.0148
LBFGS: 24 01:42:25 -342.902597* 0.0215
LBFGS: 25 01:42:25 -342.902622* 0.0159
LBFGS: 26 01:42:25 -342.902628* 0.0051
LBFGS: 27 01:42:25 -342.902660* 0.0088
LBFGS: 28 01:42:25 -342.902662* 0.0139
LBFGS: 29 01:42:25 -342.902672* 0.0115
LBFGS: 30 01:42:25 -342.902705* 0.0045
LBFGS: 31 01:42:26 -342.902683* 0.0060
LBFGS: 32 01:42:26 -342.902707* 0.0117
LBFGS: 33 01:42:26 -342.902702* 0.0114
LBFGS: 34 01:42:26 -342.902726* 0.0054
LBFGS: 35 01:42:26 -342.902690* 0.0022
LBFGS: 36 01:42:26 -342.902709* 0.0050
LBFGS: 37 01:42:26 -342.902719* 0.0051
LBFGS: 38 01:42:27 -342.902718* 0.0024
LBFGS: 39 01:42:27 -342.902727* 0.0007
Step Time Energy fmax
*Force-consistent energies used in optimization.
LBFGS: 0 01:42:27 -341.805857* 3.6304
LBFGS: 1 01:42:28 -341.905943* 4.3609
LBFGS: 2 01:42:28 -342.132912* 3.0710
LBFGS: 3 01:42:28 -342.575720* 2.0315
LBFGS: 4 01:42:28 -342.659006* 1.2385
LBFGS: 5 01:42:28 -342.830603* 0.3875
LBFGS: 6 01:42:29 -342.838912* 0.2924
LBFGS: 7 01:42:29 -342.850330* 0.2383
LBFGS: 8 01:42:29 -342.861207* 0.2834
LBFGS: 9 01:42:30 -342.864594* 0.1600
LBFGS: 10 01:42:30 -342.868136* 0.1236
LBFGS: 11 01:42:30 -342.872275* 0.2191
LBFGS: 12 01:42:30 -342.875798* 0.1969
LBFGS: 13 01:42:30 -342.877524* 0.0883
LBFGS: 14 01:42:31 -342.878510* 0.0848
LBFGS: 15 01:42:32 -342.879528* 0.1112
LBFGS: 16 01:42:32 -342.880491* 0.1151
LBFGS: 17 01:42:32 -342.881073* 0.0594
LBFGS: 18 01:42:32 -342.881386* 0.0223
LBFGS: 19 01:42:32 -342.881593* 0.0460
LBFGS: 20 01:42:32 -342.881821* 0.0534
LBFGS: 21 01:42:32 -342.881985* 0.0320
LBFGS: 22 01:42:33 -342.882093* 0.0147
LBFGS: 23 01:42:33 -342.882149* 0.0212
LBFGS: 24 01:42:33 -342.882255* 0.0235
LBFGS: 25 01:42:33 -342.882266* 0.0130
LBFGS: 26 01:42:34 -342.882320* 0.0085
LBFGS: 27 01:42:34 -342.882356* 0.0165
LBFGS: 28 01:42:34 -342.882358* 0.0183
LBFGS: 29 01:42:34 -342.882381* 0.0101
LBFGS: 30 01:42:34 -342.882421* 0.0077
LBFGS: 31 01:42:34 -342.882435* 0.0116
LBFGS: 32 01:42:34 -342.882475* 0.0149
LBFGS: 33 01:42:34 -342.882459* 0.0098
LBFGS: 34 01:42:35 -342.882441* 0.0045
LBFGS: 35 01:42:35 -342.882455* 0.0048
LBFGS: 36 01:42:35 -342.882484* 0.0062
LBFGS: 37 01:42:35 -342.882504* 0.0035
LBFGS: 38 01:42:35 -342.882483* 0.0018
LBFGS: 39 01:42:35 -342.882487* 0.0022
LBFGS: 40 01:42:35 -342.882498* 0.0032
LBFGS: 41 01:42:35 -342.882489* 0.0021
LBFGS: 42 01:42:35 -342.882478* 0.0008
[19]:
import pandas as pd
series = pd.Series({
"ontop": E_slab_mol_ontop,
"bridge": E_slab_mol_bridge,
"fcc": E_slab_mol_fcc,
"hcp": E_slab_mol_hcp,
})
E_adsorp_series = series - (E_slab + E_mol)
E_adsorp_series.plot(style="o-", xlabel="Adsorption site", ylabel="Energy [eV]", title="Adsorption energy for each adsorption site")
plt.show()
We see that the energy of the adsorption structure differs depending on the adsorption site. In this calculation, the fcc hollow site is the most stable and has the smallest adsorption energy.
The results can be compared with the calculations performed using the PBE functional described in Table 1. and Table S2. in the following reference. Adsorption energies are reported as 1.80, 1.88 and 1.89 eV for the top, bridge, and hcp hollow sites, respectively. Note that the top site is known to be the most stable in experiments, and calculations performed using PBE functional do not match with the experimental results.
[20]:
view_ngl([slab_mol_ontop, slab_mol_bridge, slab_mol_fcc, slab_mol_hcp])
FixAtoms constrains¶
All atomic coordinates were relaxed in the above calculation. But when using DFT to obtain adsorption structures, the atomic coordinates of the lower layer of slab may be fixed during structural optimization to increase the stability of the structure.
Let’s try this method by using the ASE constraint FixAtoms, we can fix the atomic coordinates at a given index.
[21]:
%%time
slab_mol_fixed = fcc111("Pt", a=a, size=(4, 4, 4), vacuum=40.0)
slab_mol_fixed.pbc = True
mol2 = molecule("CO")
add_adsorbate(slab_mol_fixed, mol2, 3.0, "bridge")
view_ngl(slab_mol_fixed)
CPU times: user 219 ms, sys: 11 ms, total: 230 ms
Wall time: 236 ms
The threshold of the z coordinate is used to specify the atoms in the bottom two layers.
You can see the z coordinate from the above visualization by hovering the cursor over it, but this time, let’s visualize and check the z-axis distribution as follows.
[22]:
import pandas as pd
atoms = slab_mol_fixed
df = pd.DataFrame({
"x": atoms.positions[:, 0],
"y": atoms.positions[:, 1],
"z": atoms.positions[:, 2],
"symbol": atoms.symbols,
})
df
[22]:
| x | y | z | symbol | |
|---|---|---|---|---|
| 0 | 0.000000 | 0.000000 | 40.000000 | Pt |
| 1 | 2.805396 | 0.000000 | 40.000000 | Pt |
| 2 | 5.610792 | 0.000000 | 40.000000 | Pt |
| 3 | 8.416188 | 0.000000 | 40.000000 | Pt |
| 4 | 1.402698 | 2.429544 | 40.000000 | Pt |
| ... | ... | ... | ... | ... |
| 61 | 7.013490 | 7.288633 | 46.871789 | Pt |
| 62 | 9.818886 | 7.288633 | 46.871789 | Pt |
| 63 | 12.624282 | 7.288633 | 46.871789 | Pt |
| 64 | 1.402698 | 0.000000 | 49.871789 | O |
| 65 | 1.402698 | 0.000000 | 48.721449 | C |
66 rows × 4 columns
[23]:
import plotly.express as px
coord = "z"
df_sorted = df.sort_values(coord).reset_index().rename({"index": "atom_index"}, axis=1)
fig = px.scatter(df_sorted, x=df_sorted.index, y=coord, color="symbol", hover_data=["x", "y", "z", "atom_index"])
fig.show()
The boundary between the second and third layers was found to be between 42.29 - 44.58.
The FixAtoms can be used by either setting an index of atoms to fix in the indices or setting an array of bool to fix or not to fix each atom in the mask.
We will use mask in this example.
[24]:
from ase.constraints import FixAtoms
thresh = 43.0
constraint = FixAtoms(mask=slab_mol_fixed.positions[:, 2] < thresh)
constraint
[24]:
FixAtoms(indices=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31])
The created constraint can be set with the atoms.set_constraint method.
[25]:
slab_mol_fixed.set_constraint(constraint)
Visualize and check that the intended location is fixed.
[26]:
v = view_ngl(slab_mol_fixed)
v.show_fix_atoms_constraint()
v
[27]:
E_slab_mol_fixed = get_opt_energy(slab_mol_fixed, fmax=1e-3)
Step Time Energy fmax
*Force-consistent energies used in optimization.
LBFGS: 0 01:48:07 -341.994172* 3.7071
LBFGS: 1 01:48:08 -342.107705* 4.4264
LBFGS: 2 01:48:09 -342.373263* 3.3174
LBFGS: 3 01:48:09 -342.579870* 2.0002
LBFGS: 4 01:48:09 -342.756605* 0.8708
LBFGS: 5 01:48:09 -342.790083* 0.3747
LBFGS: 6 01:48:09 -342.809386* 0.2911
LBFGS: 7 01:48:09 -342.823247* 0.2666
LBFGS: 8 01:48:10 -342.833935* 0.1241
LBFGS: 9 01:48:10 -342.835557* 0.0883
LBFGS: 10 01:48:10 -342.836829* 0.0753
LBFGS: 11 01:48:10 -342.838312* 0.0687
LBFGS: 12 01:48:10 -342.839152* 0.0649
LBFGS: 13 01:48:10 -342.839639* 0.0642
LBFGS: 14 01:48:11 -342.840021* 0.0514
LBFGS: 15 01:48:11 -342.840415* 0.0501
LBFGS: 16 01:48:11 -342.840788* 0.0328
LBFGS: 17 01:48:12 -342.841010* 0.0303
LBFGS: 18 01:48:12 -342.841144* 0.0279
LBFGS: 19 01:48:12 -342.841276* 0.0326
LBFGS: 20 01:48:12 -342.841395* 0.0268
LBFGS: 21 01:48:13 -342.841455* 0.0133
LBFGS: 22 01:48:13 -342.841516* 0.0147
LBFGS: 23 01:48:13 -342.841542* 0.0196
LBFGS: 24 01:48:13 -342.841588* 0.0174
LBFGS: 25 01:48:13 -342.841639* 0.0091
LBFGS: 26 01:48:13 -342.841665* 0.0096
LBFGS: 27 01:48:13 -342.841685* 0.0117
LBFGS: 28 01:48:13 -342.841731* 0.0141
LBFGS: 29 01:48:14 -342.841752* 0.0098
LBFGS: 30 01:48:15 -342.841774* 0.0066
LBFGS: 31 01:48:16 -342.841795* 0.0080
LBFGS: 32 01:48:16 -342.841815* 0.0105
LBFGS: 33 01:48:16 -342.841819* 0.0098
LBFGS: 34 01:48:16 -342.841840* 0.0072
LBFGS: 35 01:48:16 -342.841853* 0.0062
LBFGS: 36 01:48:16 -342.841868* 0.0096
LBFGS: 37 01:48:17 -342.841881* 0.0105
LBFGS: 38 01:48:17 -342.841906* 0.0075
LBFGS: 39 01:48:17 -342.841916* 0.0067
LBFGS: 40 01:48:18 -342.841925* 0.0086
LBFGS: 41 01:48:18 -342.841943* 0.0107
LBFGS: 42 01:48:18 -342.841955* 0.0084
LBFGS: 43 01:48:18 -342.841981* 0.0054
LBFGS: 44 01:48:18 -342.841993* 0.0048
LBFGS: 45 01:48:18 -342.841983* 0.0062
LBFGS: 46 01:48:18 -342.841989* 0.0056
LBFGS: 47 01:48:19 -342.841993* 0.0026
LBFGS: 48 01:48:19 -342.842015* 0.0022
LBFGS: 49 01:48:20 -342.842003* 0.0026
LBFGS: 50 01:48:20 -342.842009* 0.0026
LBFGS: 51 01:48:20 -342.842011* 0.0016
LBFGS: 52 01:48:21 -342.842007* 0.0011
LBFGS: 53 01:48:21 -342.841999* 0.0009
[28]:
E_adsorp = E_slab_mol_fixed - (E_slab + E_mol)
print(f"E_adsorp {E_adsorp}, E_slab_mol {E_slab_mol_fixed}, E_slab {E_slab}, E_mol {E_mol}")
E_adsorp -1.8265767995379179, E_slab_mol -342.8419988932891, E_slab -329.51522411375305, E_mol -11.500197979998095
When the lower layer is fixed, an adsorption energy of -1.83 eV was obtained.
Reference
CO on Pd¶
Let’s calculate CO adsorption on Pd in the same manner.
[29]:
bulk_atoms = bulk("Pd", cubic=True)
np.mean(np.diag(bulk_atoms.cell))
[29]:
3.89
[30]:
bulk_atoms.calc = calculator
E_bulk = get_opt_energy(bulk_atoms, fmax=1e-4, opt_mode="scale")
a_pd = np.mean(np.diag(bulk_atoms.cell))
a_pd
Step Time Energy fmax
*Force-consistent energies used in optimization.
LBFGS: 0 01:48:21 -14.872185* 4.6900
LBFGS: 1 01:48:21 -14.761184* 7.1380
LBFGS: 2 01:48:22 -14.923393* 0.7071
LBFGS: 3 01:48:23 -14.924605* 0.1525
LBFGS: 4 01:48:23 -14.924674* 0.0011
LBFGS: 5 01:48:23 -14.924674* 0.0000
[30]:
3.9414398359309444
[31]:
from ase.build import fcc111, molecule, add_adsorbate
slab = fcc111("Pd", a=a_pd, size=(4, 4, 4), vacuum=40.0, periodic=True)
slab.calc = calculator
E_slab = get_opt_energy(slab, fmax=1e-4, opt_mode="normal")
Step Time Energy fmax
*Force-consistent energies used in optimization.
LBFGS: 0 01:48:23 -220.564914* 0.1150
LBFGS: 1 01:48:24 -220.571361* 0.0928
LBFGS: 2 01:48:24 -220.582676* 0.0110
LBFGS: 3 01:48:24 -220.582765* 0.0013
LBFGS: 4 01:48:24 -220.582755* 0.0001
[32]:
mol = molecule("CO")
mol.calc = calculator
E_mol = get_opt_energy(mol, fmax=1e-4)
Step Time Energy fmax
*Force-consistent energies used in optimization.
LBFGS: 0 01:48:25 -11.497237* 0.8202
LBFGS: 1 01:48:25 -11.484526* 1.9590
LBFGS: 2 01:48:25 -11.500196* 0.0298
LBFGS: 3 01:48:26 -11.500203* 0.0011
LBFGS: 4 01:48:26 -11.500200* 0.0000
[33]:
%%time
slab_mol = fcc111("Pd", a=a_pd, size=(4, 4, 4), vacuum=40.0)
slab_mol.pbc = True
mol2 = molecule("CO")
add_adsorbate(slab_mol, mol2, 3.0, "ontop")
E_slab_mol = get_opt_energy(slab_mol, fmax=1e-3)
Step Time Energy fmax
*Force-consistent energies used in optimization.
LBFGS: 0 01:48:26 -233.487975* 1.0792
LBFGS: 1 01:48:26 -233.491409* 1.8944
LBFGS: 2 01:48:26 -233.515945* 0.6717
LBFGS: 3 01:48:26 -233.535180* 0.5258
LBFGS: 4 01:48:27 -233.540240* 0.3398
LBFGS: 5 01:48:27 -233.543963* 0.1764
LBFGS: 6 01:48:27 -233.548140* 0.3449
LBFGS: 7 01:48:27 -233.552462* 0.4205
LBFGS: 8 01:48:27 -233.555429* 0.2889
LBFGS: 9 01:48:28 -233.557796* 0.1417
LBFGS: 10 01:48:28 -233.560578* 0.2719
LBFGS: 11 01:48:28 -233.564315* 0.3132
LBFGS: 12 01:48:29 -233.568304* 0.2848
LBFGS: 13 01:48:29 -233.571393* 0.1970
LBFGS: 14 01:48:29 -233.573304* 0.1138
LBFGS: 15 01:48:30 -233.574655* 0.1490
LBFGS: 16 01:48:30 -233.575624* 0.1528
LBFGS: 17 01:48:30 -233.576389* 0.0986
LBFGS: 18 01:48:30 -233.576915* 0.0838
LBFGS: 19 01:48:31 -233.577257* 0.0699
LBFGS: 20 01:48:32 -233.577498* 0.0479
LBFGS: 21 01:48:32 -233.577660* 0.0408
LBFGS: 22 01:48:32 -233.577832* 0.0678
LBFGS: 23 01:48:32 -233.577997* 0.0695
LBFGS: 24 01:48:32 -233.578120* 0.0350
LBFGS: 25 01:48:32 -233.578186* 0.0196
LBFGS: 26 01:48:32 -233.578239* 0.0378
LBFGS: 27 01:48:33 -233.578351* 0.0555
LBFGS: 28 01:48:33 -233.578392* 0.0412
LBFGS: 29 01:48:33 -233.578434* 0.0142
LBFGS: 30 01:48:33 -233.578474* 0.0100
LBFGS: 31 01:48:34 -233.578455* 0.0204
LBFGS: 32 01:48:35 -233.578489* 0.0195
LBFGS: 33 01:48:35 -233.578505* 0.0084
LBFGS: 34 01:48:36 -233.578475* 0.0024
LBFGS: 35 01:48:36 -233.578514* 0.0080
LBFGS: 36 01:48:37 -233.578505* 0.0085
LBFGS: 37 01:48:37 -233.578498* 0.0040
LBFGS: 38 01:48:38 -233.578482* 0.0009
CPU times: user 369 ms, sys: 103 ms, total: 472 ms
Wall time: 11.9 s
[34]:
E_adsorp = E_slab + E_mol - E_slab_mol
print(f"E_adsorp {E_adsorp}, E_slab_mol {E_slab_mol}, E_slab {E_slab}, E_mol {E_mol}")
E_adsorp 1.4955274062837418, E_slab_mol -233.57848222102731, E_slab -220.58275492892844, E_mol -11.50019988581515
[35]:
from pfcc_extras.visualize.view import view_ngl
view_ngl(slab_mol)
In this way, adsorption energies can be calculated for various surface and molecular combinations. This is an essential process in the search for efficient catalysts, which will be discussed in chapter 5.
Coverage ratio¶
In the example discussed in this tutorial, only one molecule is adsorbed on a surface. In reality, there are cases where many molecules react with the surface at the same time, where multiple molecules are adsorbed at the same time, or where the surface is almost completely covered with adsorbed molecules.
In this case, we are interested in the adsorption energy of a single molecule adsorbed on a slab where multiple molecules are already adsorbed.
The analysis of the coverage ratio dependency of the adsorption energy of CO molecules on the Pd surface is reported in the following literature
Adsorption to multiple locations¶
Another application topic not covered here is adsorption to more than one location. When molecules are large and complex, there are some cases that become stable by adsorbing to multiple locations on the surface.