Slab + mol 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
吸着エネルギー¶
吸着エネルギー \(E_{\rm{adsorption}}\) は、ある分子が表面構造に吸着する前後のエネルギー差を表します。
\[E_{\rm{adsorption}} = E_{\rm{slab+mol}} - (E_{\rm{slab}} + E_{\rm{mol}})\]
\(E_{\rm{slab}}\)がSlab単体のエネルギー、\(E_{\rm{mol}}\)が分子単体のエネルギー、\(E_{\rm{slab+mol}}\)が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()
まずは適切なPt のセルサイズを求めます。
Bulkを用意して最適化することでセルサイズを求めます。
[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
ASE デフォルトの3.92A からは少し値が変わっていることがわかります。このセルサイズを用いて以下、吸着構造を作成していきます。
[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
吸着構造の作成¶
吸着構造は add_adsorbate 関数を用いることで、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
吸着エネルギーとして-1.87 eVが得られました。
[12]:
from pfcc_extras.visualize.view import view_ngl
view_ngl(slab_mol)
吸着サイトについて¶
add_adsorbate 関数の4番目に指定する position 引数として、座標値の代わりに吸着サイト名で指定することができます。
[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
この吸着サイトの名前は、 .info を表示することで “sites” 以下に定義されている値として確認ができます。
[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}}
今回の例fcc111だと、 ontop, bridge, fcc, hcp があるようです。それぞれ実際に見てみましょう。
[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()
上記のように、
ontop: ある１つの原子の上bridge: ある２つの原子の間fcc: ３つの原子の間で、2層目には原子が存在しない位置hcp: ３つの原子の間で、2層目の原子の直上の位置
という形となっています。
分子が表面構造のどの部分に吸着するのか(吸着サイト)は、表面・分子系によって異なり自明ではありません。 吸着エネルギーが最小(吸着構造のエネルギーが最小)となる場所を探す必要があります。
[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()
このように、吸着サイトによって吸着構造のエネルギーが異なることがわかります。 今回の計算では、fcc hollowサイトへの吸着構造が一番安定で、吸着エネルギーが一番小さいという結果が得られました。
この結果を以下の文献のTable 1.やTable S2. に記載のPBE 交換相関汎関数を用いて行われた計算結果と比較してみましょう。吸着エネルギーがTop、Bridge、HCP Hollowサイトでそれぞれ1.80, 1.88, 1.89 eV と報告されています。 (ちなみに、実験ではTopサイトが最安定となることが知られており、PBEを用いて行われた計算では実験結果と一致する結果が出ないと報告されています。)
[20]:
view_ngl([slab_mol_ontop, slab_mol_bridge, slab_mol_fcc, slab_mol_hcp])
FixAtoms constraints を使用する例¶
上記の計算例では、すべての原子座標を構造緩和しましたが、DFTを用いて吸着構造を求める場合、構造の安定性を高めるためにSlabの下層の原子座標を固定して構造緩和を行うことがあります。
この方法での計算も行ってみましょう。ASEのconstraintであるFixAtomsを用いることで、指定した 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
今回は下から2層の原子を指定する際に、z座標でしきい値を設定します。
上記の可視化結果をマウスホバーしてもわかりますが、今回は以下のようにz軸の分布を可視化して確認してみます。
[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()
2層目と3層目の境界が、42.29 - 44.58 の間であることがわかりました。
FixAtoms は indicesに固定したい原子のindexを指定するか、maskにそれぞれの原子を固定するかしないかの配列を指定します。 今回はmaskを使用してみます。
[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])
作成したconstraintは、 atoms.set_constraint methodで指定することで制約がかけられます。
[25]:
slab_mol_fixed.set_constraint(constraint)
可視化して、意図した場所を固定できているか確認してみます。
[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
下層を固定する方法を用いた場合、吸着エネルギーとして-1.83 eVが得られました。
参考文献
CO on Pd¶
同様に、Pd上でのCO吸着に関しても計算してみましょう。
[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)
このように、様々な表面や分子の組み合わせに対して吸着エネルギーが計算できます。 これは5章で説明する、効率の良い触媒探索に際して欠かせないプロセスとなります。
被覆率¶
本Tutorialで取り上げた例では、まっさらな表面構造上に、分子が１つだけ吸着するというものでした。 現実では、たくさんの分子が同時に表面と反応を起こし、同時に複数の分子が吸着しているような場合や、ほとんど吸着した分子で覆われているような場合も存在します。
その際には、分子がすでに複数個吸着している表面構造から、あらたに1つの分子が吸着する場合の吸着エネルギーなどにも興味が出てきます。
上記で上げたPd表面上でのCO分子の吸着エネルギーの被覆率依存性の解析に関しては以下の文献で報告されています。
複数箇所への吸着¶
また、今回扱えていない応用トピックとして、2箇所以上に吸着する例があります。 分子が大きく、複雑になってくると複数箇所が表面上に吸着することで安定となるようなものも存在します。