Developments of Interorbital–band Interaction Analysis and Embedded Cluster Model Incorporating Periodic Electrostatic Potential for Supported Metal Catalysts

Abstract

担体表面上に金属微粒子が高分散した担持金属触媒において，金属微粒子と担体表面との間の金属–表面相互作用は，金属微粒子を担体表面上に安定に分散担持するとともに，金属と表面間の電荷移動などの相互作用により電子状態を変化させ，触媒活性に重要な影響を及ぼす．このような金属微粒子と担体表面との界面における電子状態の関与する相互作用を，実験の測定のみから解明することは困難であり，その本質の理解・予測のためには理論計算による検討が有用である．しかし，従来のスラブモデルを用いた平面波DFT法には，軌道間相互作用に基づく解析手法が乏しい，また，高精度電子状態計算の実行がコスト的に困難，という問題がある．このような問題を解決するために我々は，「射影状態密度を用いた軌道–バンド間相互作用解析手法」と「周期的静電ポテンシャルへの埋め込みクラスターモデル」の開発を行ってきた．本総説では，これらの手法の概要を述べるとともに，Rh₂/AlPO₄，Rh₂/Al₂O₃への適用例を紹介する．

1 はじめに

担持金属触媒は金属微粒子を担体表面に分散・担持した触媒であり，工業的にひろく利用されている．特に自動車排気ガス浄化三元触媒として，金属にPtやRh，Pdなどの白金族を，担体にAl₂O₃やCeO₂などの金属酸化物を用いたものが使用されているが，希少な白金族金属の産出量の大半がその用途に消費されていることから，貴金属の使用量低減及び代替が喫緊の課題となっている．そのためのアプローチの一つとして，金属微粒子と担体表面との間の金属–表面相互作用を利用する方法が提案されている [1]．金属–表面相互作用は金属微粒子を担体表面上に安定に担持することによって活性低下の原因となる粒子成長を抑制するとともに，金属と表面間の電荷移動などによりHOMO，LUMOのような反応に直接関与する電子状態を変化させ，触媒活性に影響を与える [2]．実際に，担体を従来のAl₂O₃からAlPO₄に変えることによりRhの粒子成長が抑制され，三元触媒としての活性が向上することが報告されている [3,4,5,6,7]．このような金属微粒子と担体表面との界面において電子状態が関与する相互作用を，実験的手法のみから解明することには限界があり，その本質の理解・予測のためには理論計算と組み合わせた検討が必要不可欠である．例えば，Pt/γ-Al₂O₃においては実験的な測定と密度汎関数理論 (DFT) 計算を組み合わせることによって，表面に配位不飽和Alが存在することがPtの安定担持に寄与していることが明らかにされた [8]．またPt/CeO₂系，Pt/CeO₂/TiO₂系においても同様に，担体との相互作用によるPtの5d殻の電子密度の減少がH₂O吸着能，O–H結合解離能の増大に寄与することが示されている [9]．

このような表面系の理論計算では，スラブモデルを用いた平面波DFT法によるバンド計算が一般的に採用されている．真空層をsuper cell中に設けることにより表面を無限周期的に記述するスラブモデルと，系全体に非局在化した電子状態を記述できる平面波DFT法を採用することによって，巨視的な表面全体の影響を構造や電子状態の理論計算に取り入れることが可能である．しかし平面波DFT法には，電子状態の解析手法が乏しく，特に軌道間相互作用に基づく解析手法がない，さらに，LDAやGGAなどの交換相関汎関数に比べ，電子間のパウリ反発や多体問題をより高精度に記述することが可能なhybrid汎関数の使用やpost Hartree-Fock計算の実行がコスト的に困難である，という大きな問題がある．第一原理計算によって得られた電子状態の解析に関して，分子科学分野ではMullikenの電荷移動錯体の研究 [10] に代表されるように，電荷移動などの相互作用を軌道間相互作用に基づいて解析する手法が多数提案されている．だが結晶系ではその1電子軌道の集合であるバンドが系全体に非局在化し，波動関数がその並進対称性を利用して逆空間で表現されることから，分子科学分野での解析手法はこれまでほとんど利用されてこなかった．しかし，担持金属触媒の金属–表面相互作用，さらに触媒化学反応では，表面と，金属，反応物質間のフロンティア軌道間の相互作用が重要となることから，軌道間相互作用に基づく解析手法は必要不可欠である．また，高精度な電子相関を取り込んだ電子状態計算に関しては，周期的境界条件下での固体の計算に取り入れる試みはなされているが [11,12,13,14,15,16,17,18,19,20,21,22]，やはり計算コストの問題で，小さなunit cellを対象にしたものに限られており，隣接cell間でのartificialな相互作用を避けるためにある程度以上大きいunit cellを必要とする担持金属触媒のスラブモデルへの適用は困難である．しかし，担持金属触媒の触媒反応ではNO還元などの開殻系分子の寄与する複雑な遷移状態を考慮しなければならず，高精度な電子状態計算手法の適用もまた必要不可欠である．

我々はこれまで，結晶系における軌道間相互作用の解析手法の開発を行なってきたが [23]，この手法を表面系での射影状態密度に応用することにより，担持金属触媒における金属と表面間の軌道–バンド間相互作用の解析を可能とした [24] ．また，巨視的な表面系において高精度電子状態計算を行う手法として，周期的静電ポテンシャルへの埋め込みクラスターモデルを提案し，開発を行なってきた [25]．本総説では，これらの手法について解説するとともに，自動車三元触媒Rh/AlPO₄系およびRh/Al₂O₃に適用した研究例を紹介する．

2 軌道–バンド間相互作用解析

2.1 理論的背景

平面波基底を用いた電子状態計算の解析手法としては，電荷密度の空間分割に基づくBader電荷解析 [26,27,28,29,30] がよく用いられている．これは，電荷密度を原子間 (より正確には電荷密度極大点間) の最小点を横切る垂直面によって多面体に分割し，その多面体の空間内の電荷の総和を求めることによって原子の形式電荷を見積もる手法である．この手法は，Mulliken電子分布解析 [31] などと異なり基底の種類によらずに用いることが可能であるが，一方で軌道の描像を用いないので，原子間，または分子間での波動関数の重なりに基づいて相互作用を評価することはできない．結晶系における軌道間相互作用に基づく解析手法としては，COOP (Crystal-Orbital Overlap Population) [32]，COHP (Crystal-Orbital Hamiltonian Population) [33] 解析があるが，これらは局在化基底に対して開発されており，無限周期的に非局在化した平面波基底では使用できない．筆者らの一人は，平面波DFT法において軌道間相互作用に基づく解析を行うために，Bloch軌道線形結合 (LCBO; Linear Combination of Bloch Orbitals) を用いる手法の導出と開発を行った [23]．ここではこの手法の概要と，担持金属触媒の金属–表面相互作用の解析への応用に関して解説する．

(A) LCBO展開の利用

