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.19.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 07:54:49 -21.749341 6.135577
LBFGS: 1 07:54:49 -21.296372 15.553638
/tmp/ipykernel_8561/1510506507.py:7: DeprecationWarning: Please use atoms.calc = calc
atoms.set_calculator(calculator)
/tmp/ipykernel_8561/1510506507.py:9: FutureWarning: Import StrainFilter from ase.filters
opt1 = LBFGS(StrainFilter(atoms, mask=[1, 1, 1, 0, 0, 0]))
LBFGS: 2 07:54:49 -21.809599 1.021796
LBFGS: 3 07:54:49 -21.811497 0.156920
LBFGS: 4 07:54:49 -21.811544 0.002066
LBFGS: 5 07:54:49 -21.811544 0.000010
[4]:
-21.811543944246235
[5]:
bulk_atoms.cell
[5]:
Cell([3.9674289205307174, 3.967428532149144, 3.967428615755973])
[6]:
import numpy as np
a = np.mean(np.diag(bulk_atoms.cell))
a
[6]:
3.967428689478611
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 07:54:49 -329.391234 0.322358
LBFGS: 1 07:54:49 -329.444709 0.243355
/tmp/ipykernel_8561/1510506507.py:7: DeprecationWarning: Please use atoms.calc = calc
atoms.set_calculator(calculator)
LBFGS: 2 07:54:50 -329.511078 0.043360
LBFGS: 3 07:54:50 -329.512193 0.022307
LBFGS: 4 07:54:50 -329.512651 0.022152
LBFGS: 5 07:54:50 -329.514980 0.010735
LBFGS: 6 07:54:50 -329.515248 0.005379
LBFGS: 7 07:54:50 -329.515195 0.000678
LBFGS: 8 07:54:50 -329.515152 0.000035
[9]:
mol = molecule("CO")
mol.calc = calculator
E_mol = get_opt_energy(mol, fmax=1e-4)
Step Time Energy fmax
LBFGS: 0 07:54:50 -11.497234 0.820163
LBFGS: 1 07:54:50 -11.484524 1.958990
LBFGS: 2 07:54:50 -11.500194 0.029776
/tmp/ipykernel_8561/1510506507.py:7: DeprecationWarning: Please use atoms.calc = calc
atoms.set_calculator(calculator)
LBFGS: 3 07:54:50 -11.500197 0.001082
LBFGS: 4 07:54:50 -11.500199 0.000007
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 07:54:50 -341.994046 3.707106
LBFGS: 1 07:54:51 -342.131358 4.426330
/tmp/ipykernel_8561/1510506507.py:7: DeprecationWarning: Please use atoms.calc = calc
atoms.set_calculator(calculator)
LBFGS: 2 07:54:51 -342.429021 3.370984
LBFGS: 3 07:54:51 -342.724384 0.943089
LBFGS: 4 07:54:51 -342.800038 0.793622
LBFGS: 5 07:54:51 -342.828912 0.437039
LBFGS: 6 07:54:51 -342.872906 0.237297
LBFGS: 7 07:54:51 -342.876484 0.111745
LBFGS: 8 07:54:51 -342.880215 0.075044
LBFGS: 9 07:54:51 -342.883296 0.095750
LBFGS: 10 07:54:51 -342.885186 0.100860
LBFGS: 11 07:54:51 -342.886260 0.065886
LBFGS: 12 07:54:51 -342.887206 0.042198
LBFGS: 13 07:54:51 -342.888046 0.035738
LBFGS: 14 07:54:51 -342.888617 0.043613
LBFGS: 15 07:54:52 -342.888869 0.031778
LBFGS: 16 07:54:52 -342.889159 0.022924
LBFGS: 17 07:54:52 -342.889444 0.022108
LBFGS: 18 07:54:52 -342.889710 0.026686
LBFGS: 19 07:54:52 -342.889830 0.021686
LBFGS: 20 07:54:52 -342.889972 0.017799
LBFGS: 21 07:54:52 -342.890118 0.016525
LBFGS: 22 07:54:52 -342.890202 0.016837
LBFGS: 23 07:54:52 -342.890292 0.011782
LBFGS: 24 07:54:52 -342.890313 0.010002
LBFGS: 25 07:54:52 -342.890382 0.010341
LBFGS: 26 07:54:52 -342.890426 0.012521
LBFGS: 27 07:54:53 -342.890447 0.009512
LBFGS: 28 07:54:53 -342.890484 0.007368
LBFGS: 29 07:54:53 -342.890478 0.007693
LBFGS: 30 07:54:53 -342.890493 0.008405
LBFGS: 31 07:54:53 -342.890537 0.007140
LBFGS: 32 07:54:53 -342.890515 0.006148
LBFGS: 33 07:54:53 -342.890559 0.005694
LBFGS: 34 07:54:53 -342.890563 0.005676
LBFGS: 35 07:54:53 -342.890609 0.005642
LBFGS: 36 07:54:53 -342.890605 0.005324
LBFGS: 37 07:54:53 -342.890591 0.005027
LBFGS: 38 07:54:53 -342.890633 0.004830
LBFGS: 39 07:54:53 -342.890661 0.004811
LBFGS: 40 07:54:54 -342.890648 0.004381
LBFGS: 41 07:54:54 -342.890633 0.004414
LBFGS: 42 07:54:54 -342.890678 0.004662
LBFGS: 43 07:54:54 -342.890645 0.005618
LBFGS: 44 07:54:54 -342.890670 0.004685
LBFGS: 45 07:54:54 -342.890705 0.003669
LBFGS: 46 07:54:54 -342.890699 0.003658
LBFGS: 47 07:54:54 -342.890704 0.005095
LBFGS: 48 07:54:54 -342.890697 0.004729
LBFGS: 49 07:54:54 -342.890707 0.003384
LBFGS: 50 07:54:54 -342.890703 0.003407
LBFGS: 51 07:54:54 -342.890725 0.003284
LBFGS: 52 07:54:54 -342.890736 0.003033
LBFGS: 53 07:54:54 -342.890737 0.002020
LBFGS: 54 07:54:55 -342.890725 0.001996
LBFGS: 55 07:54:55 -342.890712 0.001847
LBFGS: 56 07:54:55 -342.890735 0.001679
LBFGS: 57 07:54:55 -342.890737 0.000997
CPU times: user 377 ms, sys: 31.7 ms, total: 408 ms
Wall time: 4.37 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.8753863462617915, E_slab_mol -342.89073734826064, E_slab -329.5151524229501, E_mol -11.500198579048762
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 07:55:06 -342.399910 0.902054
LBFGS: 1 07:55:06 -342.463817 1.697185
/tmp/ipykernel_8561/1510506507.py:7: DeprecationWarning: Please use atoms.calc = calc
atoms.set_calculator(calculator)
LBFGS: 2 07:55:06 -342.595038 0.461127
LBFGS: 3 07:55:06 -342.609559 0.374306
LBFGS: 4 07:55:06 -342.655779 0.619674
LBFGS: 5 07:55:06 -342.671951 0.635450
LBFGS: 6 07:55:06 -342.685604 0.410959
LBFGS: 7 07:55:06 -342.695938 0.185124
LBFGS: 8 07:55:06 -342.707631 0.303364
LBFGS: 9 07:55:06 -342.720453 0.393410
LBFGS: 10 07:55:06 -342.732852 0.282378
LBFGS: 11 07:55:07 -342.739609 0.113854
LBFGS: 12 07:55:07 -342.744490 0.145889
LBFGS: 13 07:55:07 -342.748745 0.254531
LBFGS: 14 07:55:07 -342.751913 0.220942
LBFGS: 15 07:55:07 -342.754120 0.085357
LBFGS: 16 07:55:07 -342.755419 0.053168
LBFGS: 17 07:55:07 -342.756486 0.118864
LBFGS: 18 07:55:07 -342.757801 0.145427
LBFGS: 19 07:55:07 -342.758954 0.095398
LBFGS: 20 07:55:07 -342.759748 0.036356
LBFGS: 21 07:55:07 -342.760303 0.043884
LBFGS: 22 07:55:08 -342.760789 0.069188
LBFGS: 23 07:55:08 -342.761185 0.053738
LBFGS: 24 07:55:08 -342.761434 0.028786
LBFGS: 25 07:55:08 -342.761635 0.022932
LBFGS: 26 07:55:08 -342.761780 0.031662
LBFGS: 27 07:55:08 -342.761886 0.023381
LBFGS: 28 07:55:08 -342.761952 0.010756
LBFGS: 29 07:55:08 -342.762002 0.011717
LBFGS: 30 07:55:08 -342.762022 0.017857
LBFGS: 31 07:55:08 -342.762073 0.012374
LBFGS: 32 07:55:08 -342.762078 0.005001
LBFGS: 33 07:55:08 -342.762066 0.004472
LBFGS: 34 07:55:08 -342.762079 0.006191
LBFGS: 35 07:55:09 -342.762057 0.003670
LBFGS: 36 07:55:09 -342.762067 0.001881
LBFGS: 37 07:55:09 -342.762076 0.002791
LBFGS: 38 07:55:09 -342.762088 0.003626
LBFGS: 39 07:55:09 -342.762083 0.002242
LBFGS: 40 07:55:09 -342.762098 0.001410
LBFGS: 41 07:55:09 -342.762075 0.001218
LBFGS: 42 07:55:09 -342.762077 0.001632
LBFGS: 43 07:55:09 -342.762120 0.001008
LBFGS: 44 07:55:09 -342.762066 0.000690
LBFGS: 45 07:55:09 -342.762096 0.000912
LBFGS: 46 07:55:10 -342.762103 0.001544
LBFGS: 47 07:55:10 -342.762073 0.001479
LBFGS: 48 07:55:10 -342.762093 0.000675
LBFGS: 49 07:55:10 -342.762060 0.000407
LBFGS: 50 07:55:10 -342.762065 0.000547
LBFGS: 51 07:55:10 -342.762073 0.000606
LBFGS: 52 07:55:10 -342.762052 0.000299
LBFGS: 53 07:55:10 -342.762062 0.000176
LBFGS: 54 07:55:10 -342.762070 0.000289
LBFGS: 55 07:55:10 -342.762041 0.000301
LBFGS: 56 07:55:10 -342.762044 0.000155
LBFGS: 57 07:55:10 -342.762075 0.000072
CPU times: user 381 ms, sys: 47.8 ms, total: 429 ms
Wall time: 4.87 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.80539573, 0. ],
[1.40269787, 2.42954397]]),
'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 07:55:12 -342.399913 0.902056
LBFGS: 1 07:55:12 -342.463811 1.697190
LBFGS: 2 07:55:12 -342.595038 0.461123
/tmp/ipykernel_8561/1510506507.py:7: DeprecationWarning: Please use atoms.calc = calc
atoms.set_calculator(calculator)
LBFGS: 3 07:55:13 -342.609555 0.374318
LBFGS: 4 07:55:13 -342.655791 0.619672
LBFGS: 5 07:55:13 -342.671952 0.635448
LBFGS: 6 07:55:13 -342.685588 0.410953
LBFGS: 7 07:55:13 -342.695934 0.185126
LBFGS: 8 07:55:13 -342.707676 0.303394
LBFGS: 9 07:55:13 -342.720449 0.393400
LBFGS: 10 07:55:13 -342.732849 0.282384
LBFGS: 11 07:55:13 -342.739609 0.113871
LBFGS: 12 07:55:13 -342.744498 0.145904
LBFGS: 13 07:55:13 -342.748734 0.254533
LBFGS: 14 07:55:13 -342.751916 0.220931
LBFGS: 15 07:55:13 -342.754138 0.085339
LBFGS: 16 07:55:14 -342.755417 0.053165
LBFGS: 17 07:55:14 -342.756478 0.118874
LBFGS: 18 07:55:14 -342.757792 0.145426
LBFGS: 19 07:55:14 -342.758954 0.095402
LBFGS: 20 07:55:14 -342.759746 0.036355
LBFGS: 21 07:55:14 -342.760329 0.043886
LBFGS: 22 07:55:14 -342.760766 0.069188
LBFGS: 23 07:55:14 -342.761198 0.053732
LBFGS: 24 07:55:14 -342.761433 0.028790
LBFGS: 25 07:55:14 -342.761631 0.022934
LBFGS: 26 07:55:14 -342.761783 0.031655
LBFGS: 27 07:55:14 -342.761892 0.023374
LBFGS: 28 07:55:14 -342.761949 0.010753
LBFGS: 29 07:55:15 -342.761996 0.011726
LBFGS: 30 07:55:15 -342.762020 0.017852
LBFGS: 31 07:55:15 -342.762072 0.012363
LBFGS: 32 07:55:15 -342.762076 0.004995
LBFGS: 33 07:55:15 -342.762060 0.004476
LBFGS: 34 07:55:15 -342.762074 0.006189
LBFGS: 35 07:55:15 -342.762062 0.003668
LBFGS: 36 07:55:15 -342.762069 0.001881
LBFGS: 37 07:55:15 -342.762082 0.002786
LBFGS: 38 07:55:15 -342.762088 0.003632
LBFGS: 39 07:55:15 -342.762079 0.002243
LBFGS: 40 07:55:15 -342.762090 0.001409
LBFGS: 41 07:55:16 -342.762078 0.001220
LBFGS: 42 07:55:16 -342.762077 0.001632
LBFGS: 43 07:55:16 -342.762085 0.001006
LBFGS: 44 07:55:16 -342.762064 0.000690
Step Time Energy fmax
LBFGS: 0 07:55:16 -341.994047 3.707106
LBFGS: 1 07:55:16 -342.131367 4.426340
LBFGS: 2 07:55:16 -342.429021 3.370986
LBFGS: 3 07:55:16 -342.724371 0.943083
LBFGS: 4 07:55:16 -342.800045 0.793617
LBFGS: 5 07:55:16 -342.828909 0.437020
LBFGS: 6 07:55:21 -342.872912 0.237315
LBFGS: 7 07:55:21 -342.876481 0.111741
LBFGS: 8 07:55:22 -342.880221 0.075046
LBFGS: 9 07:55:22 -342.883296 0.095753
LBFGS: 10 07:55:22 -342.885190 0.100857
LBFGS: 11 07:55:22 -342.886260 0.065888
LBFGS: 12 07:55:22 -342.887198 0.042198
LBFGS: 13 07:55:22 -342.888047 0.035742
LBFGS: 14 07:55:22 -342.888611 0.043609
LBFGS: 15 07:55:22 -342.888874 0.031778
LBFGS: 16 07:55:22 -342.889163 0.022925
LBFGS: 17 07:55:22 -342.889449 0.022104
LBFGS: 18 07:55:22 -342.889715 0.026693
LBFGS: 19 07:55:22 -342.889827 0.021689
LBFGS: 20 07:55:22 -342.889971 0.017798
LBFGS: 21 07:55:23 -342.890115 0.016524
LBFGS: 22 07:55:23 -342.890204 0.016821
LBFGS: 23 07:55:23 -342.890311 0.011784
LBFGS: 24 07:55:23 -342.890312 0.010002
LBFGS: 25 07:55:23 -342.890379 0.010338
LBFGS: 26 07:55:23 -342.890426 0.012522
LBFGS: 27 07:55:23 -342.890446 0.009515
LBFGS: 28 07:55:23 -342.890470 0.007352
LBFGS: 29 07:55:23 -342.890488 0.007702
LBFGS: 30 07:55:23 -342.890495 0.008391
LBFGS: 31 07:55:23 -342.890538 0.007104
LBFGS: 32 07:55:23 -342.890515 0.006136
LBFGS: 33 07:55:23 -342.890560 0.005694
LBFGS: 34 07:55:24 -342.890564 0.005681
LBFGS: 35 07:55:24 -342.890595 0.005638
LBFGS: 36 07:55:24 -342.890613 0.005325
LBFGS: 37 07:55:24 -342.890593 0.005027
LBFGS: 38 07:55:24 -342.890638 0.004830
LBFGS: 39 07:55:24 -342.890628 0.004811
LBFGS: 40 07:55:24 -342.890655 0.004379
LBFGS: 41 07:55:24 -342.890638 0.004417
LBFGS: 42 07:55:24 -342.890677 0.004661
LBFGS: 43 07:55:24 -342.890642 0.005589
LBFGS: 44 07:55:24 -342.890674 0.004677
LBFGS: 45 07:55:24 -342.890695 0.003663
LBFGS: 46 07:55:24 -342.890697 0.003655
LBFGS: 47 07:55:24 -342.890697 0.005094
LBFGS: 48 07:55:25 -342.890694 0.004730
LBFGS: 49 07:55:25 -342.890706 0.003387
LBFGS: 50 07:55:25 -342.890704 0.003405
LBFGS: 51 07:55:25 -342.890722 0.003285
LBFGS: 52 07:55:25 -342.890737 0.003030
LBFGS: 53 07:55:25 -342.890744 0.002017
LBFGS: 54 07:55:25 -342.890721 0.001996
LBFGS: 55 07:55:25 -342.890721 0.001847
LBFGS: 56 07:55:25 -342.890743 0.001673
LBFGS: 57 07:55:25 -342.890738 0.000996
Step Time Energy fmax
LBFGS: 0 07:55:25 -341.816867 3.637716
LBFGS: 1 07:55:25 -341.917219 4.362194
LBFGS: 2 07:55:26 -342.143875 3.080004
LBFGS: 3 07:55:26 -342.591196 2.065531
LBFGS: 4 07:55:26 -342.676578 1.254161
LBFGS: 5 07:55:26 -342.853274 0.396045
LBFGS: 6 07:55:26 -342.861501 0.281914
LBFGS: 7 07:55:26 -342.875726 0.297499
LBFGS: 8 07:55:26 -342.883727 0.291233
LBFGS: 9 07:55:26 -342.887155 0.137761
LBFGS: 10 07:55:26 -342.890183 0.139746
LBFGS: 11 07:55:26 -342.894227 0.239277
LBFGS: 12 07:55:26 -342.897257 0.203881
LBFGS: 13 07:55:26 -342.898728 0.084899
LBFGS: 14 07:55:26 -342.899574 0.060572
LBFGS: 15 07:55:27 -342.900480 0.124482
LBFGS: 16 07:55:27 -342.901278 0.125108
LBFGS: 17 07:55:27 -342.901797 0.063454
LBFGS: 18 07:55:27 -342.902055 0.022635
LBFGS: 19 07:55:27 -342.902177 0.051499
LBFGS: 20 07:55:27 -342.902406 0.063136
LBFGS: 21 07:55:27 -342.902540 0.040770
LBFGS: 22 07:55:27 -342.902595 0.010407
LBFGS: 23 07:55:27 -342.902655 0.014829
LBFGS: 24 07:55:27 -342.902664 0.021484
LBFGS: 25 07:55:27 -342.902694 0.015904
LBFGS: 26 07:55:27 -342.902721 0.005062
LBFGS: 27 07:55:27 -342.902700 0.008799
LBFGS: 28 07:55:27 -342.902728 0.013920
LBFGS: 29 07:55:28 -342.902768 0.011519
LBFGS: 30 07:55:28 -342.902771 0.004520
LBFGS: 31 07:55:28 -342.902798 0.005936
LBFGS: 32 07:55:28 -342.902790 0.011713
LBFGS: 33 07:55:28 -342.902779 0.011379
LBFGS: 34 07:55:28 -342.902786 0.005375
LBFGS: 35 07:55:28 -342.902779 0.002169
LBFGS: 36 07:55:28 -342.902765 0.004994
LBFGS: 37 07:55:28 -342.902789 0.005137
LBFGS: 38 07:55:28 -342.902790 0.002454
LBFGS: 39 07:55:28 -342.902776 0.000694
Step Time Energy fmax
LBFGS: 0 07:55:28 -341.805736 3.630377
LBFGS: 1 07:55:28 -341.906011 4.360864
LBFGS: 2 07:55:28 -342.132995 3.070995
LBFGS: 3 07:55:29 -342.575795 2.031545
LBFGS: 4 07:55:29 -342.659040 1.238490
LBFGS: 5 07:55:29 -342.830706 0.387560
LBFGS: 6 07:55:29 -342.838962 0.292390
LBFGS: 7 07:55:29 -342.850377 0.238325
LBFGS: 8 07:55:29 -342.861258 0.283475
LBFGS: 9 07:55:29 -342.864656 0.159977
LBFGS: 10 07:55:29 -342.868207 0.123694
LBFGS: 11 07:55:29 -342.872319 0.219206
LBFGS: 12 07:55:29 -342.875880 0.196974
LBFGS: 13 07:55:29 -342.877610 0.088268
LBFGS: 14 07:55:29 -342.878570 0.084757
LBFGS: 15 07:55:30 -342.879534 0.111292
LBFGS: 16 07:55:30 -342.880512 0.115116
LBFGS: 17 07:55:30 -342.881138 0.059364
LBFGS: 18 07:55:30 -342.881456 0.022330
LBFGS: 19 07:55:30 -342.881648 0.046044
LBFGS: 20 07:55:30 -342.881885 0.053383
LBFGS: 21 07:55:30 -342.882048 0.032017
LBFGS: 22 07:55:30 -342.882155 0.014687
LBFGS: 23 07:55:30 -342.882240 0.021256
LBFGS: 24 07:55:30 -342.882318 0.023533
LBFGS: 25 07:55:30 -342.882350 0.013021
LBFGS: 26 07:55:30 -342.882363 0.008467
LBFGS: 27 07:55:30 -342.882427 0.016556
LBFGS: 28 07:55:30 -342.882411 0.018341
LBFGS: 29 07:55:31 -342.882437 0.010138
LBFGS: 30 07:55:31 -342.882491 0.007715
LBFGS: 31 07:55:31 -342.882510 0.011601
LBFGS: 32 07:55:31 -342.882539 0.014853
LBFGS: 33 07:55:31 -342.882517 0.009812
LBFGS: 34 07:55:31 -342.882496 0.004526
LBFGS: 35 07:55:31 -342.882508 0.004830
LBFGS: 36 07:55:31 -342.882536 0.006192
LBFGS: 37 07:55:31 -342.882526 0.003496
LBFGS: 38 07:55:31 -342.882528 0.001782
LBFGS: 39 07:55:36 -342.882551 0.002238
LBFGS: 40 07:55:37 -342.882544 0.003163
LBFGS: 41 07:55:37 -342.882551 0.002128
LBFGS: 42 07:55:37 -342.882529 0.000808
[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 104 ms, sys: 3.72 ms, total: 108 ms
Wall time: 141 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 07:55:39 -341.994038 3.707107
/tmp/ipykernel_8561/1510506507.py:7: DeprecationWarning:
Please use atoms.calc = calc
LBFGS: 1 07:55:39 -342.107696 4.426356
LBFGS: 2 07:55:39 -342.373290 3.317400
LBFGS: 3 07:55:40 -342.579909 2.000243
LBFGS: 4 07:55:40 -342.756639 0.870773
LBFGS: 5 07:55:40 -342.790093 0.374673
LBFGS: 6 07:55:40 -342.809422 0.291128
LBFGS: 7 07:55:40 -342.823264 0.266580
LBFGS: 8 07:55:40 -342.833950 0.124070
LBFGS: 9 07:55:40 -342.835572 0.088253
LBFGS: 10 07:55:40 -342.836857 0.075296
LBFGS: 11 07:55:40 -342.838347 0.068667
LBFGS: 12 07:55:40 -342.839204 0.064898
LBFGS: 13 07:55:41 -342.839661 0.064189
LBFGS: 14 07:55:41 -342.840044 0.051451
LBFGS: 15 07:55:41 -342.840505 0.050056
LBFGS: 16 07:55:41 -342.840820 0.032795
LBFGS: 17 07:55:41 -342.841054 0.030308
LBFGS: 18 07:55:41 -342.841195 0.027926
LBFGS: 19 07:55:41 -342.841300 0.032563
LBFGS: 20 07:55:41 -342.841448 0.026804
LBFGS: 21 07:55:41 -342.841483 0.013311
LBFGS: 22 07:55:41 -342.841538 0.014673
LBFGS: 23 07:55:41 -342.841598 0.019561
LBFGS: 24 07:55:41 -342.841640 0.017394
LBFGS: 25 07:55:41 -342.841670 0.009116
LBFGS: 26 07:55:41 -342.841695 0.009561
LBFGS: 27 07:55:42 -342.841708 0.011721
LBFGS: 28 07:55:42 -342.841772 0.014140
LBFGS: 29 07:55:42 -342.841791 0.009850
LBFGS: 30 07:55:42 -342.841842 0.006582
LBFGS: 31 07:55:42 -342.841835 0.008020
LBFGS: 32 07:55:42 -342.841845 0.010491
LBFGS: 33 07:55:42 -342.841836 0.009795
LBFGS: 34 07:55:42 -342.841877 0.007181
LBFGS: 35 07:55:42 -342.841859 0.006186
LBFGS: 36 07:55:42 -342.841893 0.009587
LBFGS: 37 07:55:42 -342.841898 0.010452
LBFGS: 38 07:55:42 -342.841930 0.007519
LBFGS: 39 07:55:43 -342.841952 0.006661
LBFGS: 40 07:55:43 -342.841949 0.008566
LBFGS: 41 07:55:43 -342.841968 0.010709
LBFGS: 42 07:55:43 -342.842006 0.008396
LBFGS: 43 07:55:43 -342.842020 0.005409
LBFGS: 44 07:55:43 -342.842028 0.004772
LBFGS: 45 07:55:43 -342.842030 0.006176
LBFGS: 46 07:55:43 -342.842019 0.005571
LBFGS: 47 07:55:43 -342.842023 0.002580
LBFGS: 48 07:55:43 -342.842029 0.002196
LBFGS: 49 07:55:43 -342.842009 0.002595
LBFGS: 50 07:55:43 -342.842016 0.002606
LBFGS: 51 07:55:43 -342.842025 0.001621
LBFGS: 52 07:55:44 -342.842046 0.001115
LBFGS: 53 07:55:44 -342.842026 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.8266751227401983, E_slab_mol -342.84202612473905, E_slab -329.5151524229501, E_mol -11.500198579048762
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 07:55:44 -14.872189 4.690001
LBFGS: 1 07:55:44 -14.761179 7.138004
/tmp/ipykernel_8561/1510506507.py:7: DeprecationWarning:
Please use atoms.calc = calc
/tmp/ipykernel_8561/1510506507.py:9: FutureWarning:
Import StrainFilter from ase.filters
LBFGS: 2 07:55:44 -14.923395 0.707045
LBFGS: 3 07:55:44 -14.924615 0.152445
LBFGS: 4 07:55:44 -14.924664 0.001121
LBFGS: 5 07:55:44 -14.924666 0.000014
[30]:
3.9414397766391507
[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 07:55:44 -220.564856 0.115019
LBFGS: 1 07:55:44 -220.571327 0.092847
/tmp/ipykernel_8561/1510506507.py:7: DeprecationWarning:
Please use atoms.calc = calc
LBFGS: 2 07:55:44 -220.582634 0.010980
LBFGS: 3 07:55:44 -220.582743 0.001272
LBFGS: 4 07:55:44 -220.582707 0.000075
[32]:
mol = molecule("CO")
mol.calc = calculator
E_mol = get_opt_energy(mol, fmax=1e-4)
Step Time Energy fmax
LBFGS: 0 07:55:44 -11.497234 0.820158
LBFGS: 1 07:55:45 -11.484527 1.958979
LBFGS: 2 07:55:45 -11.500196 0.029784
/tmp/ipykernel_8561/1510506507.py:7: DeprecationWarning:
Please use atoms.calc = calc
LBFGS: 3 07:55:45 -11.500199 0.001079
LBFGS: 4 07:55:45 -11.500198 0.000010
[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 07:55:45 -233.487857 1.079220
LBFGS: 1 07:55:45 -233.491375 1.894393
/tmp/ipykernel_8561/1510506507.py:7: DeprecationWarning:
Please use atoms.calc = calc
LBFGS: 2 07:55:45 -233.515930 0.671719
LBFGS: 3 07:55:45 -233.535165 0.525847
LBFGS: 4 07:55:45 -233.540258 0.339795
LBFGS: 5 07:55:45 -233.543925 0.176446
LBFGS: 6 07:55:45 -233.548123 0.344875
LBFGS: 7 07:55:45 -233.552475 0.420478
LBFGS: 8 07:55:45 -233.555398 0.288876
LBFGS: 9 07:55:46 -233.557775 0.141699
LBFGS: 10 07:55:46 -233.560570 0.271896
LBFGS: 11 07:55:46 -233.564281 0.313258
LBFGS: 12 07:55:46 -233.568309 0.284822
LBFGS: 13 07:55:46 -233.571348 0.197042
LBFGS: 14 07:55:46 -233.573282 0.113801
LBFGS: 15 07:55:46 -233.574606 0.148960
LBFGS: 16 07:55:46 -233.575601 0.152804
LBFGS: 17 07:55:46 -233.576371 0.098615
LBFGS: 18 07:55:46 -233.576925 0.083781
LBFGS: 19 07:55:46 -233.577258 0.069899
LBFGS: 20 07:55:46 -233.577481 0.047896
LBFGS: 21 07:55:46 -233.577654 0.040826
LBFGS: 22 07:55:46 -233.577835 0.067757
LBFGS: 23 07:55:47 -233.578010 0.069439
LBFGS: 24 07:55:47 -233.578084 0.034979
LBFGS: 25 07:55:47 -233.578156 0.019567
LBFGS: 26 07:55:47 -233.578243 0.037782
LBFGS: 27 07:55:47 -233.578312 0.055524
LBFGS: 28 07:55:47 -233.578400 0.041176
LBFGS: 29 07:55:47 -233.578460 0.014200
LBFGS: 30 07:55:47 -233.578464 0.009982
LBFGS: 31 07:55:47 -233.578483 0.020392
LBFGS: 32 07:55:47 -233.578480 0.019457
LBFGS: 33 07:55:47 -233.578483 0.008402
LBFGS: 34 07:55:47 -233.578493 0.002447
LBFGS: 35 07:55:52 -233.578524 0.008048
LBFGS: 36 07:55:53 -233.578501 0.008466
LBFGS: 37 07:55:53 -233.578520 0.004002
LBFGS: 38 07:55:53 -233.578512 0.000913
CPU times: user 284 ms, sys: 30.7 ms, total: 315 ms
Wall time: 7.95 s
[34]:
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.4956079430033924, E_slab_mol -233.57851244191073, E_slab -220.58270687257885, E_mol -11.500197626328491
[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.