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_U0, model_version="v3.0.0")
calculator = ASECalculator(estimator)
pfp_api_client: 1.18.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
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

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

吸着構造の作成

吸着構造は 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
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

吸着エネルギーとして-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
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

この吸着サイトの名前は、 .info を表示することで “sites” 以下に定義されている値として確認ができます。

[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}}

今回の例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()
_images/3_3_slab_adsorption_energy_26_0.png

上記のように、

  • ontop: ある1つの原子の上

  • bridge: ある2つの原子の間

  • fcc: 3つの原子の間で、2層目には原子が存在しない位置

  • hcp: 3つの原子の間で、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
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()
_images/3_3_slab_adsorption_energy_29_0.png

このように、吸着サイトによって吸着構造のエネルギーが異なることがわかります。 今回の計算では、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 97.7 ms, sys: 7.32 ms, total: 105 ms
Wall time: 112 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.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()

2層目と3層目の境界が、42.29 - 44.58 の間であることがわかりました。

FixAtomsindicesに固定したい原子の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
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

下層を固定する方法を用いた場合、吸着エネルギーとして-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
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)

このように、様々な表面や分子の組み合わせに対して吸着エネルギーが計算できます。 これは5章で説明する、効率の良い触媒探索に際して欠かせないプロセスとなります。

被覆率

本Tutorialで取り上げた例では、まっさらな表面構造上に、分子が1つだけ吸着するというものでした。 現実では、たくさんの分子が同時に表面と反応を起こし、同時に複数の分子が吸着しているような場合や、ほとんど吸着した分子で覆われているような場合も存在します。

その際には、分子がすでに複数個吸着している表面構造から、あらたに1つの分子が吸着する場合の吸着エネルギーなどにも興味が出てきます。

上記で上げたPd表面上でのCO分子の吸着エネルギーの被覆率依存性の解析に関しては以下の文献で報告されています。

複数箇所への吸着

また、今回扱えていない応用トピックとして、2箇所以上に吸着する例があります。 分子が大きく、複雑になってくると複数箇所が表面上に吸着することで安定となるようなものも存在します。