LCBOを導出するために，super cellアプローチを導入する．ここで，super cellの格子ベクトルは十分に大きく，隣接super cell間でのBloch軌道波動関数の重なりは無視できるとして，周期的境界条件下での孤立系を定義する．このとき，結晶系のBloch軌道はその分割系である孤立系のBloch軌道の線形結合として以下のように表される:   

${\psi }_{i,k}\left(r\right)=\sum _{R^{\text{I}}}\sum _{F,n}{\varphi }_{Fn,k^I}\left(r-R^{\text{I}}\right)C_{Fn i,k}{e}^{i\left(k^{\text{I}}+K\right)\cdot R^{\text{I}}}$

(2.1)

ここでψ_i,k，ϕ_Fn,kIは各々結晶系，孤立系Fのi，n番目のBloch軌道波動関数，R^Iはsuper cell中の孤立系の格子点，C_{Fn i,k}はLCBO展開係数，k，k^Iは各々結晶系，孤立系の波数ベクトルで，k = k^I + Kである．平面波基底においてこのように全系のBloch軌道を分割系に局在化した関数で表す手法としては，他に局在化原子基底関数への射影を行う手法 [34,35,36,37,38] や，局在化Wannier関数へと変換する手法 [39] が存在するが，今回のLCBOを用いた手法ではそれらと異なり，軌道の占有・非占有状態を明確に定義可能である．以上のように全系の結晶系Bloch軌道をその分割系の孤立系Bloch軌道で展開できたので，軌道間の重なり積分やHamiltonian積分が評価可能となり，分子系において開発された電荷分解解析 [40] のような軌道間相互作用に基づく手法が平面波DFT法によるバンド計算にも適用可能となった．より詳細は参考文献23を参照されたい．

(B) p-DOSによる解析

担持金属触媒の金属と表面間の相互作用を軌道–バンド間相互作用に基づいて解析するために，この開発した手法を射影状態密度 (p-DOS; projected density of states) に応用した [24]．全系バンドのDOSにおける分割系Fの軌道nへのp-DOSは，全系軌道と分割系軌道の重なり積分を用いて以下のように表される:   

${\rho }_{Fn}\left(\epsilon \right)=\sum _k{\omega }_k\sum _i{\left|⟨{\varphi }_{Fn,k^I}\text{|}{\psi }_{i,k}⟩\right|}^2\delta \left(\epsilon -{\epsilon }_{i,k}\right)$

(2.2)

ω_kはk-point samplingの重み，ε_i,kは全系軌道の軌道エネルギーである．ここで，全系として金属クラスターが担体表面に吸着した系を，その分割系として，金属クラスター部分と担体表面部分とに分けて考える．このとき全系軌道は分割系軌道の線形結合で表されるので，全系バンドのDOSにおける分割系軌道へのp-DOSは，全系バンド構造への分割系軌道の寄与を表す．よって，例えば担体表面の導電バンド (非占有軌道) へのp-DOSが吸着系の価電子バンド (占有軌道) のDOSに現れたら，それは担体表面導電バンドが金属クラスターの吸着により電子を受け取ったことを示している．このように，軌道–バンド間相互作用に基づくp-DOSによる解析は，担持金属触媒での金属クラスターと担体表面との電荷移動を伴う金属–表面相互作用を調べるうえで有用な手法である．

2.2 適用例の紹介

このp-DOSによる手法を用いて，Rh/AlPO₄系とRh/Al₂O₃系の金属–表面相互作用を解析し，両者の比較により表面の特性評価を行った [24]．Rhクラスターとしては最小のRh₂を採用した．AlPO₄は110面をtridymite構造 [41] より切り出し，Al₂O₃は100面を脱水boehmite構造 [42] より切り出した．Rh₂/Al₂O₃吸着構造はShi and Sholl [43] を参考にした．計算にはVASP [44,45,46] を用いた．その他の計算条件の詳細については参考文献24を参照されたい．

Figure 1にRh₂/AlPO₄とRh₂/Al₂O₃の吸着構造を示した．AlPO₄では，Rh₂は2個のPOサイトに架橋して吸着し，その近傍に3配位Alが存在した(Figure 1A)．Al₂O₃では，Rh₂は1原子が表面Oに吸着し，もう1原子のRhは表面からやや離れて位置しており，近傍に5配位Alが存在する (Figure 1B)．Table 1にRh₂の担体表面への吸着エネルギー (E_ad)，Rh₂と担体表面の吸着に伴う構造変化エネルギー (E_df[Rh₂], E_df[Support])，Rh₂と担体表面間の相互作用エネルギー (E_int)をまとめた．ここで，E_adは平衡構造のRh₂と担体表面と，吸着系とのエネルギー差であり，E_df[Rh₂]とE_df[Support]は各々Rh₂と担体表面の平衡構造から吸着構造への変化に伴う不安定化エネルギー，E_intは吸着により構造変化したRh₂と担体表面と，吸着系とのエネルギー差であって，以下のように表される:    E_ad = E_eq[Rh₂/Support] − E_eq[Rh₂] − E_eq[Support]       E_df[Rh₂] = E_dis[Rh₂] − E_eq[Rh₂]       E_df[Supprt] = E_dis[Support] − E_eq[Support]       E_int = E_eq[Rh₂/Support] − E_dis[Rh2] − E_dis[Support]    ここで添え字のeqは平衡構造を，disは吸着構造変化により歪んだ構造を示す．このとき，E_adはE_intとE_dfの寄与に分解されて表される:    E_ad = E_int + E_df[Rh₂] + E_df[Supprt]   

Figure 1.

 Structures of (A) Rh₂/AlPO₄ and (B) Rh₂/Al₂O₃ systems: Green, yellow, red, and black balls indicate Al, P, O, and Rh atoms, respectively. The second and lower layers are represented by pale colors. (Adapted with permission from M. Matsui, M. Machida, and S. Sakaki, J. Phys. Chem. C 2015, 119, 19752–19762. Copyright (2018) American Chemical Society.)

Table 1  Adsorption, deformation, and interaction energies together with important geometrical parameters of Rh₂/AlPO₄ and Rh₂/Al₂O₃ systems. (Adapted with permission from M. Matsui, M. Machida, and S. Sakaki, J. Phys. Chem. C 2015, 119, 19752–19762. Copyright (2018) American Chemical Society)

Table 1において，吸着安定化を示すE_adはRh₂/AlPO₄ (−3.63 eV) の方がRh₂/Al₂O₃ (−2.41 eV) よりも大きく，これはアンカー効果がRh/AlPO₄の方がRh/Al₂O₃よりも大きいという実験結果 [3,4,5,6,7] と一致している．両者でE_df[Rh₂]の差はほとんどないが，E_df[Support]の差は大きい．一方でE_intの差も小さいので，この吸着安定化の違いには，吸着に伴う担体表面の構造変化の寄与が大きいことがわかる．

