113

量子化学計算技術の紹介と今後の展望〈入口紹広紀・工藤孝夫〉

介1

量子化学計算技術の紹介と今後の展望
中央技術研究所 ソリューションセンター
シミュレーショングループ
いりぐち

入口

中央技術研究所 ソリューションセンター
シミュレーショングループ

ひろき

くどう広紀

たかお工藤

1. はじめに孝夫

やソフトウェアの発展により、企業での研究においても活用しやすい状況になってきている。

量子化学計算は分子シミュレーション技術の一つであ

当社では燃料、触媒、潤滑油、機能材などエネルギー・る。分子シミュレーションは化学物質の構造をもとに物理素材に関する研究開発において分子シミュレーション技学の原理に従い計算する技術で、量子化学計算、分子動術の活用を進めている。実験および各種の分析に加えて、力学計算や粗視化分子動力学計算などがある。量子化学

分子シミュレーションを相補的に活用することで実験や分計算は原子や分子の構造や性質を電子状態から解析する

析では捉えることが困難なミクロな領域での現象について手法であり、経験的なパラメータを用いない場合は第一原も知見を拡充することができる。本稿では量子化学計算

理計算とも呼ばれる。分子動力学計算では原子や分子のの理論背景概略、石油精製プロセスで使われる脱硫触媒運動を力学法則に基づいて解析する。粗視化動力学計算についての解析事例を紹介する。

は高分子など分子サイズがより大きな場合に、あるまとまった原子の集団を原子集合体として簡略化して解析する手

2. 量子化学（第一原理）計算の理論背景概略 1 〜 3）法である。各計算手法は図 1 に示すように扱う空間・時間スケールが異なり、解析したい現象に応じて計算手法を量子化学計算はシュレディンガー方程式選択することになる。

（2 − 1）

H は計算する系を表現を議論の出発点とする。このとき＾するハミルトニアンであり、これに対応する電子の固有状
ψ〉と E である。したがっ—態およびそのエネルギー値が—て ψ〉は系を構成する全電子数× 3（スピンは無視している）の次元を持つ領域で定義される波動関数となり、そのままの形で厳密な解を得ることはほぼ不可能に近い。このため、実際に計算機でシュレディンガー方程式を解くためには何らかの近似を導入する必要がある。今日の基本的なアプローチとしてハートリー - フォック法と密度汎関数法（DFT：Density Functional Theory）の二つが挙げられる。
ハートリー - フォック法は電子相関のない系を仮定として置く。このとき波動関数は
図 1 分子シミュレーションの各手法と扱うことのできる領域（2 − 2）をのように個々の電子スピンを含む 1 電子波動関数φ（xij）が規要素としたスレーター行列式の形で表される。φ（xij）

について解析することが可能であり、物質の物性、反応、格直交性を持つという制約を課し、 ψ〉の期待値 E が最—

本稿で紹介する量子化学計算は扱える原子数は限りがあるものの、物性の起因たる要素の一つである電子状態そして物質間の相互作用を詳細に評価することができる。

H を探索するアプローチを小になるようなハミルトニアン＾

物質の特性の本質を知ることができ、メカニズムの理解や取ることによって解くべきシュレディンガー方程式は新規の材料開発に有効な技術と期待される。このようなシミュレーションは多大な計算量のため多くの計算リソースを必要とするが、近年のコンピュータの飛躍的な高性能化− 37 −（2 − 3）
JXTG Technical Review・第60巻第3号（2018年11月）

量子化学計算技術の紹介と今後の展望〈入口広紀・工藤114孝夫〉という形になる（ハートリー - フォック方程式）。r は空間座した触媒を作製した。事実、水素化脱硫反応活性は向上標である。式の左側の各項は順に電子の運動エネルギー、したが、さらなる高活性化を目指すべく現象の知見拡充と外場ポテンシャル、一電子のクーロンポテンシャル、交換相互作用の評価を目的に DFT 計算を行った。

ポテンシャルを表す演算子である。この方程式によって電子の相関作用以外は正確に描写することが可能となった。

3. 1 計算モデル設定

