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()