Figure 2に，平衡構造の担体表面，及び吸着により構造変化した担体表面の最低非占有バンド (以後，LUMOバンドと表記) を，吸着差電子密度とともに示す．Rh₂/AlPO₄では，吸着前のLUMOバンドは表面全体に非局在化しているが (Figure 2A-1)，構造変化により吸着位置近傍の3配位Al上に局在化する (Figure 2A-2)．差電子密度では，Rh₂と吸着位置近傍の3配位Alとの間での電子密度の増大が明らかであり (Figure 2A-3)，このLUMOバンドへのRh₂からの電荷移動が示唆されている．Rh₂/Al₂O₃では，同じように平衡構造でのLUMOバンドは表面全体に非局在化し，構造変化により吸着位置近傍の5配位Al上に局在化する (Figure 2B-1, 2)．差電子密度では，増大領域と減少領域は複雑に存在しており，Rh₂と担体表面間の電荷移動は明確ではない (Figure 2B-3)．

Figure 2.

 Orbital plots of the lowest unoccupied band (LUMO band) of equilibrium and deformed surfaces, and difference electron density map for AlPO₄ and Al₂O₃ systems. The red/blue colors indicate the different sign of the wave function. The cyan area represents an increase in electron density and the orange one does its decrease. (Adapted with permission from M. Matsui, M. Machida, and S. Sakaki, J. Phys. Chem. C 2015, 119, 19752–19762. Copyright (2018) American Chemical Society.)

この表面構造変化により生じた局在化LUMOバンドが，吸着系の電子状態にどのように影響しているかを明らかにするために，p-DOSを用いた解析を行った．平衡構造の担体表面，及び吸着により構造変化した担体表面のDOSと，吸着系のDOS，及びLUMOバンドのp-DOSをFigure 3に示した．ここでバンドエネルギーの基準は，スラブモデルの真空層中心の静電ポテンシャルの値を0として決めた．AlPO₄のDOSにおいて，平衡構造ではバンドギャップが大きく開いていたものが，構造変化によりFermi準位のすぐ上に孤立したバンドが生成しており (Figure 3A-1)，これが3配位Al上局在化LUMOバンドである．吸着系Rh₂/AlPO₄のDOSにおいて，LUMOバンドのp-DOSは，ほぼ全てが価電子バンドに現れている (Figure 3A-2)．これはLUMOバンドがRh₂吸着後は電子に占有されたことを意味しており，すなわちRh₂から3配位Al上局在化LUMOバンドへの電荷移動が起きたことを意味している．Rh₂/Al₂O₃においても，吸着構造変化によりFermi準位より上に孤立LUMOバンドが生成するが，そのバンドエネルギーはAlPO₄よりも高い (Figure 3B-1)．Rh₂/Al₂O₃のDOSにおいてLUMOバンドのp-DOSは，価電子バンドと導電バンドの両方に現れている (Figure 3B-2)．これは，Rh₂/Al₂O₃においてもRh₂吸着により5配位Al上局在化LUMOバンドへの電荷移動が起こるが，その程度はRh₂/AlPO₄よりも小さいことを意味している． Mullikenの電荷移動錯体の研究 [10] で示されたように，非占有バンドのエネルギーが高いほどRh₂からの電荷移動は起こりにくいので，LUMOバンドエネルギーが高いことがAl₂O₃の方がAlPO₄よりも電荷移動が起こりにくいことの要因である．以上のように，Rh₂吸着に伴う構造変化により軌道準位の低下した局在化LUMOバンドが生成し，このLUMOバンドへのRh₂からの電荷移動が起きていること，また，Rh₂/AlPO₄の方がRh₂/Al₂O₃よりも電荷移動が大きいことが，p-DOSを用いた解析より明らかになった．実験的にも，Rh/AlPO₄の方がRh/Al₂O₃よりRhから担体表面への電荷移動が大きいということが示唆されており [6]，今回の結果はこの実験結果をよく説明するものとなっている．

Figure 3.

 (A-1) Density of states (DOS) of equilibrium and deformed surfaces of AlPO₄, and (A-2) DOS of Rh₂/AlPO₄ and p-DOS onto the LUMO band of deformed AlPO₄ surface. (B-1) DOS of equilibrium and deformed ones of Al₂O₃, and (B-2) DOS of Rh₂/Al₂O₃ and p-DOS onto the LUMO band of deformed Al₂O₃ one. Vertical pink line represents the Fermi level. The vacuum level is set to be zero. The negative density of states corresponds to the density of β spin states. (Adapted with permission from M. Matsui, M. Machida, and S. Sakaki, J. Phys. Chem. C 2015, 119, 19752–19762. Copyright (2018) American Chemical Society.)

これらのことから，Rh₂/AlPO₄，Rh₂/Al₂O₃における金属–表面相互作用は以下のような機構によって生成していると理解できる．Rh₂吸着により担体表面の構造変化が起き，軌道準位が低下した局在化LUMOバンドが生成する．そして，そのLUMOバンドへのRh₂からの電荷移動が起きるが，AlPO₄の方がLUMOバンドのエネルギーが低いため，Al₂O₃よりも電荷移動が大きい．それとともに，AlPO₄の方が吸着構造変化に伴う不安定化エネルギーが小さいため，全体の吸着安定化が大きい．これらの表面構造変化に対する挙動の違いは，AlPO₄とAl₂O₃の結合の性質から以下のように理解できる: AlPO₄は共有結合性が強くAl₂O₃はイオン結合性が強い．よって表面構造変化エネルギーについては，共有結合性が強いAlPO₄の方が柔軟に変化可能であり不安定化が小さい．また軌道エネルギーについては，イオン結合性が強いAl₂O₃では元々軌道エネルギーが高く，一方でAlPO₄は表面配位不飽和性が高いので軌道エネルギーが低い．

このように担体の結合の性質や表面特性から，系統的に金属–表面相互作用の予測ができれば，金属–表面相互作用を利用した担持金属触媒の設計に有用な指針を与えることが可能であり，現在，p-DOSによる解析手法を用いて，検討を進めている．一方で，今回のp-DOSによる解析は電荷移動のみを対象としており，より高度な解析を行うためには，他の静電，分極，交換反発といった相互作用にも解析対象を拡大する必要があり，また，より定量的な評価も必要である．分子系においてはこのような相互作用の定量的な解析手法はすでに存在する [40]．よって，今回開発した手法を，金属クラスター軌道と無限周期的な表面バンドとの相互作用の定量的解析手法へと発展させることが今後は望まれる．

3 埋め込みクラスターモデル

3.1 理論的背景

担持金属触媒のような大規模系における高精度計算の実行という問題を克服するアプローチとして，埋め込みクラスターモデルがある [47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65]．これは大規模系の中の着目する一部分のみの構造をクラスターとして切り出して，残りの系全体の効果を外部ポテンシャルなどで表現し，それに埋め込んだクラスターの高精度分子軌道計算を行う手法である．この埋め込みクラスターモデルを担持金属触媒に適用するにあたっては，クラスターと担体表面との静電相互作用の評価が重要となる．何故なら，担体として用いられる金属酸化物はイオン結晶性が強いために分極が大きく，表面全体との長距離静電相互作用を打ち切ることなく考慮する必要があるからである．このような長距離静電相互作用を考慮した埋め込みクラスターモデルとしては，密度汎関数埋め込み理論 (DFET; Density Functional Embedding Theory) を周期系に応用した手法 [54] と，周期的高速多重極子法 [66,67] を用いた周期的静電埋め込みクラスター法 (PEECM; Periodic Electrostatic Embedding Cluster Method) [59] が存在するが，最近我々は新たに，super cell以外の近似を用いずにより厳密に周期的静電ポテンシャルの1電子積分を評価する，周期的静電ポテンシャル埋め込みクラスターモデルの開発を行なった [25]．ここでは，これらの周期的静電ポテンシャルへの波動関数の埋め込み手法について解説する．