一方で DFT は個々の電子を解く手法を取らず、系の電計算を実行するにあたり必要となるものが、先述の汎関子分布を電子密度 n（r ）によって表現する。具体的には、数・基底関数等の設定に加え、分子の構造である。今回n（r ）によってハミルトニアンのポテンシャル項が一意に定まが求まる。この対ることで固有基底状態波動関数 ψ（i r）
応はホーエンバーグ - コーンの定理により保障される。このときのハミルトニアンの表式はの計算対象は触媒担体およびその表面にある活性金属触媒であるため、担体表面に触媒クラスターが吸着されている状態のモデルを作成した。担体としてγ構造 Al 2O（100）3面とアナターゼ型 TiO2（001）面を用意し、これらを基にTiO2 積層 Al 2O3 モデルも作成した。これら担体モデルに（2 − 4）

担持金属種として Mo 6 S12 クラスターを吸着させたモデルを計算した（図 2）。表面を再現するために周期境界構造と表せる（コーン - シャム方程式）。一見この表式はハートの厚さ方向に 10 Å 程度の真空層を導入し、担体の下半分リー - フォック方程式と同じように見えるが、電子の交換がバルクを再現するように構造を固定して計算を行った。

V XC（r ）はこの定理では概念的に相関を表すポテンシャル＾
V XC（r ）はのみ語られ具体式の形を持たない。ここで ＾（2 − 5）と交換相関エネルギー汎関数 E XC［ （n r）］によって与えられ、DFT は E XC［ n（r ）］を改良していくことで完全かつ具体的に計算可能な形のシュレディンガー方程式を完成させる可Mo 6 S12 クラスターの担体吸着モデル
担体は a）γ構造 Al2O 3（100）面、b）アナターゼ型 TiO 2（001）面、 c）TiO 2 積層 Al2O 3
図2能性を将来に託している。
実際の計算においては以上からもわかる通り電子の交換相関を説明するポテンシャル、すなわち汎関数もしくは計算手法を適切に設定することが肝要となる。特に DFT

以上のモデル計算は DFT 計算パッケージ VASP によっ

では汎関数によって適切な計算が可能な分子系が異なるため、検討したい現象を正しく把握することが求められる。て行われ、交換相関汎関数部分には一般化密度勾配近似

また、ハートリー - フォックないしコーン - シャム方程式を（GGA：Generalized Gradient Approximation）の一種で具体的に解くにあたって波動関数は既知の基底関数の組ある PBE 汎関数（Perdew, Burke, Ernzerhof ）を、基底関数に平面波基底を用いた。そして得られた全電子状態（2 − 6）エネルギーを用いて吸着エネルギーで展開して解く（ローターン方程式）。基底関数の数 m が（3 − 1）多いほど軌道表現は正確となるが、計算数が大きくなっていく。このように、計算手法および基底関数は、精度とコを算出してこれを担体−担 持 金 属 相互作用として評 価ストのトレードオフ関係を考慮して選ぶことになる。

した。ここで E support は担体表面モデルの全エネルギー、
E cluster は担持金属種クラスターの全エネルギー、E total はクラスターが担体表面に吸着したモデルの全エネルギーである。

3. 解析事例：水素化脱硫触媒の担体 - 担持金属相互作用
量子化学計算を活用して従来の研究に提言を行った事
例として、石油精製プロセスに用いられる水素化脱硫（HDS

3. 2 結果
各吸 着モデルにおける吸 着エネルギーの計 算結果を：Hydrodesulfurization）触媒に関する DFT 計算を紹介する 。本検討の動機は高性能な触媒設計指針の確立に

表 1 に 示 す。 まず Al 2O3 と TiO2 の Mo 6 S12 吸 着 モデルあり、触媒の担体 - 担持金属相互作用に着目した。こので比較したところ、Al 2O3 では吸 着エネルギーが -0.0094）相互作用は触媒活性と相関があると言われており、当社中

eV、TiO2 における吸着エネルギーは -0.125 eV となった。

央技術研究所では相互作用を制御するという考えのもとで

