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_U0, model_version="v3.0.0")
calculator = ASECalculator(estimator)
pfp_api_client: 1.18.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
LBFGS: 0 14:21:21 -21.749341 6.135578
LBFGS: 1 14:21:21 -21.296372 15.553629
LBFGS: 2 14:21:21 -21.809601 1.021754
/tmp/ipykernel_2402/1510506507.py:7: DeprecationWarning: Please use atoms.calc = calc
atoms.set_calculator(calculator)
/tmp/ipykernel_2402/1510506507.py:9: FutureWarning: Import StrainFilter from ase.filters
opt1 = LBFGS(StrainFilter(atoms, mask=[1, 1, 1, 0, 0, 0]))
LBFGS: 3 14:21:21 -21.811500 0.156843
LBFGS: 4 14:21:21 -21.811549 0.002061
LBFGS: 5 14:21:21 -21.811545 0.000003
[4]:
-21.81154493913019
[5]:
bulk_atoms.cell
[5]:
Cell([3.967428746092142, 3.967428368434641, 3.9674287321083397])
[6]:
import numpy as np
a = np.mean(np.diag(bulk_atoms.cell))
a
[6]:
3.9674286155450407
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
LBFGS: 0 14:21:21 -329.391240 0.322358
LBFGS: 1 14:21:21 -329.444738 0.243357
LBFGS: 2 14:21:21 -329.511075 0.043367
/tmp/ipykernel_2402/1510506507.py:7: DeprecationWarning: Please use atoms.calc = calc
atoms.set_calculator(calculator)
LBFGS: 3 14:21:21 -329.512257 0.022308
LBFGS: 4 14:21:21 -329.512616 0.022152
LBFGS: 5 14:21:21 -329.514956 0.010743
LBFGS: 6 14:21:21 -329.515157 0.005390
LBFGS: 7 14:21:21 -329.515235 0.000673
LBFGS: 8 14:21:22 -329.515122 0.000039
[9]:
mol = molecule("CO")
mol.calc = calculator
E_mol = get_opt_energy(mol, fmax=1e-4)
Step Time Energy fmax
LBFGS: 0 14:21:22 -11.497233 0.820161
LBFGS: 1 14:21:22 -11.484526 1.958979
LBFGS: 2 14:21:22 -11.500194 0.029776
/tmp/ipykernel_2402/1510506507.py:7: DeprecationWarning: Please use atoms.calc = calc
atoms.set_calculator(calculator)
LBFGS: 3 14:21:22 -11.500199 0.001079
LBFGS: 4 14:21:22 -11.500198 0.000010
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
LBFGS: 0 14:21:22 -341.994067 3.707105
LBFGS: 1 14:21:22 -342.131356 4.426341
/tmp/ipykernel_2402/1510506507.py:7: DeprecationWarning: Please use atoms.calc = calc
atoms.set_calculator(calculator)
LBFGS: 2 14:21:22 -342.429018 3.370975
LBFGS: 3 14:21:22 -342.724332 0.943148
LBFGS: 4 14:21:22 -342.800055 0.793638
LBFGS: 5 14:21:22 -342.828916 0.437042
LBFGS: 6 14:21:22 -342.872909 0.237324
LBFGS: 7 14:21:22 -342.876486 0.111755
LBFGS: 8 14:21:22 -342.880222 0.075047
LBFGS: 9 14:21:23 -342.883305 0.095742
LBFGS: 10 14:21:23 -342.885177 0.100864
LBFGS: 11 14:21:23 -342.886265 0.065901
LBFGS: 12 14:21:23 -342.887208 0.042197
LBFGS: 13 14:21:23 -342.888047 0.035733
LBFGS: 14 14:21:23 -342.888613 0.043612
LBFGS: 15 14:21:23 -342.888868 0.031783
LBFGS: 16 14:21:23 -342.889168 0.022924
LBFGS: 17 14:21:23 -342.889450 0.022106
LBFGS: 18 14:21:23 -342.889721 0.026692
LBFGS: 19 14:21:23 -342.889832 0.021688
LBFGS: 20 14:21:23 -342.889968 0.017798
LBFGS: 21 14:21:23 -342.890106 0.016531
LBFGS: 22 14:21:23 -342.890210 0.016840
LBFGS: 23 14:21:23 -342.890299 0.011784
LBFGS: 24 14:21:24 -342.890323 0.010000
LBFGS: 25 14:21:24 -342.890392 0.010334
LBFGS: 26 14:21:24 -342.890425 0.012523
LBFGS: 27 14:21:24 -342.890447 0.009516
LBFGS: 28 14:21:24 -342.890467 0.007353
LBFGS: 29 14:21:24 -342.890488 0.007703
LBFGS: 30 14:21:24 -342.890494 0.008397
LBFGS: 31 14:21:24 -342.890531 0.007112
LBFGS: 32 14:21:24 -342.890548 0.006135
LBFGS: 33 14:21:24 -342.890565 0.005694
LBFGS: 34 14:21:24 -342.890575 0.005681
LBFGS: 35 14:21:24 -342.890603 0.005636
LBFGS: 36 14:21:24 -342.890626 0.005332
LBFGS: 37 14:21:24 -342.890614 0.005024
LBFGS: 38 14:21:24 -342.890641 0.004830
LBFGS: 39 14:21:24 -342.890622 0.004810
LBFGS: 40 14:21:25 -342.890665 0.004385
LBFGS: 41 14:21:25 -342.890639 0.004397
LBFGS: 42 14:21:25 -342.890684 0.004667
LBFGS: 43 14:21:25 -342.890663 0.005615
LBFGS: 44 14:21:25 -342.890678 0.004676
LBFGS: 45 14:21:25 -342.890695 0.003660
LBFGS: 46 14:21:25 -342.890719 0.003655
LBFGS: 47 14:21:25 -342.890701 0.005097
LBFGS: 48 14:21:25 -342.890696 0.004733
LBFGS: 49 14:21:25 -342.890708 0.003382
LBFGS: 50 14:21:25 -342.890722 0.003401
LBFGS: 51 14:21:25 -342.890728 0.003294
LBFGS: 52 14:21:25 -342.890725 0.003050
LBFGS: 53 14:21:25 -342.890738 0.002019
LBFGS: 54 14:21:25 -342.890747 0.001996
LBFGS: 55 14:21:25 -342.890733 0.001843
LBFGS: 56 14:21:26 -342.890749 0.001675
LBFGS: 57 14:21:26 -342.890751 0.001001
LBFGS: 58 14:21:26 -342.890732 0.000938
CPU times: user 347 ms, sys: 49 ms, total: 396 ms
Wall time: 3.82 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.8754126953839432, E_slab_mol -342.89073213009203, E_slab -329.5151218083796, E_mol -11.500197626328491
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
LBFGS: 0 14:21:33 -342.399955 0.902052
LBFGS: 1 14:21:33 -342.463820 1.697183
LBFGS: 2 14:21:33 -342.595051 0.461154
/tmp/ipykernel_2402/1510506507.py:7: DeprecationWarning: Please use atoms.calc = calc
atoms.set_calculator(calculator)
LBFGS: 3 14:21:33 -342.609561 0.374312
LBFGS: 4 14:21:33 -342.655787 0.619668
LBFGS: 5 14:21:33 -342.671953 0.635467
LBFGS: 6 14:21:33 -342.685593 0.410986
LBFGS: 7 14:21:33 -342.695922 0.185128
LBFGS: 8 14:21:33 -342.707634 0.303415
LBFGS: 9 14:21:33 -342.720433 0.393416
LBFGS: 10 14:21:34 -342.732847 0.282414
LBFGS: 11 14:21:34 -342.739607 0.113855
LBFGS: 12 14:21:34 -342.744493 0.145870
LBFGS: 13 14:21:34 -342.748729 0.254530
LBFGS: 14 14:21:34 -342.751903 0.220945
LBFGS: 15 14:21:34 -342.754126 0.085359
LBFGS: 16 14:21:34 -342.755420 0.053167
LBFGS: 17 14:21:34 -342.756482 0.118871
LBFGS: 18 14:21:34 -342.757791 0.145426
LBFGS: 19 14:21:34 -342.758975 0.095405
LBFGS: 20 14:21:34 -342.759756 0.036352
LBFGS: 21 14:21:34 -342.760324 0.043871
LBFGS: 22 14:21:34 -342.760776 0.069186
LBFGS: 23 14:21:34 -342.761192 0.053748
LBFGS: 24 14:21:35 -342.761441 0.028788
LBFGS: 25 14:21:35 -342.761644 0.022931
LBFGS: 26 14:21:35 -342.761787 0.031662
LBFGS: 27 14:21:35 -342.761886 0.023382
LBFGS: 28 14:21:35 -342.761962 0.010754
LBFGS: 29 14:21:35 -342.762006 0.011718
LBFGS: 30 14:21:35 -342.762019 0.017848
LBFGS: 31 14:21:35 -342.762079 0.012378
LBFGS: 32 14:21:35 -342.762086 0.005001
LBFGS: 33 14:21:35 -342.762066 0.004476
LBFGS: 34 14:21:35 -342.762081 0.006191
LBFGS: 35 14:21:35 -342.762064 0.003673
LBFGS: 36 14:21:35 -342.762064 0.001884
LBFGS: 37 14:21:35 -342.762081 0.002779
LBFGS: 38 14:21:35 -342.762090 0.003624
LBFGS: 39 14:21:35 -342.762081 0.002235
LBFGS: 40 14:21:35 -342.762095 0.001411
LBFGS: 41 14:21:36 -342.762083 0.001223
LBFGS: 42 14:21:36 -342.762076 0.001633
LBFGS: 43 14:21:36 -342.762115 0.000993
LBFGS: 44 14:21:36 -342.762071 0.000690
LBFGS: 45 14:21:36 -342.762092 0.000912
LBFGS: 46 14:21:36 -342.762073 0.001545
LBFGS: 47 14:21:36 -342.762066 0.001473
LBFGS: 48 14:21:36 -342.762067 0.000673
LBFGS: 49 14:21:36 -342.762086 0.000399
LBFGS: 50 14:21:36 -342.762065 0.000534
LBFGS: 51 14:21:36 -342.762078 0.000613
LBFGS: 52 14:21:36 -342.762057 0.000294
LBFGS: 53 14:21:36 -342.762074 0.000169
LBFGS: 54 14:21:36 -342.762085 0.000286
LBFGS: 55 14:21:36 -342.762057 0.000305
LBFGS: 56 14:21:36 -342.762042 0.000144
LBFGS: 57 14:21:37 -342.762050 0.000072
CPU times: user 342 ms, sys: 64.6 ms, total: 406 ms
Wall time: 3.63 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.80539568, 0. ],
[1.40269784, 2.42954392]]),
'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
LBFGS: 0 14:21:38 -342.399929 0.902054
LBFGS: 1 14:21:38 -342.463815 1.697184
LBFGS: 2 14:21:38 -342.595049 0.461122
/tmp/ipykernel_2402/1510506507.py:7: DeprecationWarning: Please use atoms.calc = calc
atoms.set_calculator(calculator)
LBFGS: 3 14:21:39 -342.609558 0.374303
LBFGS: 4 14:21:39 -342.655783 0.619656
LBFGS: 5 14:21:39 -342.671954 0.635447
LBFGS: 6 14:21:39 -342.685590 0.410967
LBFGS: 7 14:21:39 -342.695941 0.185127
LBFGS: 8 14:21:39 -342.707641 0.303349
LBFGS: 9 14:21:39 -342.720469 0.393403
LBFGS: 10 14:21:39 -342.732855 0.282399
LBFGS: 11 14:21:39 -342.739609 0.113855
LBFGS: 12 14:21:39 -342.744494 0.145903
LBFGS: 13 14:21:39 -342.748741 0.254528
LBFGS: 14 14:21:39 -342.751915 0.220938
LBFGS: 15 14:21:39 -342.754123 0.085354
LBFGS: 16 14:21:39 -342.755443 0.053166
LBFGS: 17 14:21:39 -342.756487 0.118871
LBFGS: 18 14:21:39 -342.757801 0.145417
LBFGS: 19 14:21:39 -342.758962 0.095414
LBFGS: 20 14:21:40 -342.759785 0.036355
LBFGS: 21 14:21:40 -342.760316 0.043894
LBFGS: 22 14:21:40 -342.760790 0.069189
LBFGS: 23 14:21:40 -342.761194 0.053739
LBFGS: 24 14:21:40 -342.761435 0.028788
LBFGS: 25 14:21:40 -342.761633 0.022932
LBFGS: 26 14:21:40 -342.761793 0.031664
LBFGS: 27 14:21:40 -342.761891 0.023381
LBFGS: 28 14:21:40 -342.761960 0.010756
LBFGS: 29 14:21:40 -342.762002 0.011715
LBFGS: 30 14:21:40 -342.762004 0.017841
LBFGS: 31 14:21:40 -342.762063 0.012370
LBFGS: 32 14:21:40 -342.762090 0.005001
LBFGS: 33 14:21:40 -342.762061 0.004469
LBFGS: 34 14:21:40 -342.762082 0.006188
LBFGS: 35 14:21:40 -342.762069 0.003673
LBFGS: 36 14:21:40 -342.762063 0.001883
LBFGS: 37 14:21:40 -342.762081 0.002777
LBFGS: 38 14:21:41 -342.762091 0.003629
LBFGS: 39 14:21:41 -342.762059 0.002238
LBFGS: 40 14:21:41 -342.762091 0.001410
LBFGS: 41 14:21:41 -342.762085 0.001217
LBFGS: 42 14:21:41 -342.762082 0.001625
LBFGS: 43 14:21:41 -342.762108 0.001000
LBFGS: 44 14:21:41 -342.762070 0.000691
Step Time Energy fmax
LBFGS: 0 14:21:41 -341.994068 3.707104
LBFGS: 1 14:21:42 -342.131361 4.426339
LBFGS: 2 14:21:42 -342.429013 3.370969
LBFGS: 3 14:21:42 -342.724333 0.943137
LBFGS: 4 14:21:42 -342.800057 0.793617
LBFGS: 5 14:21:42 -342.828916 0.437051
LBFGS: 6 14:21:42 -342.872930 0.237322
LBFGS: 7 14:21:42 -342.876489 0.111753
LBFGS: 8 14:21:42 -342.880220 0.075050
LBFGS: 9 14:21:42 -342.883296 0.095742
LBFGS: 10 14:21:42 -342.885185 0.100861
LBFGS: 11 14:21:42 -342.886274 0.065890
LBFGS: 12 14:21:42 -342.887206 0.042198
LBFGS: 13 14:21:42 -342.888036 0.035742
LBFGS: 14 14:21:42 -342.888611 0.043613
LBFGS: 15 14:21:42 -342.888880 0.031777
LBFGS: 16 14:21:42 -342.889164 0.022923
LBFGS: 17 14:21:42 -342.889449 0.022100
LBFGS: 18 14:21:42 -342.889717 0.026692
LBFGS: 19 14:21:43 -342.889824 0.021687
LBFGS: 20 14:21:43 -342.889971 0.017798
LBFGS: 21 14:21:43 -342.890121 0.016483
LBFGS: 22 14:21:43 -342.890207 0.016868
LBFGS: 23 14:21:43 -342.890301 0.011804
LBFGS: 24 14:21:43 -342.890314 0.009989
LBFGS: 25 14:21:43 -342.890384 0.010333
LBFGS: 26 14:21:43 -342.890434 0.012522
LBFGS: 27 14:21:43 -342.890448 0.009516
LBFGS: 28 14:21:43 -342.890462 0.007353
LBFGS: 29 14:21:43 -342.890484 0.007708
LBFGS: 30 14:21:43 -342.890493 0.008407
LBFGS: 31 14:21:43 -342.890530 0.007109
LBFGS: 32 14:21:43 -342.890519 0.006141
LBFGS: 33 14:21:43 -342.890566 0.005695
LBFGS: 34 14:21:43 -342.890566 0.005676
LBFGS: 35 14:21:43 -342.890602 0.005644
LBFGS: 36 14:21:44 -342.890621 0.005324
LBFGS: 37 14:21:44 -342.890595 0.005026
LBFGS: 38 14:21:44 -342.890627 0.004829
LBFGS: 39 14:21:44 -342.890629 0.004812
LBFGS: 40 14:21:44 -342.890648 0.004384
LBFGS: 41 14:21:44 -342.890642 0.004406
LBFGS: 42 14:21:44 -342.890680 0.004658
LBFGS: 43 14:21:44 -342.890655 0.005631
LBFGS: 44 14:21:44 -342.890677 0.004679
LBFGS: 45 14:21:44 -342.890685 0.003663
LBFGS: 46 14:21:44 -342.890701 0.003660
LBFGS: 47 14:21:44 -342.890699 0.005087
LBFGS: 48 14:21:44 -342.890715 0.004735
LBFGS: 49 14:21:44 -342.890732 0.003381
LBFGS: 50 14:21:44 -342.890705 0.003406
LBFGS: 51 14:21:44 -342.890719 0.003280
LBFGS: 52 14:21:45 -342.890731 0.003047
LBFGS: 53 14:21:45 -342.890739 0.002019
LBFGS: 54 14:21:45 -342.890735 0.001995
LBFGS: 55 14:21:45 -342.890736 0.001850
LBFGS: 56 14:21:45 -342.890743 0.001676
LBFGS: 57 14:21:45 -342.890749 0.000997
Step Time Energy fmax
LBFGS: 0 14:21:45 -341.816890 3.637717
LBFGS: 1 14:21:45 -341.917226 4.362195
LBFGS: 2 14:21:45 -342.143871 3.080002
LBFGS: 3 14:21:45 -342.591186 2.065524
LBFGS: 4 14:21:45 -342.676575 1.254153
LBFGS: 5 14:21:45 -342.853280 0.396113
LBFGS: 6 14:21:45 -342.861514 0.281871
LBFGS: 7 14:21:45 -342.875709 0.297696
LBFGS: 8 14:21:45 -342.883728 0.291208
LBFGS: 9 14:21:45 -342.887138 0.137633
LBFGS: 10 14:21:46 -342.890184 0.139875
LBFGS: 11 14:21:46 -342.894233 0.239341
LBFGS: 12 14:21:46 -342.897241 0.203828
LBFGS: 13 14:21:46 -342.898725 0.084824
LBFGS: 14 14:21:46 -342.899589 0.060598
LBFGS: 15 14:21:46 -342.900452 0.124498
LBFGS: 16 14:21:46 -342.901272 0.125059
LBFGS: 17 14:21:46 -342.901798 0.063416
LBFGS: 18 14:21:46 -342.902036 0.022631
LBFGS: 19 14:21:46 -342.902185 0.051523
LBFGS: 20 14:21:46 -342.902402 0.063152
LBFGS: 21 14:21:46 -342.902525 0.040769
LBFGS: 22 14:21:46 -342.902591 0.010405
LBFGS: 23 14:21:46 -342.902639 0.014826
LBFGS: 24 14:21:46 -342.902679 0.021462
LBFGS: 25 14:21:46 -342.902700 0.015898
LBFGS: 26 14:21:46 -342.902718 0.005060
LBFGS: 27 14:21:47 -342.902699 0.008795
LBFGS: 28 14:21:47 -342.902717 0.013945
LBFGS: 29 14:21:47 -342.902768 0.011543
LBFGS: 30 14:21:47 -342.902767 0.004522
LBFGS: 31 14:21:47 -342.902777 0.005967
LBFGS: 32 14:21:47 -342.902800 0.011696
LBFGS: 33 14:21:47 -342.902785 0.011374
LBFGS: 34 14:21:47 -342.902802 0.005357
LBFGS: 35 14:21:47 -342.902779 0.002168
LBFGS: 36 14:21:47 -342.902777 0.004994
LBFGS: 37 14:21:47 -342.902793 0.005139
LBFGS: 38 14:21:47 -342.902792 0.002448
LBFGS: 39 14:21:47 -342.902766 0.000696
Step Time Energy fmax
LBFGS: 0 14:21:47 -341.805788 3.630381
LBFGS: 1 14:21:47 -341.906010 4.360870
LBFGS: 2 14:21:47 -342.133001 3.070985
LBFGS: 3 14:21:48 -342.575791 2.031607
LBFGS: 4 14:21:48 -342.659042 1.238475
LBFGS: 5 14:21:48 -342.830708 0.387585
LBFGS: 6 14:21:48 -342.838979 0.292372
LBFGS: 7 14:21:48 -342.850407 0.238356
LBFGS: 8 14:21:48 -342.861280 0.283465
LBFGS: 9 14:21:48 -342.864647 0.159965
LBFGS: 10 14:21:48 -342.868204 0.123570
LBFGS: 11 14:21:48 -342.872314 0.219109
LBFGS: 12 14:21:48 -342.875857 0.196918
LBFGS: 13 14:21:48 -342.877613 0.088332
LBFGS: 14 14:21:48 -342.878571 0.084757
LBFGS: 15 14:21:48 -342.879536 0.111184
LBFGS: 16 14:21:48 -342.880518 0.115103
LBFGS: 17 14:21:48 -342.881125 0.059387
LBFGS: 18 14:21:48 -342.881471 0.022328
LBFGS: 19 14:21:49 -342.881656 0.046029
LBFGS: 20 14:21:49 -342.881875 0.053349
LBFGS: 21 14:21:49 -342.882041 0.032034
LBFGS: 22 14:21:49 -342.882160 0.014692
LBFGS: 23 14:21:49 -342.882236 0.021242
LBFGS: 24 14:21:49 -342.882321 0.023529
LBFGS: 25 14:21:49 -342.882347 0.013020
LBFGS: 26 14:21:49 -342.882364 0.008465
LBFGS: 27 14:21:49 -342.882424 0.016525
LBFGS: 28 14:21:49 -342.882419 0.018357
LBFGS: 29 14:21:49 -342.882435 0.010142
LBFGS: 30 14:21:49 -342.882489 0.007714
LBFGS: 31 14:21:49 -342.882508 0.011590
LBFGS: 32 14:21:49 -342.882547 0.014856
LBFGS: 33 14:21:49 -342.882513 0.009807
LBFGS: 34 14:21:49 -342.882499 0.004526
LBFGS: 35 14:21:49 -342.882504 0.004809
LBFGS: 36 14:21:50 -342.882515 0.006189
LBFGS: 37 14:21:50 -342.882532 0.003508
LBFGS: 38 14:21:50 -342.882564 0.001783
LBFGS: 39 14:21:50 -342.882545 0.002227
LBFGS: 40 14:21:50 -342.882540 0.003167
LBFGS: 41 14:21:50 -342.882549 0.002123
LBFGS: 42 14:21:50 -342.882524 0.000812
[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 97.7 ms, sys: 7.32 ms, total: 105 ms
Wall time: 112 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.610791 | 0.000000 | 40.000000 | Pt |
| 3 | 8.416187 | 0.000000 | 40.000000 | Pt |
| 4 | 1.402698 | 2.429544 | 40.000000 | Pt |
| ... | ... | ... | ... | ... |
| 61 | 7.013489 | 7.288632 | 46.871788 | Pt |
| 62 | 9.818885 | 7.288632 | 46.871788 | Pt |
| 63 | 12.624281 | 7.288632 | 46.871788 | Pt |
| 64 | 1.402698 | 0.000000 | 49.871788 | O |
| 65 | 1.402698 | 0.000000 | 48.721448 | 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
LBFGS: 0 14:21:52 -341.994066 3.707103
LBFGS: 1 14:21:52 -342.107708 4.426348
LBFGS: 2 14:21:52 -342.373294 3.317411
/tmp/ipykernel_2402/1510506507.py:7: DeprecationWarning:
Please use atoms.calc = calc
LBFGS: 3 14:21:52 -342.579885 2.000330
LBFGS: 4 14:21:52 -342.756645 0.870795
LBFGS: 5 14:21:52 -342.790102 0.374663
LBFGS: 6 14:21:52 -342.809421 0.291134
LBFGS: 7 14:21:52 -342.823274 0.266582
LBFGS: 8 14:21:52 -342.833961 0.124068
LBFGS: 9 14:21:52 -342.835576 0.088252
LBFGS: 10 14:21:53 -342.836858 0.075297
LBFGS: 11 14:21:53 -342.838391 0.068676
LBFGS: 12 14:21:53 -342.839206 0.064896
LBFGS: 13 14:21:53 -342.839665 0.064187
LBFGS: 14 14:21:53 -342.840034 0.051451
LBFGS: 15 14:21:53 -342.840451 0.050059
LBFGS: 16 14:21:53 -342.840821 0.032788
LBFGS: 17 14:21:53 -342.841032 0.030302
LBFGS: 18 14:21:53 -342.841183 0.027928
LBFGS: 19 14:21:53 -342.841316 0.032598
LBFGS: 20 14:21:53 -342.841461 0.026782
LBFGS: 21 14:21:53 -342.841503 0.013326
LBFGS: 22 14:21:53 -342.841527 0.014671
LBFGS: 23 14:21:53 -342.841584 0.019564
LBFGS: 24 14:21:53 -342.841645 0.017394
LBFGS: 25 14:21:54 -342.841674 0.009111
LBFGS: 26 14:21:54 -342.841700 0.009561
LBFGS: 27 14:21:54 -342.841723 0.011716
LBFGS: 28 14:21:54 -342.841777 0.014141
LBFGS: 29 14:21:54 -342.841805 0.009851
LBFGS: 30 14:21:54 -342.841836 0.006582
LBFGS: 31 14:21:54 -342.841845 0.008036
LBFGS: 32 14:21:54 -342.841855 0.010523
LBFGS: 33 14:21:54 -342.841858 0.009743
LBFGS: 34 14:21:54 -342.841871 0.007184
LBFGS: 35 14:21:54 -342.841879 0.006200
LBFGS: 36 14:21:54 -342.841899 0.009582
LBFGS: 37 14:21:54 -342.841923 0.010455
LBFGS: 38 14:21:54 -342.841943 0.007517
LBFGS: 39 14:21:54 -342.841954 0.006664
LBFGS: 40 14:21:54 -342.841970 0.008600
LBFGS: 41 14:21:54 -342.841977 0.010671
LBFGS: 42 14:21:54 -342.842022 0.008441
LBFGS: 43 14:21:55 -342.842018 0.005408
LBFGS: 44 14:21:55 -342.842035 0.004768
LBFGS: 45 14:21:55 -342.842032 0.006172
LBFGS: 46 14:21:55 -342.842023 0.005571
LBFGS: 47 14:21:55 -342.842038 0.002580
LBFGS: 48 14:21:55 -342.842036 0.002194
LBFGS: 49 14:21:55 -342.842026 0.002596
LBFGS: 50 14:21:55 -342.842015 0.002593
LBFGS: 51 14:21:55 -342.842041 0.001618
LBFGS: 52 14:21:55 -342.842062 0.001117
LBFGS: 53 14:21:55 -342.842030 0.000947
[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.8267101219207689, E_slab_mol -342.84202955662886, E_slab -329.5151218083796, E_mol -11.500197626328491
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
LBFGS: 0 14:21:55 -14.872188 4.690002
LBFGS: 1 14:21:55 -14.761178 7.138005
LBFGS: 2 14:21:55 -14.923395 0.707045
/tmp/ipykernel_2402/1510506507.py:7: DeprecationWarning:
Please use atoms.calc = calc
/tmp/ipykernel_2402/1510506507.py:9: FutureWarning:
Import StrainFilter from ase.filters
LBFGS: 3 14:21:55 -14.924614 0.152437
LBFGS: 4 14:21:55 -14.924663 0.001121
LBFGS: 5 14:21:56 -14.924665 0.000015
[30]:
3.9414396825202953
[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
LBFGS: 0 14:21:56 -220.564855 0.115019
LBFGS: 1 14:21:56 -220.571326 0.092846
LBFGS: 2 14:21:56 -220.582642 0.010979
/tmp/ipykernel_2402/1510506507.py:7: DeprecationWarning:
Please use atoms.calc = calc
LBFGS: 3 14:21:56 -220.582741 0.001269
LBFGS: 4 14:21:56 -220.582703 0.000075
[32]:
mol = molecule("CO")
mol.calc = calculator
E_mol = get_opt_energy(mol, fmax=1e-4)
Step Time Energy fmax
LBFGS: 0 14:21:56 -11.497233 0.820161
LBFGS: 1 14:21:56 -11.484526 1.958979
LBFGS: 2 14:21:56 -11.500195 0.029777
LBFGS: 3 14:21:56 -11.500199 0.001079
/tmp/ipykernel_2402/1510506507.py:7: DeprecationWarning:
Please use atoms.calc = calc
LBFGS: 4 14:21:56 -11.500198 0.000009
[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
LBFGS: 0 14:21:56 -233.487893 1.079224
LBFGS: 1 14:21:56 -233.491378 1.894393
LBFGS: 2 14:21:56 -233.515958 0.671717
/tmp/ipykernel_2402/1510506507.py:7: DeprecationWarning:
Please use atoms.calc = calc
LBFGS: 3 14:21:56 -233.535149 0.525830
LBFGS: 4 14:21:56 -233.540274 0.339791
LBFGS: 5 14:21:56 -233.543944 0.176446
LBFGS: 6 14:21:57 -233.548148 0.344857
LBFGS: 7 14:21:57 -233.552480 0.420469
LBFGS: 8 14:21:57 -233.555402 0.288860
LBFGS: 9 14:21:57 -233.557787 0.141675
LBFGS: 10 14:21:57 -233.560567 0.271889
LBFGS: 11 14:21:57 -233.564279 0.313259
LBFGS: 12 14:21:57 -233.568334 0.284809
LBFGS: 13 14:21:57 -233.571354 0.197038
LBFGS: 14 14:21:57 -233.573287 0.113800
LBFGS: 15 14:21:57 -233.574607 0.148964
LBFGS: 16 14:21:57 -233.575602 0.152801
LBFGS: 17 14:21:57 -233.576374 0.098614
LBFGS: 18 14:21:57 -233.576939 0.083779
LBFGS: 19 14:21:57 -233.577267 0.069898
LBFGS: 20 14:21:57 -233.577457 0.047900
LBFGS: 21 14:21:58 -233.577647 0.040828
LBFGS: 22 14:21:58 -233.577807 0.067762
LBFGS: 23 14:21:58 -233.578013 0.069454
LBFGS: 24 14:21:58 -233.578091 0.034968
LBFGS: 25 14:21:58 -233.578182 0.019564
LBFGS: 26 14:21:58 -233.578254 0.037787
LBFGS: 27 14:21:58 -233.578334 0.055525
LBFGS: 28 14:21:58 -233.578413 0.041198
LBFGS: 29 14:21:58 -233.578473 0.014187
LBFGS: 30 14:21:58 -233.578457 0.009982
LBFGS: 31 14:21:58 -233.578497 0.020393
LBFGS: 32 14:21:58 -233.578480 0.019467
LBFGS: 33 14:21:58 -233.578490 0.008406
LBFGS: 34 14:21:58 -233.578493 0.002450
LBFGS: 35 14:21:58 -233.578532 0.008077
LBFGS: 36 14:21:58 -233.578517 0.008449
LBFGS: 37 14:21:59 -233.578513 0.003999
LBFGS: 38 14:21:59 -233.578514 0.000911
CPU times: user 243 ms, sys: 27.3 ms, total: 271 ms
Wall time: 2.45 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.4956127981674285, E_slab_mol -233.5785141151784, E_slab -220.5827033337869, E_mol -11.500197983224071
[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.