3.2 周期的静電埋め込みクラスターモデルの概要

周期的静電ポテンシャル埋め込みクラスターモデルの計算手順を，模式的にFigure 4に示した．まずは，Figure 4A に示した無限周期的なスラブモデルを用いて，バンド計算や古典力場による計算を行い，表面の構造と静電ポテンシャルを求める．次に，Figure 4Bに示すように，担持金属クラスターとその近傍の表面など，着目する領域の構造をクラスターモデルとして切り出す．そしてFigure 4Cで示すように，スラブモデルによって求めた静電ポテンシャルのもとでクラスターモデルによる分子軌道計算を行う．これにより表面全体の影響を考慮した高精度電子状態計算が可能となる．

Figure 4.

 Schematic representation of (A) slab model, (B) bare cluster model, and (C) embedded cluster model incorporating electrostatic potential. (Adapted with permission from M. Matsui and S. Sakaki, J. Phys. Chem. C 2017, 121, 20242–20253. Copyright (2018) American Chemical Society.)

この埋め込みクラスターモデルにおいて，周期的静電ポテンシャルをクラスターの電子軌道に作用させるためには，一般的な分子軌道計算では1電子軌道波動関数をガウス基底により展開するので，ガウス基底における1電子積分の表式が必要となる．そのためにはまず，クラスターの軌道を周期的ポテンシャルの固有状態であるBloch軌道として表すために，周期的ガウス基底関数を定義する;   

${\varphi }_{Aa,k}^{\text{PBC}}\left(r\right)=\sum _R{\varphi }_{Aa}\left(r-R\right)e^{ik\cdot R}$

(3.1)

ここでϕ_Aaは位置r_Aにある原子のa番目の短縮ガウス基底関数であり，kは逆空間の波数ベクトル，Rは実空間の格子ベクトルである．ここでsuper cellアプローチを導入し，super cellの格子ベクトルは十分に大きく，隣接cell間での基底関数間の重なりは無視できるとする;   

${\varphi }_{An}\left(r-R\right){\varphi }_{Bm}\left(r-R^{\prime }\right)\approx 0 \left(R\ne R^{\prime }\right)$

(3.2)

このとき，クラスターの原子は全て同じunit cellに属するとする．この結果，バンド分散は消失するので，逆空間の積分はΓ点; k = (0, 0, 0) のみを考慮すれば良い．周期的静電ポテンシャルV^ES(r)の周期的ガウス基底関数での1電子積分の表式は，Poisson方程式とフーリエ変換，及び式 3.2を使って，以下のように表される:   

$⟨{\varphi }_{Aa}^{\text{PBC}}\text{|}{\hat{V}}^{\text{ES}}\text{|}{\varphi }_{Bb}^{\text{PBC}}⟩=\sum _G{V}^{\text{ES}}\left(G\right){\int }_{all}dr {\varphi }_{Aa}^{}\left(r\right){\varphi }_{Bb}^{}\left(r\right)e^{iG\cdot r}$

(3.3)

ここでGは逆空間の格子ベクトル，V^ES(G)はV^ES(r)のフーリエ変換である．右式の積分は，ガウス基底関数の積のフーリエ変換であり，これが得られれば周期的静電ポテンシャルの1電子積分が評価できる．ガウス基底関数が (単純) ガウス関数の場合， ガウス関数の積はガウス関数であり，ガウス関数のフーリエ変換もまたガウス関数であるので，これは解析的に評価可能である．一方で，ガウス基底関数が，角運動量を持つ原子軌道を表すカーテシアンガウス関数の場合は，その表式と評価方法は，先に紹介した各々の埋め込みクラスターモデルによって異なる．また，点電荷の作る周期的静電ポテンシャルに関して，逆空間における高波数成分が原因で，式3.3のフーリエ展開が収束しないことが知られているが，この問題への対応も各々の埋め込みクラスターモデルにより異なる．以下では，これらの各モデルでの違いについて簡潔に説明する．

周期系でのDFETによる埋め込みクラスターモデルにおいては，式3.3に現れるカーテシアンガウス関数の積を含む積分に関して，1中心の場合はフーリエ変換を再帰的表式を用いて評価し (参考文献54の式46) ，2中心の場合は実空間のgridでの数値積分によって評価している．点電荷静電ポテンシャルに関しては，DFETにおいては外部ポテンシャルに擬ポテンシャルと電荷密度によるポテンシャルも考慮することにより高波数成分が消えるので，低波数成分のみの評価を行っている．

PEECMによる埋め込みクラスターモデルでは，カーテシアンガウス関数の積は，球面調和関数で展開されて表される (参考文献59の式5，7，10)．点電荷静電ポテンシャルの1電子積分は，周期的高速多重極子法 [66,67] を用いて評価される．

我々の提案した埋め込みクラスターモデルにおいては，二項展開を用いることにより，2中心のカーテシアンガウス関数の積のフーリエ変換も再帰的表式により評価可能とした．導出の詳細については，参考文献25のSupporting Information (Appendix A) を参照されたい．点電荷静電ポテンシャルに関しては，Ewaldの方法を，周期的点電荷分布と1電子軌道との静電相互作用の評価に応用した．このとき，誤差関数の実空間での1電子積分を求める必要があるが，誤差関数はガウス電荷分布の静電ポテンシャルであることを利用して，2電子積分の表式を用いて評価した．導出の詳細については，上述の文献25 (Appendix B) を参照されたい．このようにして，super cell以外の近似を用いない周期的静電ポテンシャルの1電子積分の厳密な表式を得ることにより，無限周期的な表面とクラスターの1電子軌道との長距離静電相互作用の高速・高精度な評価が可能になった．

3.3 適用例の紹介

クラスターモデルの構造に関しては，Rh₂/Al₂O₃, Rh₂/AlPO₄のスラブモデル (参考文献24参照) から， Rh₂/(Al₂O₃)₁₂, Rh₂/(AlPO₄)₁₅を切り出して構築した．これらの構造をFigure 5に示す．静電ポテンシャルに関しては，スラブモデル計算によって得られたBader電荷を用いた周期的点電荷分布から求めた．長距離静電相互作用の影響を検討するうえで，今回の無限周期的な静電ポテンシャルを考慮した埋め込みクラスターモデルと，有限の点電荷を用いて静電ポテンシャルを表現した埋め込みクラスターモデルとを比較することは有用である．そこで非常に多数の点電荷をクラスターを取り囲んで遠距離にまで配置することにより表面の静電ポテンシャルの影響を取り入れる，有限点電荷埋め込みモデルも構築した．以後，有限点電荷埋め込みクラスターモデルをPC (Point Charge) モデル，周期的静電ポテンシャル埋め込みクラスターモデルをPE (Periodic Electrostatic) モデルと表記する．PEモデルの計算はgamess [68] に1電子積分計算法を実装して行い，PCモデルの計算はGaussian09 [69] を用いて行なった．その他の計算条件の詳細については，参考文献25を参照されたい．