Al 2O3 では相互作用が認められない一方で TiO2 では物担体表面に TiO2 薄膜をコーティングし、活性金属を担持理吸着レベルの相互作用がみられる値となった。この相互JXTG Technical Review・第 60巻第3号（ 2018年11月）− 38 −

115

量子化学計算技術の紹介と今後の展望〈入口作用の序列は実験研究者の経験および既往の理論研究5,6）と一致したことで計算の妥当性が示されたと考える。

Al 2O32.7TiO（1layer）/Al 2O32-0.1042.6TiO（2layer）/Al 2O32-0.263

2.4TiO2-0.125

3.2

べると、結晶系とアモルファスで相互作用の劇的な変化は示す可能性が示唆された。アモルファス構造では原子の分布が不均一であるため、相互作用が最大となる箇所が存在するようになったと考えられる。

Adsorption Energy/eV Distance/Å
-0.009孝夫〉

ないものの、アモルファス化で結晶系より強い相互作用を

表 1 各担体の Mo6S12 吸着エネルギーおよび距離の見積値
Support広紀・工藤次に、Al 2O3 に TiO2 を 1 層および 2 層積 層した担体の吸着エネルギーは図 3 のポテンシャルエネルギー曲線（PES：Potential Energy Surface）に示される。ここで、
計 算負荷の都合上本モデルでは構造最適化を行わず担図4

体 - クラスター間距離を変えてそれぞれ一点計算を行い、PES を描くことで最安定な吸着距離とエネルギーを見積a）Mo 6 S12 が吸着した TiO2 1 層積層 Al2O 3モデルと b）その PESもった。この計算方法で得られた見積値は表 1に示される。
このとき、1 層積層モデルでは吸着距離 2.6 Å で吸着エネ

3. 3 考察と指針

ルギー -0.104 eV の最小値をとり、2 層積層モデルでは吸本 研究の計 算で得られた相互作用と触媒活性を考察

着距離 2.4 Å で吸着エネルギー -0.263 eV の最小値をとっする。 担体 種の違いで比 較 すると、 実 験 では Al 2O3 よた。すなわち、2 層積層モデルの吸着エネルギーは TiO2り TiO2 が高い触媒活性を示し、量子化学計算では担持単体における吸着エネルギー -0.125 eV よりも大きな値と金属との相互作用は Al 2O3 より TiO2 が強かった。また、なった。Al 2O3 担体を TiO2 積層数 0、TiO2 担体を TiO2

Al 2O3 に TiO2 をコーティングした担体も Al 2O3 より触媒積層数が十分多いモデルとみなせば、TiO2 積層数と吸着活性は高く 7）、相互作用の計算値は大きかった。さらに、エネルギーの関係に極小値があることがわかり、相互作用

TiO2 コーティングの積層回数で比較すると、薄膜コーティを最大にする適切な積層数の存在が示唆された。

ング回数を増やすごとに反応活性が高まり、相互作用も強まった。すなわち、今回計算したモデルの範囲では担体 担持金属相互作用と触媒活性は正の相関が示された。これは相互作用が弱ければ HDS 活性が高くなるという従来の考察と一見矛盾するように見える。しかしながら、今回の計算モデルは担体原子と担持金属原子が O を介して結合していないこと、および、吸着エネルギーが物理吸着レベルであることから、担体との相互作用が弱いタイプの触媒であると解釈できる。つまり、このタイプの触媒では担体 - 担持金属相互作用が強いほど触媒活性が向上する領域が存在することを示している。以上のことから、ただ相互作用が弱いほど触媒活性が良い訳ではなく、適切な相互作用の強さになるよう制御することが重要であることが分かった。本結果は従来経験則的に行われてきた触媒設

図3計に対して理論的に提言をするものであり、このような解TiO2 積層 Al2O 3 の Mo 6 S12 吸着に関するポテンシャルエネルギー曲線（PES）

析の蓄積が将来の触媒設計に繋がると考える。

以上は結晶系を基にした触媒担体モデルであったが、4. 今後の展望実験で合成される構造はアモルファスであるという要請からアモルファス Al 2O3 に TiO2 を 1 層コーティングした担体