Figure 5.

 Structures of (A) slab model, (B) bare cluster model, and (C) embedded cluster model for Rh₂/Al₂O₃, and (D) slab model, (E) bare cluster model, and (F) embedded cluster model for Rh₂/AlPO₄: The green, yellow, red, and black balls indicate Al, P, O, and Rh atoms, respectively. The layers below the first and second layers are represented by sticks in slab models. The atoms in external region of the cluster are represented by pale colors in embedded cluster models. (Reprinted with permission from M. Matsui and S. Sakaki, J. Phys. Chem. C 2017, 121, 20242–20253. Copyright (2018) American Chemical Society.)

前章で述べたように，Rh₂/Al₂O₃, Rh₂/AlPO₄の金属–表面相互作用では，Rh₂吸着に伴う表面構造変化により吸着サイトに軌道準位の低下した局在化LUMOバンドが生成すること，そのLUMOバンドへのRh₂からの電荷移動が起こることが重要であった．今回開発した埋め込みクラスターモデルが，このような金属–表面相互作用を十分に記述可能なものであるためには，(1) Rh₂–担体表面間相互作用エネルギー，(2) 担体表面のHOMO，LUMO分布とエネルギーギャップ，そして (3) Fermi準位付近の状態密度 (DOS) が，スラブモデルの結果を再現する必要がある．

最初に，Rh₂–担体表面間相互作用エネルギーの検討結果を紹介する．Table 2にRh₂/Al₂O₃, Rh₂/AlPO₄のスラブモデル，静電ポテンシャルを考慮しない孤立クラスターモデル (以後，bareモデルと表記)，PCモデル，及びPEモデルを用いて評価したRh₂–担体表面間相互作用エネルギー (E_int) をまとめた．PBE汎関数を使用したRh₂/Al₂O₃のスラブモデルの計算では，E_intは-5.44 eVになるが，bareモデルでは-7.18 eVと大きく過大評価した．それに対して約1000個から10⁶個までの点電荷を用いたPCモデルでは，点電荷数が増大するにしたがいE_intは一定値に近づき，最終的にBSSE補正前は-5.99 eV，BSSE補正後は-4.93 eVとなり，スラブモデルの結果に近い値となった．このとき，約1000個のPCモデルではクラスターから約25 Åまでの範囲の点電荷を，約10⁶個のPCモデルでは同様に約1000 Åまでの範囲の点電荷を考慮したことに対応する．PEモデルのE_intは，BSSE補正前もBSSE補正後も，最も多くの点電荷数を用いたPCモデルの結果とほぼ一致した．Rh₂/AlPO₄では点電荷数が少ないとPCモデルの電子状態のSCF計算が収束しなかったが，最も多くの点電荷数を用いたPCモデルでは収束し，そのE_intはスラブモデルの値と近く，またPEモデルの結果ともほぼ一致した．これらの結果から，表面の作り出す長距離静電ポテンシャルを考慮することにより，埋め込みクラスターモデルを用いて妥当な値のRh₂–担体表面間相互作用エネルギーが得られることがわかった．さらに，Table 2に示したように，スラブモデルでは現実的に実行不可能であるhybrid汎関数のB3LYPを用いた計算も可能となった．B3LYPではPBEよりもRh₂/Al₂O₃では0.5 eV程度，Rh₂/AlPO₄では0.2 eV程度E_intを小さく見積もる傾向があったが，これは後に示すようにHOMO–LUMOギャップがPBEよりも大きくなることが関与しているのではないかと考えられる．

Table 2  Interaction energies (E_int; in eV) of Rh₂/Al₂O₃ and Rh₂/AlPO₄ slab and cluster models. (Reprinted with permission from M. Matsui and S. Sakaki, J. Phys. Chem. C 2017, 121, 20242–20253. Copyright (2018) American Chemical Society.)

^aNumbers of point charges are 1060 (24 × 28 × 15 ų), 11940 (83 × 83 × 15 ų), 1.079 × 10⁵(250 × 250 × 15 ų), 5.879×10⁵(590 × 590 × 15 ų), and 1.452 × 10⁶(920 × 920 × 15 ų) for the Al₂O₃ embedded models with very small, small, middle, large, and very large numbers of point charges (VSP, SP, MP, LP, and VLP), respectively. Those are 2934 (53 × 50 × 15 ų), 8310 (84 × 84 × 15 ų), 7.551 × 10⁵(270 × 250 × 15 ų), 4.115 × 10⁵(620 × 580 × 15 ų), and 1.016 × 10⁶(970 × 920 × 15 ų) for the AlPO₄'s with VSP, SP, MP, LP, and VLP, respectively. ^bPE indicates periodic electrostatic potential.^cRh₂/Al₂O₃ and Rh₂/Al₂O₃-L mean Rh₂/(Al₂O₃)₁₂ and Rh₂/(Al₂O₃)₁₈. Rh₂/AlPO₄ and Rh₂/AlPO₄-L mean Rh₂/(AlPO₄)₁₅ and Rh₂/(AlPO₄)₁₉.^dIn parentheses, the E_int with BSSE correction are presented.^eB3LYP calculation cannot be performed with the slab model in VASP.^fSCF calculation cannot be converged in the case of PBE functional.

次に，担体表面のHOMO，LUMO分布とエネルギーギャップを検討した結果を紹介する．Figure 6にRh₂吸着によって構造変化が生じた担体表面のHOMO–LUMOの分布を，Table 3にその軌道エネルギーとエネルギーギャップを示した．Al₂O₃表面のHOMO，LUMOは，スラブモデルでは吸着位置近傍に局在化するが，bareモデルではdangling bondの存在するクラスターの末端領域に非局在化してしまう．このとき，bareモデルのHOMO–LUMOギャップはほぼ閉じている．対照的に，点電荷数の多いPCモデルとPEモデルでは，スラブモデルとほぼ同じHOMO，LUMO分布が得られている．HOMO–LUMOギャップは最も多くの点電荷数を用いたPCモデル，PEモデルで各々1.81，1.82 eVとなり，スラブモデルの1.60 eVと近い値となった．これらの結果からRh₂/Al₂O₃においては，長距離静電相互作用を考慮した静電ポテンシャルへの埋め込みによって，クラスターの末端領域に存在するdangling bondの影響がFermi準位付近の電子状態から取り除かれ，クラスターモデルを用いても妥当なHOMO，LUMOの描像を得ることができることが示された．同様の検討をRh₂/AlPO₄においても行ったところ，HOMOが末端のPO近傍に局在化した軌道となったが，この場合もHOMO–2はスラブモデルと同じ描像の軌道となり， LUMOとHOMO–2のエネルギーギャップもスラブモデルのバンドギャップと近い値が得られた．これらの結果からRh₂/AlPO₄においても静電ポテンシャルの埋め込みにより，dangling bondの影響はFermi準位付近の電子状態に残るものの，大きく改善されていることが示された．

Figure 6.

 Frontier orbitals of deformed surfaces of Rh₂/Al₂O₃ and Rh₂/AlPO₄ for (A) slab model, (B) bare cluster model, and (C) embedded cluster models. Numbers of point charges are 1060 (24 × 28 × 15 ų), 1.079 × 10⁵(250 × 250 × 15 ų), and 1.452 × 10⁶(900 × 900 × 15 ų) for the Al₂O₃ embedded models with very small, middle, large, and very large number of point charges (VSP, MP, and VLP), respectively. Those are 2934 (53 × 50 × 15 ų), 7.551 × 10⁵(270 × 250 × 15 ų), and 1.016 × 10⁶(970 × 920 × 15 ų) for the AlPO₄'s with VSP, MP, and VLP, respectively. PE means periodic electrostatic potential. (Reprinted with permission from M. Matsui and S. Sakaki, J. Phys. Chem. C 2017, 121, 20242–20253. Copyright (2018) American Chemical Society.)

Table 3  Frontier orbital energies (ε_HOMO and ε_LUMO relative to vacuum level;^a eV) and band gaps of deformed Al₂O₃ and AlPO₄ surfaces. (Reprinted with permission from M. Matsui and S. Sakaki, J. Phys. Chem. C 2017, 121, 20242–20253. Copyright (2018) American Chemical Society.)

^aIn the slab model, the electrostatic potential at the middle point between two surfaces is taken as a standard (vacuum level).^bThese geometries were taken to be the same as the corresponding moiety of Rh₂/Al₂O₃ and Rh₂/AlPO₄ optimized by the slab calculations. ^cNumbers of point charges are 1060 (24 × 28 × 15 ų), 11940 (83 × 83 × 15 ų), 1.079 × 10⁵(250 × 250 × 15 ų), 5.879 × 10⁵(590 × 590 × 15 ų), and 1.452 × 10⁶(920 × 920 × 15 ų) for the Al₂O₃ embedded models with very small, small, middle, large, and very large number of point charges (VSP, SP, MP, LP, and VLP), respectively. Those are 2934 (53 × 50 × 15 ų), 8310 (88 × 84 × 15 ų), 7.551 × 10⁵(270 × 250 × 15 ų), 4.115 × 10⁵(620 × 580 × 15 ų), and 1.016 × 10⁶(970 × 920 × 15 ų) for the AlPO₄'s with VSP, SP, MP, LP, and VLP, respectively. ^dPE indicates periodic electrostatic potential. ^eBand gap = ε_LUMO − ε_HOMO unless no caution is presented. ^fB3LYP calculation cannot be performed with the slab model in VASP. ^gFrontier orbitals close to HOMO or LUMO.^hSCF calculation cannot be converged in PBE calculations. ⁱHOMO−2. ^jε_LUMO − ε_HOMO–2. ^kLUMO+1. ^lHOMO−4.

^mε_LUMO+1 − ε_HOMO–4.

金属–表面相互作用の解析には，フロンティア軌道の寄与するFermi準位付近の状態密度が重要となる．スラブモデルとPEモデルのFermi準位付近のDOSと，LUMOバンドへのp-DOSをFigure 7に示した．スラブモデルにおいては，前章で示したたように，Rh₂/AlPO₄ではLUMOバンドのp-DOSはほぼ全てが価電子バンドのDOSに現れ (Figure 7A-2)，Rh₂/Al₂O₃では価電子，導電バンドの両方に現れている (Figure 7A-1)．Rh₂/(AlPO₄)₁₅，Rh₂/(Al₂O₃)₁₂のPEモデルにおいても，いずれもFermi準位付近においてスラブモデルと同じような分布となった (各々Figure 7B-2, 7B-1)．以上のように，今回開発した静電ポテンシャル埋め込みクラスターモデルは，スラブモデルの電荷移動の挙動を再現可能なものであり，金属–表面相互作用を議論するうえで有用なものであることが示された．

Figure 7.

 DOSs of the Rh₂/Al₂O₃ and Rh₂/AlPO₄, p-DOSs of HOMO and LUMO of the deformed Al₂O₃ and AlPO₄ surfaces of (A) slab model, (B) embedded cluster model with the periodic electrostatic potential (PE model). Black, red, and blue lines represent total DOS, p-DOS onto the LUMO band of the deformed Al₂O₃ or AlPO₄ surface, and p-DOS onto the HOMO band of the deformed Al₂O₃ or AlPO₄ surface. Vertical pink line represents the Fermi level. The vacuum level is set to be zero. No level is found in the cyan stripe in the case of slab model. The negative density of states corresponds to the density of β spin states. (Adapted with permission from M. Matsui and S. Sakaki, J. Phys. Chem. C 2017, 121, 20242–20253. Copyright (2018) American Chemical Society.)

一方で，スラブモデルと埋め込みクラスターモデルのDOSには違いも見られる．Fermi準位より高エネルギー部分 (Figure 7の図中に背景色で示した領域) においては，スラブモデルではDOSが存在しないが，PEモデルではDOSが存在する．これは，クラスター末端領域に存在する軌道によるものである．このように，今回の開発で行ったような長距離静電相互作用の考慮のみでは，クラスターモデルのdangling bondの影響を完全に取り除くことはできない．この影響は，Fermi準位付近にd軌道の関与する複雑な電子状態の存在する遷移金属酸化物担体を対象にした際に，より深刻な問題を引き起こすことが懸念される．この問題への対応方法としては，外部ポテンシャルにパウリ反発などの量子的な効果を有効的に取り入れることがこれまでに検討されている．先に紹介したPEECMにおいては，クラスター周囲に存在するカチオンサイトに基底を用いない有効擬ポテンシャルを配置することで，強すぎる外部点電荷ポテンシャルにクラスターの電子状態が束縛されることを抑制している [59]．また，DFETでは，電荷密度に基づく埋め込みポテンシャルを汎関数微分を利用して求めるより高度で複雑な手法を開発しており，現在も発展中である [55]．このように，担持金属触媒のための埋め込みクラスターモデルの開発には，今回行った表面の長距離静電ポテンシャルの考慮の他に，さらに量子的な影響を有効的に取り入れることが，課題として存在している．

4 まとめ