4. 1 多分子系シミュレーションとの連携モデルを第一原理分子動力学法（ab-initio MD）により作このように量子化学計 算は実験で見ることが難しいス

成し、同様に Mo 6 S12 クラスターの吸着計算を行った。得ケールから考察を行い、理論的に現象を説明することでられた構造と PES を図 4 に示す。PES の吸着エネルギー実験考察の深化を促すことができる。しかしながら、物最小値は距離 d = 2.4 Å のときに - 0.126 eV となった。結質の最小単位に基づくシミュレーションゆえに現象の完全晶系 TiO2 1layer-Al 2O3 吸着エネルギーの -0.104 eVと比な再現は難しいことが課題となる。例えば、多量体効果− 39 −JXTG Technical Review・第60巻第3号（2018年11月）

量子化学計算技術の紹介と今後の展望〈入口広紀・工藤116孝夫〉を考察するために高分子の凝集構造をシミュレーションすることを考えると、高分子単体は量子化学計算で構造最書店 , 20094）I riguchi H., Uratani H., Kudo T., Seki H., Yamashita

適化できるだろうが、複数個が独立して存在するモデルにK.,; J. Jpn. Petrol. Inst ., 61, 5, 288 - 293（2018）

なると計算量が大幅に増加して実質計算不可能になるの

5）Joshi, Y.V., Ghosh, P., Daage, M., Delgass, W.N.; J.

は想像に難くない。こうした系を計算するためには図 1 をCatal ., 257, 71（2008）.
6）A rrouvel, C., Breysse, M., Toulhoat, H., Raybaud, P.;
J. Catal ., 232, 161（2005）.
7）S egawa K., Takahashi, K., Satoh, S.; Catal. Today ,63, 123（2000）.

参照すれば分子動力学によって計算することが適切である
ことがわかる。しかしながら、分子動力学は運動方程式に基づくシミュレーションであることから、そのままでは凝集の要因たる相互作用を取り入れることが難しい。ここで量子化学計算によって先に相互作用に関するパラメータを得ることで、それを分子動力学の理論に導入可能なパラメータとすることで双方の良い点を取り入れることができ信頼性が高まる。複数のシミュレーション領域を自在に横断するマルチスケールな解析を行うことによって量子化学計算の活用範囲がより広くなるだろう。
4. 2 機械学習との連携
また、これまでの量子化学計算は特定の実験・現象報告を追従して原因を考察することが多かったが、本来は理論が先行して実験を提案することが理想である。他方、近年では「データ駆動型化学」の確立を提唱するマテリアルズ・インフォマティクス（MI）という概念が注目されている。この分野の理論背景である機械学習・統計解析は一度の計算で多種類のアウトプットを出せる量子化学計算と相性がいいと考えられる。MI という概念を取り込むことによって、計算のみによって新規材料開発を達成することが期待される。さらには強化学習なども活用することによって、これまで実験では考えられなかった領域の構造も探索され、量子化学計算そのものもより効率的になることが期待される。そのためにも、大量の計 算ジョブを自動で実行すること、一度の計算でなるべく多くの情報量を得ることといったワークフローの確立が求められる。
5. おわりに高性能な触媒設計指針の確立を目的に、触媒活性と相関があると言われる担体 - 担持金属相互作用に着目し、量子化学計算を適用した検討事例を紹介した。また、これから求められていくだろう技術分野についても述べた。
当社中央技術研究所では研究活動を加速・革新的技術を開拓する手法として分子シミュレーション技術の向上を目指すものである。
− 引用文献 −1）S z a b o A . , O st lu nd N . S . ; Mo der n Q u a nt u mChemistry: Introduction to Advanced ElectronicStructure Theory , Dover, Minneola, NY, 1996
2）M artin R.M.,;Electronic Structure: Basic Theoryand Practical Methods , Cambridge, UK, 19893）高
 田康民 ;『朝倉物理学体系 15 多体問題特論』, 朝倉
JXTG Technical Review・第 60巻第3号（ 2018年11月）− 40 −