金属微粒子と担体表面との界面の電子状態が関与し，触媒活性の決定に重要な役割を担う担持金属触媒の金属–表面相互作用の本質を理解・予測するためには，金属微粒子の軌道と無限周期的な表面バンドとの相互作用の電子論的な解析が不可欠であり，また，巨視的な表面全体の電子状態を求めながら，同時に微視的な触媒反応の進行する金属微粒子近傍領域の電子状態を高精度に記述する第一原理計算手法が要求される．本総説で紹介した軌道–バンド間相互作用に基づくp-DOSによる解析手法と，長距離静電ポテンシャルを考慮した埋め込みクラスターモデルは，これまでに述べたように，そのような解析手法，計算手法として有用，有望なものである． 一方で，軌道–バンド間相互作用に基づく解析手法には，電荷移動以外の静電，分極，交換反発といった相互作用への適用や定量的評価がなされていない，また周期的静電ポテンシャル埋め込みクラスターモデルには，長距離静電ポテンシャルの考慮だけではdangling bondの影響を完全には取り除けないなど，まだ課題が存在する．しかし，相互作用解析に関しては分子科学分野で開発された手法のアイディアを取り入れる，埋め込みクラスターモデルにおいては有効的に量子効果を外部ポテンシャルに取り入れるなど，課題を克服するための改善方法はあり，我々も現在検討を進めているところである．このような解析手法と計算手法の発展によって，金属–表面相互作用を利用した貴金属フリー担持金属触媒の開発設計において，理論研究がさらに重要な寄与を示すことが期待される．

謝辞

本総説で紹介した研究は，文部科学省の委託 (元素戦略拠点形成型プロジェクト) と，JSPS科研費JP17K05750の助成を受けて実施したものである．また，計算の一部は自然科学研究機構・計算科学研究センターのスーパーコンピュータを利用して行った．

参考文献

• [1]   M. A. Newton, Chem. Soc. Rev., 37, 2644 (2008). doi:10.1039/b707746g19020678
• [2]   C. T. Campbell, Nat. Chem., 4, 597 (2012). doi:10.1038/nchem.141222824888
• [3]   M. Machida, K. Murakami, S. Hinokuma, K. Uemura, K. Ikeue, M. Matsuda, M. Chai, Y. Nakahara, T. Sato, Chem. Mater., 21, 1796 (2009). DOI:10.1021/cm9005844
• [4]   K. Ikeue, K. Murakami, S. Hinokuma, K. Uemura, D. Zhang, M. Machida, Bull. Chem. Soc. Jpn., 83, 291 (2010). doi:10.1246/bcsj.20090256
• [5]   M. Machida, S. Minami, K. Ikeue, S. Hinokuma, Y. Nagao, T. Sato, Y. Nakahara, Chem. Mater., 26, 5799 (2014). doi:10.1021/cm503061g
• [6]   M. Machida, S. Minami, S. Hinokuma, H. Yoshida, Y. Nagao, T. Sato, Y. Nakahara, J. Phys. Chem. C, 119, 373 (2015). doi:10.1021/jp509649r
• [7]   M. Machida, T. Eidome, S. Minami, H. P. Buwono, S. Hinokuma, Y. Nagao, Y. Nakahara, J. Phys. Chem. C, 119, 11653 (2015). doi:10.1021/acs.jpcc.5b01846
• [8]   J. H. Kwak, J. Hu, D. Mei, C. W. Yi, D. H. Kim, C. H. F. Peden, L. F. Allard, J. Szanyi, Science, 325, 1670 (2009). doi:10.1126/science.117674519779194
• [9]   A. Bruix, J. A. Rodriguez, P. J. Ramírez, S. D. Senanayake, J. Evans, J. B. Park, D. Stacchiola, P. Liu, J. Hrbek, F. Illas, J. Am. Chem. Soc., 134, 8968 (2012). doi:10.1021/ja302070k22563752
• [10]   R. S. Mulliken, J. Am. Chem. Soc., 74, 811 (1952). doi:10.1021/ja01123a067
• [11]   G. H. Booth, A. Grüneis, G. Kresse, A. Alavi, Nature, 493, 365 (2012). doi:10.1038/nature1177023254929
• [12]   C. Müller, B. Paulus, Phys. Chem. Chem. Phys., 14, 7605 (2012). doi:10.1039/c2cp24020c22373600
• [13]   L. Maschio, D. Usvyat, F. R. Manby, S. Casassa, C. Pisani, M. Schütz, Phys. Rev. B, 76, 075101 (2007). doi:10.1103/PhysRevB.76.075101
• [14]   S. Casassa, M. Halo, L. Maschio, C. Roetti, C. Pisani, Theor. Chem. Acc., 117, 781 (2007). doi:10.1007/s00214-006-0198-x
• [15]   C. Pisani, M. Schütz, S. Casassa, D. Usvyat, L. Maschio, M. Lorenz, A. Erba, Phys. Chem. Chem. Phys., 14, 7615 (2012). doi:10.1039/c2cp23927b22334044
• [16]   P. Y. Ayala, K. N. Kudin, G. E. Scuseria, J. Chem. Phys., 115, 9698 (2001). doi:10.1063/1.1414369
• [17]   D. Usvyat, B. Civalleri, L. Maschio, R. Dovesi, C. Pisani, and M. Schütz, J. Chem. Phys., 134, 214105 (2011). PMID:21663342 DOI:10.1063/1.3595514
• [18]   M. Marsman, A. Grüneis, J. Paier, G. Kresse, J. Chem. Phys., 130, 184103 (2009). doi:10.1063/1.312624919449904
• [19]   A. Grüneis, G. H. Booth, M. Marsman, J. Spencer, A. Alavi, G. Kresse, J. Chem. Theory Comput., 7, 2780 (2011). doi:10.1021/ct200263g26605469
• [20]   S. J. Nolan, M. J. Gillan, D. Alfè, N. L. Allan, F. R. Manby, Phys. Rev. B, 80, 165109 (2009). doi:10.1103/PhysRevB.80.165109
• [21]   S. J. Nolan, P. J. Bygrave, N. L. Allan, F. R. Manby, J. Phys. Condens. Matter, 22, 074201 (2010). doi:10.1088/0953-8984/22/7/07420121386379
• [22]   T. Shiozaki, S. Hirata, J. Chem. Phys., 132, 151101 (2010). doi:10.1063/1.339607920423161
• [23]   M. Matsui, J. Phys. Chem. C, 118, 19294 (2014). doi:10.1021/jp5047084
• [24]   M. Matsui, M. Machida, S. Sakaki, J. Phys. Chem. C, 119, 19752 (2015). doi:10.1021/acs.jpcc.5b02691
• [25]   M. Matsui, S. Sakaki, J. Phys. Chem. C, 121, 20242 (2017). doi:10.1021/acs.jpcc.7b04126
• [26]   R. F. W. Bader, Atoms in Molecules: A Quantum Theory (Oxford University Press), 1990.
• [27]   G. Henkelman, A. Arnaldsson, H. Jónsson, Comput. Mater. Sci., 36, 354 (2006). doi:10.1016/j.commatsci.2005.04.010
• [28]   E. Sanville, S. D. Kenny, R. Smith, G. Henkelman, J. Comput. Chem., 28, 899 (2007). doi:10.1002/jcc.2057517238168
• [29]   W. Tang, E. Sanville, G. Henkelman, J. Phys. Condens. Matter, 21, 084204 (2009). doi:10.1088/0953-8984/21/8/08420421817356
• [30]   M. Yu, D. R. Trinkle, J. Chem. Phys., 134, 064111 (2011). doi:10.1063/1.355371621322665
• [31]   R. S. Mulliken, J. Chem. Phys., 23, 1833 (1955). doi:10.1063/1.1740588
• [32]   T. Hughbanks, R. Hoffmann, J. Am. Chem. Soc., 105, 3528 (1983). doi:10.1021/ja00349a027
• [33]   R. Dronskowski, P. E. Blöchl, J. Phys. Chem., 97, 8617 (1993). doi:10.1021/j100135a014
• [34]   D. Sańchez-Portal, E. Artacho, J. M. Soler, Solid State Commun., 95, 685 (1995). doi:10.1016/0038-1098(95)00341-X
• [35]   D. Sánchez-Portal, E. Artacho, J. M. Soler, J. Phys. Condens. Matter, 8, 3859 (1996). doi:10.1088/0953-8984/8/21/012
• [36]   M. D. Segall, C. J. Pickard, R. Shah, M. C. Payne, Mol. Phys., 89, 571 (1996). doi:10.1080/002689796173912
• [37]   M. D. Segall, R. Shah, C. J. Pickard, M. C. Payne, Phys. Rev. B, 54, 16317 (1996). doi:10.1103/PhysRevB.54.16317
• [38]   V. L. Deringer, A. L. Tchougréeff, R. Dronskowski, J. Phys. Chem. A, 115, 5461 (2011). doi:10.1021/jp202489s21548594
• [39]   N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, D. Vanderbilt, Rev. Mod. Phys., 84, 1419 (2012). doi:10.1103/RevModPhys.84.1419
• [40]   S. Dapprich, G. Frenking, J. Phys. Chem., 99, 9352 (1995). doi:10.1021/j100023a009
• [41]   H. A. Graetsch, Acta Crystallogr. C, 57, 665 (2001). doi:10.1107/S010827010100513311408662
• [42]   M. Digne, P. Sautet, P. Raybaud, P. Euzen, H. Toulhoat, J. Catal., 226, 54 (2004). doi:10.1016/j.jcat.2004.04.020
• [43]   X. R. Shi, D. S. Sholl, J. Phys. Chem. C, 116, 10623 (2012). doi:10.1021/jp301114n
• [44]   G. Kresse, J. Hafner, Phys. Rev. B, 47, 558 (1993). doi:10.1103/PhysRevB.47.558
• [45]   G. Kresse, J. Furthmüller, Comput. Mater. Sci., 6, 15 (1996). DOI:10.1016/0927-0256(96)00008-0
• [46]   G. Kresse, J. Furthmüller, Phys. Rev. B, 54, 11169 (1996). doi:10.1103/PhysRevB.54.11169
• [47]   F. Libisch, C. Huang, E. A. Carter, Acc. Chem. Res., 47, 2768 (2014). doi:10.1021/ar500086h24873211
• [48]   N. Govind, Y. A. Wang, A. J. R. da Silva, E. A. Carter, Chem. Phys. Lett., 295, 129 (1998). doi:10.1016/S0009-2614(98)00939-7
• [49]   N. Govind, Y. A. Wang, E. A. Carter, J. Chem. Phys., 110, 7677 (1999). doi:10.1063/1.478679
• [50]   P. Huang, E. A. Carter, J. Chem. Phys., 125, 084102 (2006). doi:10.1063/1.233642816964996
• [51]   C. Huang, M. Pavone, E. A. Carter, J. Chem. Phys., 134, 154110 (2011). doi:10.1063/1.357751621513378
• [52]   F. Libisch, C. Huang, P. Liao, M. Pavone, E. A. Carter, Phys. Rev. Lett., 109, 198303 (2012). doi:10.1103/PhysRevLett.109.19830323215432
• [53]   S. Mukherjee, F. Libisch, N. Large, O. Neumann, L. V. Brown, J. Cheng, J. B. Lassiter, E. A. Carter, P. Nordlander, N. J. Halas, Nano Lett., 13, 240 (2013). doi:10.1021/nl303940z23194158
• [54]   K. Yu, F. Libisch, E. A. Carter, J. Chem. Phys., 143, 102806 (2015). doi:10.1063/1.492226026373999
• [55]   K. Yu, E. A. Carter, Proc. Natl. Acad. Sci. USA, 114, E10861 (2017). doi:10.1073/pnas.171261111429203675
• [56]   J. D. Goodpaster, T. A. Barnes, F. R. Manby, T. F. Miller, III, J. Chem. Phys., 137, 224113 (2012). doi:10.1063/1.477022623248993
• [57]   F. R. Manby, M. Stella, J. D. Goodpaster, T. F. Miller, III, J. Chem. Theory Comput., 8, 2564 (2012). doi:10.1021/ct300544e22904692
• [58]   G. Knizia, G. K. L. Chan, Phys. Rev. Lett., 109, 186404 (2012). doi:10.1103/PhysRevLett.109.18640423215304
• [59]   A. M. Burow, M. Sierka, J. Döbler, J. Sauer, J. Chem. Phys., 130, 174710 (2009). doi:10.1063/1.312352719425801
• [60]   N. A. Richter, S. Sicolo, S. V. Levchenko, J. Sauer, M. Scheffler, Phys. Rev. Lett., 111, 045502 (2013). doi:10.1103/PhysRevLett.111.04550223931382
• [61]   C. Tuma, J. Sauer, Chem. Phys. Lett., 387, 388 (2004). DOI:10.1016/j.cplett.2004.02.056
• [62]   C. Tuma, J. Sauer, Phys. Chem. Chem. Phys., 8, 3955 (2006). doi:10.1039/B608262A17028686
• [63]   S. Tosoni, J. Sauer, Phys. Chem. Chem. Phys., 12, 14330 (2010). doi:10.1039/c0cp01261k20886145
• [64]   A. D. Boese, J. Sauer, Phys. Chem. Chem. Phys., 15, 16481 (2013). doi:10.1039/c3cp52321g23949344
• [65]   A. Mazheika, S. V. Levchenko, J. Phys. Chem. C, 120, 26934 (2016). doi:10.1021/acs.jpcc.6b09505
• [66]   K. N. Kudin, G. E. Scuseria, Chem. Phys. Lett., 283, 61 (1998). doi:10.1016/S0009-2614(97)01329-8
• [67]   K. N. Kudin, G. E. Scuseria, J. Chem. Phys., 121, 2886 (2004). doi:10.1063/1.177163415291598
• [68]   M. W. Schmidt, K. K. Baldridge, J. A. Boatz, S. T. Elbert, M. S. Gordon, J. H. Jensen, S. Koseki, N. Matsunaga, K. A. Nguyen, S. Su, T. L. Windus, M. Dupuis, J. A. Montgomery, J. Comput. Chem., 14, 1347 (1993). doi:10.1002/jcc.540141112
• [69]   M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, T. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman, D. J. Fox, Gaussian, Inc., Wallingford CT, 2016.