# Fileset [MDR_202407_SS.pdf](https://mdr.nims.go.jp/filesets/477a993d-a6b3-4017-90d1-9c8331c13bb0/download) ## Creator [末原 茂](https://orcid.org/0000-0001-7423-2830) ## Rights © 2024 無機マテリアル学会[In Copyright](http://rightsstatements.org/vocab/InC/1.0/) ## Other metadata [材料研究のツールとして第一原理計算を使うために](https://mdr.nims.go.jp/datasets/9f4a71ed-d33b-42fb-9871-eb8e60738779) ## Fulltext 材料研究のツールとして第⼀原理計算を使うために* 末原 茂** １ はじめに 実⽤材料の研究開発では、様々な環境下での特性予測や評価の⼀部に，実験に基づく経験的なパラメータを⽤いた数値シミュレーシ ョンが使われている．⼀⽅、⽐較的単純な組成の物質・材料研究では，基礎的な物理法則（原理）に基づく「第⼀原理電⼦状態計算」（第⼀原理計算）1),2)を⽤いた“物質⾃体”のシミュレーションから，材料の特性を導き出すことも珍しくなくなってきている．従来，第⼀原理計算は，固体電⼦論や理論化学を学んだ研究者が，理論計算の専⾨家として⾏うものであった．この計算法の中⾝をキチンと知るためには，量⼦⼒学に基礎を置いた量⼦化学や固体物理学の内容を理解し，さらに電⼦状態の数値的な解き⽅まで学ぶ必要があったからである．電⼦状態計算のプロフェッショナルを⽬指すのであれば，これら全てを学ぶことに時間を掛けても惜しくないだろう．しかし，こうしたシミュレーションの結果に期待を寄せ，とにかく計算を実⾏して，⾃分たちの研究に役⽴てられるかどうかを試してみたいという研究技術者も多いはずである．幸いなことに，現在，有償・無償を問わず，世界的に広く利⽤されている第⼀原理計算プログラムの多くは，⾮常に使い勝⼿が良くなっている（VASP, Quantum-ESPRESSO, Gaussian, 他）3)~5). 今後は，計算法の理論的な背景の理解は“そこそこ”に⽌め，数ある研究開発ツールの１つとして，第⼀原理計算をブラックボックス的に利⽤していくスタイルが主流となっていくと予想される．実験装置の使い⽅の要点だけを押さえておけば，電気回路や機械部品の詳しい働きや仕組みを知らずとも，正しく実験できるのと同様に，専⾨家でなくとも第⼀原理計算を使った研究に取り組める時代になってきたということである．以下では，そのようなスタイルで計算研究に臨む場合にポイントとなる部分を，筆者が過去に取り扱った無機物質（マイカ，層状ホウ素，TaO3 ナノメッシュ，テルライト光学材料）の計算例と共に解説する．⽬を三⾓にしなくても読み進められるように，可能な限り数式を⽤いないようにした．表現に不正確さや不適切さを感じる読者もいるかもしれないが，ご容赦いただきたい． ２ 第⼀原理計算のアウトライン 2.1 第⼀原理計算は（も）近似計算である 厳密な第⼀原理計算法では，故意に実験結果を再現するような仕組みは完全に排除されている．⾔い換えれば，⽤いられる計算式の係数などの決定に実験データは関わっていない．そのため“⾮経験的”な電⼦状態計算法と呼ばれることもある．ただし，“第⼀原理”，“⾮経験的”（欧⽂では“first-principles”, “ab initio”, “non-empirical”などのワードが使われる）という響きから，例えば「第⼀原理計算は前提のない完璧な理論を採⽤しており，計測誤差がなければ，実験値は計算結果と⼀致するはずだ」などと期待を膨らませてはいけない．物質・材料研究分野で⽤いられる第⼀原理計算では，⼀部の物理法則しか採⽤しておらず（例えば，⼀番⾝近な万有引⼒の法則は含まれていない），電⼦状態を求めるために使われる理論も⼤胆な近似そのものである．こうした理論レベルでの近似による現実と差異は，計算環境（メモリや記憶媒体の量，CPU 速度など）の都合で，物質の構造モデルをやむを得ず⼩さくしたことに起因する誤差とは別の種類のものである．つまり，第⼀原理計算は根っこの部分から，すでに複数の近似が導⼊された計算法である．対象とする系にもよるが，より現実に近い結果を求める場合には，第⼀原理計算に多少の経験的な補正を加えるケースも珍しくない．この場合，得られる結果は，第⼀原理計算の仕組みを使ったものではあるが，厳密な意味では第⼀原理計算のそれではない．このような計算を半経験的（semi-empirical）な計算と呼ぶ場合があるが，第⼀原理計算プログラムをツールの１つとして利⽤するのであれば，その計算が厳密な第⼀原理計算か半経験的な計算かということに拘る必要はなく，⽬的とするシミュレーションに相応しい⽅を選択すればよい．半経験的な計算による結果であっても，第⼀原理計算と同様の⼿続きで得られた電⼦状態を解析すれば，電⼦論の⽴場から材料特性や機能発現のメカニズムに深くメスを⼊れることは可能である． 2.2 物質を核と電⼦の混合物として取り扱う 物質⾃体のシミュレーションを⾏う第⼀原理計算では，対象とする物質を正負の荷電粒⼦の混合物とみなす．正負の粒⼦を，それぞれ“原⼦核”と“電⼦”に対応させると，全電⼦第⼀原理計算法と呼ばれる⽅法になるが，ここで，内殻電⼦は原⼦間の結合に重要な関与はしていないとみなして，内殻電⼦を原⼦核に押し込めた“イオン芯”と，結合に関与する“価電⼦”の組み合わせで近似することも多い．この⽅法は，偽物の原⼦（擬原⼦; pseudo atom）の核の静電場を使うという意味で，擬ポテンシャル（pseudo potential）法と呼ばれる．各元素の擬ポテンシャルは，全電⼦第⼀原理計算によって求められた孤⽴原⼦が⽰す価電⼦のエネルギー準位や，価電⼦に対するポテンシャルを参照して作られる．参照先が実験値ではなく，全電⼦第⼀原理計算の結果であるため，この⼿順で作成された擬ポテンシャルは第⼀原理的（⾮経験的）に作成されたものである．ただし，分⼦や固体の中にある原⼦は，原⼦ 1 個の状態とは異なるイオン価数を持つこともあるため，作成時に参照する全電⼦計算の原⼦の価数をいくつに設定するか，あるいは，近似精度を上げるために価電⼦状態の影響を受けやすい浅い準位の内殻電⼦（semi-core electron という）を価電⼦として陽に扱うかどうかなどは，その作成者に委ねられている．これは，第⼀原理的に作成された擬ポテンシャルでも，作成条件が異なれば，結果も多少異なるということを意味するが，昨今の計算プログラムに付随する擬ポテンシャルは，価電⼦数の変化（即ち，イオン価数の変化）に対して，全電⼦第⼀原理計算の結果を⼀定の誤差範囲内で再現するように調整して作られているため，その影響は限定的である（擬ポテンシャルに“Transferabilityがある”と表現する）．内殻電⼦の計算を省略できる第⼀原理擬ポテンシャル法は，計算コストの⾯で⾮常に有利であり，多くの原⼦が集まった巨⼤分⼦やクラスター，ユニットセルの⼤きな固体，液体などでは標準的に使われる 6)． 通常，対象となる物質の構造モデルは，真空中（図１a）かユニットセル内（図１b）に配置される．前者は，主に分⼦やクラスターを対象とした第⼀原理“分⼦軌道法”で⽤いられ，後者は固体や液体などの凝縮系を対象とした第⼀原理“バンド計算法”で⽤いられる．対象となる系やシミュレーションに求める精度にもよるが，特に後者のユニットセルを⽤いたバンド計算では，“端”や“表⾯”のない理想的な凝縮系だけでなく，真空中に浮いた分⼦や表⾯構造などを対象としたシミュレーションなども⾏うことができる．例えば，⼀辺 30Å程度の⼤きな⽴⽅体ユニットセル内に，原⼦や⼩さな分⼦モデルを配置したバンド計算では，ほぼ孤⽴した原⼦・分⼦のシミュレーション結果が得られる（ただし，図１a のような完全孤⽴系モデルを⽤いた分⼦軌道計算と⽐べると，精度的には劣る）．また，厚さ 12〜20Å程度の真空層を，板状の物質構造モデルに重ねてユニットセル内に配置すれば，積層⽅向に真空境界をもった表⾯構造モデルができあがる．これは表⾯・界⾯のシミュレーション研究で頻繁に⽤いられるスラブモデル（slab model）である． 2.3 第⼀原理計算では原⼦配置から全エネルギーと電⼦状態を求める 最もシンプルな第⼀原理計算は，⼊⼒された原⼦配置での全エネルギーと電⼦状態の計算である．冒頭に挙げた第⼀原理計算プログラムをはじめ，広く使われるプログラムパッケージでは，多種多様なシミュレーションを⾏うことができるが，静⽌した原⼦配置での系の全エネルギーと電⼦状態を求める計算が，それ以降に⾏われる全てのシミュレーションの核となる計算である．この計算は１つの固定された原⼦配置で⾏われる計算という意味で“シングルポイント”の計算と呼ばれる． 第⼀原理計算では，核系（イオン芯または原⼦核）は古典⼒学に，電⼦系は量⼦⼒学に従うとして，別々に計算が進められる．ここで，“別々に計算が進められる”の意味は，核の状態と電⼦の状態を分離して取り扱う近似を採⽤しているということである（ボルン-オッペンハイマー近似 7）という）．通常の核の運動に対して，電⼦は⼗分な速度で応答し，すぐに新たな電⼦状態を形成できるので，電⼦状態の計算時には，核は静⽌していると考えて良いという⾒⽅である．求められた全エネルギーは，⼤まかな理解としては，正負の荷電粒⼦の運動エネルギーとポテンシャルエネルギーの合計であるが，電⼦が波動性も備えた量⼦⼒学的な粒⼦であり，さらに磁性に関わる“スピン”という⾃由度も持っているため，電⼦に関するエネルギーの計算は⼀筋縄ではいかず，そこには⼤きな近似が⼊ってくる．困難さの原因は，電⼦間相互作⽤（電⼦間反発）の取り扱いであるが 8)~11)，ここでは深⼊りせず，とりあえずは，以下に説明する⼀電⼦近似の概要だけを知っておけば⼗分と考える． * First-principles calculation as a Materials research Tool ** Shigeru SUEHARA 2.4 ⼀電⼦近似を⽤いて電⼦状態を求める 多電⼦系の電⼦状態を求めるには，いわゆる“多体問題”を解かねばならないが，直接解くことが困難であるため，結論として，１つ１つの電⼦（𝑖）に注⽬した⼀電⼦の基礎⽅程式： ℎ𝜙! =ε!𝜙! (ℎ：⼀電⼦のハミルトニアン）を解き，そこで得られる⼀電⼦エネルギー（ε!）と⼀電⼦軌道関数（𝜙!）を使って，近似的に系の電⼦状態を得る．これを“⼀電⼦近似”と呼ぶ（１つの電⼦につき，１つのエネルギーと１つの軌道が対応するようなイメージとなるが，あくまでこの近似の下での描像である)． ここで，⼀電⼦の基礎⽅程式は⼀つではなく，“２流派”あることに触れておく．先の基礎⽅程式の形は変わらないが，⼀電⼦ハミルトニアン（ℎ）の内容に少し異なりがある（電⼦相関研究の専⾨家には⼤きな違いとなる部分である）．１つは，多電⼦のシュレーディンガー⽅程式を直接⼀電⼦近似で解いていく議論の中で⽣まれた“ハートリー・フォック(Hartree-Fock; HF)⽅程式”である 8)．HF法は多電⼦原⼦や分⼦で適⽤・検証され，分⼦軌道法の中核をなす近似法として⼤きな成功を収めた．もう１つの流派は，密度汎関数理論(Density Functional Theory; DFT)を，実⽤⾯から⽀えるコーン・シャム(Kohn-Sham; KS)⽅程式である 9)-11)．DFT では，波動関数ではなく，電⼦密度の観点から “真の多電⼦状態”を追いかけて議論が展開される．その中で⽣まれた KS 法で⽤いられる局所密度近似（Local Density Approximation; LDA）は，特に固体の電⼦状態計算（バンド計算法）を⾶躍的に進歩させた．以降，これらの基礎⽅程式をベースに，⼀電⼦近似の⽋点を補うような様々な近似法が開発・展開されている 11)-13)． 第⼀原理分⼦軌道計算プログラムでは，“理論レベル”2）,3)，第⼀原理バンド計算法では，“交換相関汎関数(exchange-correlation function)”4)-6)というキーワードを⽤いて，⼀電⼦近似に使う基礎⽅程式の内容（電⼦間相互作⽤の近似部分）などを指定する．この選択は，“第⼀原理的”に決まるものではなく，ユーザーが⾏うものであることに注意しよう（つまり，様々な近似の第⼀原理計算法があるということである）． 2.5 多数の基底関数を使って電⼦分布を表現する 第⼀原理計算によって求められるものは，全エネルギーと電⼦状態である．エネルギーは“数値”であるが，電⼦状態を表現する⼀電⼦軌道，あるいは，そこから導かれる電⼦密度分布は“関数”である．⼀般に，各原⼦に付随する電⼦の分布は，原⼦が単体で存在していた場合と，その原⼦が分⼦や固体などの物質を構成する⼀部となった場合で異なっているはずである．また，同じ物質でも，電場や磁場の中に置かれた場合とそうでない場合とで，その分布は異なるだろう．第⼀原理計算では，複数の単純な関数（基底関数）に係数をつけて，その和をとる格好で⼀電⼦軌道関数を近似し（図２），結果として様々な形の電⼦分布を表現する．これは，任意の関数を，⾼い変形⾃由度を持つテイラー級数展開式やフーリエ級数展開式で表現するのと同様のアプローチである．たとえば，第⼀原理分⼦軌道計算プログラムで筆頭に挙げられる Gaussianパッケージ 2）,3)では，ガウス型関数が基底関数として⽤いられる（パッケージ名の由来となっている）．そこでは，ひとまず原⼦軌道関数（s, p, d,…）をガウス型関数で展開近似し，さらに各原⼦位置に配置した原⼦軌道関数に係数をつけて⾜し合わせ，⼀電⼦軌道（分⼦軌道）関数をつくる（分⼦軌道からみれば，原⼦軌道関数も基底関数とみなせる）．この展開を⽤いる⽅法を LCAO-MO 法と呼ぶ．また，周期構造をもつ固体などを対象とした第⼀原理バンド計算法では，様々な周期（振動数）をもつ sin関数と cos関数を基底関数として係数付きで⾜し合わせて，⼀電⼦軌道（バンド軌道と呼ぶ）を表現する．この⽅法は，平⾯波法（Plane Wave法; PW法）と呼ばれる 6),10)． プログラムは，基礎⽅程式を満たすまで⼀電⼦軌道関数の形を変化させながら演算を繰り返す．⼀電⼦軌道関数を展開近似する１つのメリットは，解が得られるまで何度も繰り返し⾏われる積分演算などが，素性の良い基底関数のおかげで，短時間で簡単に済むことである．また，展開する基底関数の数を制限することで，計算精度と計算速度の調整をユーザーが⾏うことも可能になる．例えば，研究初期段階では，少ない数の基底関数を使った計算で，⼤まかな結果を素早く知り，次のステップで，数多くの基底関数を使⽤した⾼コストな⾼精度の計算を⾏うといった戦略をとることが可能となる．基底関数を多く使⽤する場合，分⼦軌道法では「⼤きな基底関数系(基底関数セット)を使う」と表現する．PW法の場合は，基底関数の数が，“カットオフエネルギー”という値で間接的に指定されるので「⾼いカットオフエネルギーを使う」などと表現する 6)． 2.6電⼦状態の近似精度と基底関数系の⼤きさ 繰り返しになるが，第⼀原理計算で得られる電⼦状態は，⼀電⼦近似のもとで解かれた⼀電⼦エネルギーと⼀電⼦軌道で近似的に表現される．ここで仮に“真の電⼦状態”がわかっていたとしよう．これまでの話から，真の電⼦状態と第⼀原理計算で得られた“電⼦状態の近似解”との差は，（１）基礎⽅程式で切り捨てられている電⼦間相互作⽤が原因となっている部分と，（２）基底関数が有限であることによる近似の悪さに起因する部分の２つから⽣じていると考えてよい．どのような電⼦の分布も表現可能な，理想的な（無限⼤の）基底関数系を⽤意して計算できたとすれば，後者（２）の差分は無くなるはずだが，適⽤された電⼦間相互作⽤の近似による差分（１）は必ず残る．これは，近似解が真の電⼦状態にどうしても届かない部分の差として認識しておくべき点である． ⼀⽅，どのような理論レベルを⽤いても，⼩さな基底関数系を⽤いれば，（２）による差分が⼤きくなり，貧しい結果しか得られない．基底関数系は⼤きいほど良いが，その⼤⼩で，精度と計算時間のトレードオフがあるため，無闇に⼤きくできない．そこで，本格的なシミュレーション計算に臨む前に，これ以上基底関数系を⼤きくしても，殆ど全エネルギーが変化せず，“⼀定値に収束した”と判断できる基底関数系の⼤きさ，すなわち，（２）の差が事実上なくなる基底関数系の条件をあらかじめ調べておくとよい．具体的には，徐々に基底関数系を⼤きくしながら（PW 法であれば，カットオフエネルギーを増加させながら），対象物質の構造は変えずに計算を⾏い，全エネルギーの変化を⾒ていく．この作業により，⾼精度計算で⽤いる基底関数の⼤きさの上限を決めることができる．また，基底関数系の⼤きさと全エネルギー関係を定量的に知ることで，⼩さめの基底関数系の選択も⾃信をもって⾏うことができるようになる． 第⼀原理計算のアウトラインは以上である．以下では，筆者の過去の研究を例にとりながら，ツールとして第⼀原理計算を使う際に，押さえておいた⽅がよいポイントのいくつかを紹介する． ３ 第⼀原理計算シミュレーションの実際 3.1 構造最適化計算 静⽌した１つの原⼦配置で第⼀原理計算を⾏うと，系の全エネルギーと電⼦状態がわかる（シングルポイントの計算）．次に，核の位置を少し移動させたシングルポイントの計算を⾏う．エネルギーが減った配置について同様の作業を繰り返していくと，最終的に初期構造に⼀番近い安定構造に辿り着く．各核に働く⼒はプログラムが出してくれるので，それに従って原⼦配置を変化させていけば効率的である．これは“構造最適化”と呼ばれる⼿続きであり，プログラムに構造最適化の指⽰をすれば，⾃動的に⾏ってくれる．また，周期構造をもつ固体の場合には，ユニットセルの格⼦定数を変化させながら，求めた全エネルギー値をプロットしていくと，エネルギーの極⼩点で格⼦定数が，さらにプロット曲線から圧縮率がわかる 1)．表１に，そのようにして求めた⾦雲⺟（マイカ）の第⼀原理バンド計算による格⼦定数を⽰した 15)．交換相関汎関数に PBE(Perdew-Burke-Ernzerhof による近似)を⽤いた計算結果は，実験値と良い⼀致を⽰しているが，van der Waals (vdW)⼒補正を⽤いた半経験的な計算結果は，さらに実験値に近づいている．PBE は，第⼀原理計算において最も頻繁に利⽤される交換相関汎関数の１つであるが，vdW ⼒のような弱い⼒の源となる相互作⽤成分が取りこぼされていることがわかっており，そのため僅かだがこのような差があらわれていると考えられる．当たり前だが，その構造を形成するために vdW ⼒が重要な役割をもつ物質では，vdW ⼒が考慮された計算を⾏わないと実際の構造を精度よく再現できない（例えば，グラファイトや，層間に⽔分⼦層が存在するような粘⼟物質，あるいは分⼦性の液状物質など）．vdW ⼒を考慮した⽅法にも経験的なものから，より第⼀原理に近づいたものまで多くの種類がある 16)．表１に⽰したマイカのケースでは，vdW⼒補正は必須ではないと考えても良いが，同じ構造で，層間がほぼvdW⼒でしか結合していないパイロフィライト（葉蝋⽯）やタルク（滑⽯）などでは，vdW⼒を考慮した計算が必須となってくる 17). 全エネルギーが原⼦位置の変化に対して極⼩となる点（安定構造）がわかると，フォノン（格⼦振動）の計算が可能となり，そこからさらにフォノン⽐熱（格⼦⽐熱）を求めることができる 18)．第⼀原理計算を熱⼒学計算に結びつける⽅法の詳細は，他稿に譲るが 19)，表１に⽰した第⼀原理計算（PBE）によるエントロピーと⽐熱の値が，実験値と近い値になっていることに，少し驚きを覚えないだろうか．特定の特性を，実験データを基に予測して得る経験的なシミュレーションではなく，物理法則だけを頼りにした“物質⾃体”のシミュレーションから出発し，材料の特性を導き出している１つの例である． 3.2 第⼀原理分⼦動⼒学計算 構造最適化を⾏う際に使った核にかかる⼒を，そのまま運動⽅程式にのせれば，原⼦が熱⼒学的に指定された条件下で動き回るような分⼦動⼒学（Molecular Dynamics; MD）シミュレーションを⾏うことができる（第⼀原理分⼦動⼒学法；以下第⼀原理MD法と略す）．図 3に，層状ホウ素の構造探索時に⾏った第⼀原理 MD 計算の結果を⽰す 20)．対象物質は，現在，ボラフェン（ボロフェン）として知られている物質の１つである．研究当初，表⾯構造や表⾯組成が明らかではなかったことから，ホウ素被覆率の異なる構造モデルを複数⽤意して MD 計算を⾏った．具体的には，2.2 節で述べたスラブモデルを採⽤し，⼤きなユニットセル内に，板状の⾦属基板（ZrB2）と真空層を重ねて配置，基板上にホウ素原⼦を均⼀に並べて初期構造モデルとしている（図３左列）．そして，実験を模した温度条件で⼀定時間保持後，冷却するような温度制御をしたMD計算によって最終構造を得ている（図３右列）．アニール実験を模したこの⼀連の MD 計算により，点線で囲まれた７つのホウ素で形成される表⾯ホウ素特有の新しい短距離秩序がわかり（図 3a, 3bの右列の点線領域），これをヒントに最終的に層状ホウ素の新しい構造を発⾒できた 20)． この研究内容はさておき，⼀般に MD 計算での温度調節は，系内の核の速度分布が指定した温度に合致するように，核の動く速度をプログラムが加減することで⾏っている 21）．そして，タイムステップ（dt）毎に，新しい原⼦配置が⽣成され，その配置で第⼀原理計算（シングルポイントの計算）が⾏われる．dt はユーザーが指定しなければならない値であるが，ここで，原⼦が動き回る結果を早くみたいからと，dt を⻑くしてしまうと，速度の⼤きな原⼦が，次のステップで遠くに⾏き過ぎてしまい MD 計算が破綻する可能性が出てくる．⼀⽅，精度の良い軌跡を得ようとして dt を必要以上に短くしてしまうと，単位時間あたりの第⼀原理計算回数が増加し，計算効率が下がる．安定した MD 計算を効率よく⾏うには，どういう dt を選べばよいだろうか．そのヒントとして，系内で⼀番早く動く原⼦の持つ振動を再現できる程度の dt を設けるというやり⽅がある．図４に振動周期を 8 分割した dt を設定した場合における第⼀原理計算のタイミングと，速度の更新の様⼦を⽰した（dt としてはかなり⻑めである）．例えば，⽔素を含むセラミクス材料の多くは OH結合が⼀番⾼い振動数（3700 cm-1程度）である．その周期は⾮常に短く，第⼀原理MD計算では，通常 dt は 0.5 fs から 1 fs程度の値が採⽤される．逆に⽔素を含まない系であれば，dt を⻑くとることも可能であるが，核の運動速度が上昇すると想定される⾼温 MD などでは，やはり dt を少し短めに設定する⽅が安全である．上に紹介した層状ホウ素の MD 計算では，最速の振動は系内で⼀番軽い元素同⼠の B-B 結合（周期は 28fs 程度）であると予想されたため，dt はその1/10 程度に設定して結果を得ていた 20)． 3.3 化学反応や拡散現象の追跡 化学反応過程や原⼦・イオンの拡散過程を，上述の第⼀原理 MD 計算を使ってシミュレーションできるだろうか．“構造が変化する”，“原⼦が移動する”という点においては，どちらも前節の研究例と同じであるため，MD計算でシミュレーションできそうに思える．しかし，そのような期待は裏切られる場合が多い．以下でその理由を説明する． 図５に⽰すような２つの状態 A，Bが，活性化エネルギーΔE で隔てられているとした場合，室温で A 状態から B 状態に系（構造）が移る⼤雑把な待ち時間は，表２のとおりである．第⼀原理 MD 計算がシミュレーションできる時間は，最近の⼤規模並列計算機やワークステーションを使⽤し，計算負荷を減らすために構造モデルを少し⼩さめに調整したとしても，せいぜい~100 psぐらいまでが限界である．どんなに頑張っても到底 1秒にも満たないのが，現在（2023年）の相場であろう．表２に⽰した A 状態から B 状態への系の変化の待ち時間は，単純な理論（⼀次元調和遷移状態理論 1)）に基づいた⼤変雑な⾒積りとなるが，例えば，活性化エネルギーとして１eV が必要な化学反応や拡散の場合，室温でそれが起こるのは，運が良くて数時間（〜10 時間）に 1 回あるかないかであり，この事象を 1秒のシミュレーションもできない第⼀原理 MD 計算によって観測しようとするのは⾮現実的である．また仮に⻑時間の第⼀原理 MD 計算ができたとしても，その期間内に必ずその事象が起きるとは限らず，終始同じような軌跡を眺めるだけの徒労に終わってしまう可能性も⾼い．温度条件にもよるが，⼀般に数 eV のエネルギーバリアで隔てられた状態間の往来は，通常の第⼀原理 MD 計算では，滅多に起きることのない事象（レアイベント）であり，その詳細を調べることは困難である．レアイベントを MD 計算で観測する⽅法として，A 状態から B状態への変化を強制的，あるいは，そのように促すような仕組みを取り⼊れた⾼度なMD計算法も存在するが 12), 21) ,22)，計算機に⾼負荷をかける第⼀原理計算を dt 毎に⾏いながら軌跡を求める点では通常の MD 計算と違いはないため，やはり⼀定の計算コストを覚悟しておかなければならない． ⼀⽅，MD計算に拘らなければ，ナッジド・エラスティック・バンド（Nudged Elastic Band; NEB）法 1), 23)という，明⽰的に反応経路や拡散経路を指⽰して探索する⼿法を選択できるケースもある．この⽅法では，A 状態から B 状態への構造変化をある程度予想して，中間構造モデル（中間イメージと呼ぶ）をユーザー⾃⾝が⽤意し，プログラムが，その中間イメージを仮想的なゴムで繋ぎ合わせて（全体を連ねたまま）構造最適化を⾏うようなことをする．この仕組みのおかげで，構造最適化中に各中間イメージがより安定な A 状態や B 状態に落ちてしまうのを避けることができる．NEB 計算で得られた最終的な中間イメージの連なりは， A 状態から B状態への最もエネルギーの低い構造変化の経路であり，最⼩エネルギー経路（Minimum Energy Path; MEP）と呼ばれる．MEP は反応経路や拡散経路の重要な候補となる．NEB法は，初期構造と最終構造が判明している場合にしか使えないが，適⽤できれば，活性化エネルギーや遷移状態構造も同時に得ることができる．計算中にいつ反応や拡散が起こるかわからない MD シミュレーションを実⾏するよりも⼤変効率的に系の変化を調べることができる．以下に，NEB 計算の例として，TaO3 ナノメッシュのイオン透過性を調べた研究を紹介する． TaO3 ナノメッシュは，⽳の空いたシート状セラミックスであり，様々な応⽤が期待される新物質である（図６a）．NEB計算を使って，この⽳をイオン半径の異なる１価のイオン（Li+, Na+, K+）が，表側から裏側に拡散移動可能かどうかを調べたみた結果，MEP エネルギーは図６b のようになった 24）．拡散移動できるかどうかは，TaO3 ナノメッシュの表⾯から裏⾯に抜ける際のエネルギーバリアの様⼦から議論できる．例えば，Li+イオンの透過パス（MEP）は，図６a の点線 1-5（1 は表側，5は裏側，2-4は中間イメージ）となるが，透過パスの中に⼤きなエネルギーバリアは⾒られず（図６b），表と裏を⾏ったり来たりできそうである．Na+イオンと K+イオンの場合は，中間イメージ３（各イオンが⽳のほぼ中央に位置したときの構造）で，エネルギーの極⼤値があらわれ，拡散移動のエネルギー障壁となっている．エネルギーバリアの値から，Li+＞Na+>K+の順番で，この⽳を通りやすいことがわかるが，透過イベントの起き る頻度を HTST（表２） に 基 づ い て考えて み る と ，Li+>Na+>>>K+と表現すべき程に K+イオンは通らないという印象になる． さて，もしこのイオンの拡散移動現象を，NEB 計算ではなく室温のMD計算でシミュレーションしたらどうなるだろうか．以下の考察では，各イオンが表⾯から離れていかないような仕掛けのある MD 計算を思考する．まず，K+イオンと Na+イオンについてであるが，⼤きなエネルギー障壁があるため，その拡散現象はレアイベントである．時間をかけて室温の第⼀原理 MD 計算を続けても，同じような状態がひたすら繰り返され，拡散現象を観測できる可能性は極めて低い．⼀⽅，Li+イオンについては，エネルギー障壁がほぼ無いため，もしかしたら，表から裏へ通り抜ける現象を，第⼀原理 MD 計算で観測できるかもしれない．実験シミュレーションで⾏われる加速試験と同様に，少し温度を⾼めにした MD シミュレーションを⾏えば，観測できる確率も少しは上げられるはずである．その理由は，先に述べた HTST（表２）から，事象の起こる頻度が活性化エネルギーと系の温度の⽐で決まっていることからもわかる．ただし，MD計算をしている限り，⽬的とするイベントを観測できるかどうかは，やはり計算の成り⾏き次第であることには違いない．つまり，この考察から学べることは，活性化エネルギー（エネルギー障壁）が⾮常に低いと期待できない限り，最初から第⼀原理MD計算を使って化学反応や拡散の過程を追跡しようとするのは⽌めておいた⽅が良さそうだということである． 第⼀原理 MD 計算は，“仮想実験”の重要な柱の１つに位置付けられるものではあるが，計算負荷が⼤きく，取り扱える原⼦数やシミュレーション時間などに事実上の制約が多い．そのため，詳細不明なイベントやダイナミクスを追う場合には，MD計算以外の別のアプローチも選択肢に⼊れておくべきである． 3.4 基底関数の⼤きさの重要性 第⼀原理計算で得られる“全エネルギー”だけで議論できることは多い．例えば，活性化エネルギーについては前述のとおりであるが，単に化学反応エネルギー（エネルギー収⽀）だけを知りたければ，反応前後の各系の構造モデルを⼊⼒し，それぞれの構造最適化後に，系間の全エネルギーの差を求めればよい．それぞれの構造モデルで，フォノン⽐熱の要素を加えた全エネルギーの⽐較を⾏えば，標準状態（有限温度）での議論も可能である 15), 19), 20）．また，異なる原⼦配置の全エネルギーの差分だけでなく，“物質⾃体”がもつ物性値の多くが，全エネルギーの変化（ストレス，電場，磁場などに対するエネルギーの変化）で与えられることを考えると，全エネルギーの数値精度は⾮常に重要である．2.6節で述べた通り，⼀度使⽤する電⼦間相互作⽤の近似理論を決めた場合には，基底関数系の⼤きさが，得られる結果の精度を決定する．その重要性がわかる例を以下に紹介する 25）． 表３に，基底関数系の⼤きさによる光学特性の理論計算値の違いを⽰した．計算対象となった物質は，⾮線形光学材料のテルライト（TeO2）系化合物である．この物質群では，⻑短２種類のTe−O 結合を持つTeO4構造単位が，⾼い分極率（α）と２次⾮線形光学効果（γ）を⽣み出している．ここで⾔う⾮線形光学材料とは，光を⼊射した際に，異なる波⻑の光を発⽣する物質である．古典電磁気学的な描像としては，核系の重⼼電荷と電⼦系の重⼼電荷で作られるマクロなダイポールモーメントが，⼊射光の電場に僅かに歪んで応答することで，この⾮線形光学効果があらわれる．電気回路で例えると，質の悪いアンプ（線形性が良くないアンプ）に交流信号を⼊⼒すると⾼調波が発⽣する現象と同じである．分⼦軌道計算法を使って⾏ったこの計算では，TeO4 構造単位モデル（クラスター）に強さの異なる電場をかけ，得られたダイポールモーメントの電場に対する変化からαとγを算出している．表３には，SBKJCと呼ばれる基底関数セットをベースに，d型， sp型（d型とは異なる形をもつ基底関数），先の d型とは異なる d型などの関数を加えた基底関数系で計算した結果が並んでいる．表３（最下⾏）の SBKJC(3d, sp)という⼤きな基底関数系を使った計算結果が，この計算例の中では⼀番精度の⾼いものであり，最も全エネルギーが低い．補強されていない SBKJC 基底関数セットでの計算結果（表３最上⾏の結果）は，SBKJC(3d,sp)基底関数セットを使った結果に⽐べ，⻑短２種類あるはずの Te-O 結合が再現できておらず，αは２割弱低い．さらにγに⾄っては 1/3 程度とかけ離れて低い値になっている．これは電場の中に置かれた TeO4構造内のわずかな電⼦分布の歪みを SBKJC 基底関数系が⼩さすぎて精度良く再現できていないからである．また，γの値がαより基底関数系の⼤きさに敏感であるのは，αが電場に対するダイポールモーメントの１階微分（全エネルギーの 2階微分）なのに対して，γは 3階微分（全エネルギーの 4 階微分）だからであろう．以上から，この計算例では，構造のみに注⽬した場合は SBKJC(d )基底関数セットで⼗分であるが，αやγをある程度精度良く求めるのであれば，少なくとも d型と sp 型の関数を補強した SBKJC(d,sp)の基底関数セットを使うべきであることがわかる． ⼀般的な⼿続きとして，構造を次々変化させていく研究初期の構造最適化（やMD計算の初期段階）では，計算効率の⾯から標準的な基底関数系か，少しだけ補強した基底関数系を⽤いて研究を先に進めるとよい．そして，異なる系の同⼠のエネルギーの⽐較や物性値などを算出する段階で，⼤きな基底関数系を⽤いた⾼精度計算を⾏い，最終結果を得ると同時に計算精度の確認を⾏うと効率的である． 3.5 光電⼦スペクトルと⼀電⼦エネルギー 実験で電⼦状態を調べると⾔った場合，その筆頭に挙げられるのは，X 線や紫外線などの光で試料の電⼦を叩き出して（イオン化して），その電⼦の束縛エネルギーを調べる光電⼦分光法であろう．特に前者は，表⾯化学の分析法として “ESCA（Electron Spectroscopy for Chemical Analysis）”という名称でも知られている．この実験スペクトルを，第⼀原理計算を使ってシミュレーションする場合は，束縛されている各電⼦についてイオン化前後の系の全エネルギーの差を求めればよいが，⼀電⼦近似の理論式上でその差分を求めると，HF近似では，なんとイオン化前の基底状態の⼀電⼦エネルギーが，そのままその電⼦のイオン化エネルギーになっている（Koopmans の定理）．⼀⽅，KS近似のもとでは，イオン化エネルギーは対象とする⼀電⼦軌道から電⼦を半分取り除いた格好で得られた⼀電⼦軌道エネルギーの値で近似される（Janak の定理）．つまり，⼀電⼦近似の基礎⽅程式の格好は同じであっても，その中⾝が違えば，基底状態の⼀電⼦エネルギーの意味・解釈は異なっている 10）．しかし，いずれの近似においても（HF 近似，KS 近似，あるいはその発展系も含めて），イオン化前の基底状態の⼀電⼦エネルギーの並びが，光電⼦スペクトルと⽐較されることは多い．本来別のものを並べて⽐較していることを知った上で，参考程度に並べて議論しているのである．ただし，結合状態やイオン価数の違いによる内殻電⼦の化学シフトを精度よく求めて議論したい場合は，イオン化過程を考慮した計算⼿続きで⾒積る必要は出てくる 10), 26）． 3.6 電荷密度解析と化学結合 第⼀原理計算で得られた電⼦状態を解析すると，原⼦間の結合がどの原⼦軌道からつくられているか，イオン価数はいくつか等を定量的に評価できる．すると，共有結合やイオン結合，σ結合，π結合，混成軌道（例えば sp3 ハイブリッド軌道など）等といった馴染み深い“化学の⾔葉”を使いながら，その物質を解釈できるようになる．こうした化学的な描像を得るには，原⼦間に跨って分布している電⼦を個々の原⼦や軌道成分に帰属させる必要がある．この作業は，“電荷密度解析”（あるいは電⼦密度解析）と呼ばれ，その代表的なものには，Mulliken 法 28），Löwdin 法 29)がある（電荷だけの帰属であれば Bader 法 30)も選択できる）．⼀電⼦軌道に原⼦軌道様の関数を⽤いる分⼦軌道計算では電荷密度解析に時間は掛からないが，⾮局在型の基底関数(sin 関数,cos 関数)を使⽤している PW法では，原⼦軌道を⽤いた電⼦分布の表現に変換する計算⼿続きを経る必要がある 27）． 第⼀原理計算では，基底関数の数が多いほど電⼦状態をより正確に計算できることは，既に述べた．そのため，⾼精度計算で得られた電⼦状態の電荷密度解析を⾏うと，１つの結合に多くの軌道成分が関わってくることも珍しいことではない．例えば，sp3 混成軌道型の結合をもつとされる分⼦の電荷密度解析を⾏うと，d軌道も混成していたという場合もあるかもしれない．また，イオン価数も従来の化学結合モデルと符合しないことはよくあることである（イオン結合性とされる結晶でも電荷密度解析を⾏うと，正負イオンともに，想定以上に中性側に寄っている場合も多い）．もし，従来の化学結合モデルと，第⼀原理計算の電荷密度解析の結果に差異があると感じられた場合には，それをどう捉えるべきだろうか．第⼀原理計算が与えてくれるのは，物理法則を基礎にして算出した電⼦状態に関わるデータであり，その計算内訳も明快である．電荷密度解析も，その⼿法の選択により値が多少は異なるが，定量的な値を提供してくれる．しかし，第⼀原理計算も，その結果を使った電荷密度解析も，そこまでのツールであり，解釈まで与えてくれるものではない．⾔い換えれば，計算結果の最終的な解釈はユーザーに委ねられている．例えば，計算結果を⽤いて従来の化学結合モデルを再解釈し，その過程で部分的に新しい解釈を与えるようなこともできるかもしれないが，計算結果の詳細にまで化学的な視点での解釈を追い求めても意味がないような場合もある．⽐較対象として並べて良いものではないかもしれないが，化学結合という描像も，第⼀原理計算も，いわば近似モデルである．それぞれが違った視点から物質を眺めているので，共通する解釈ができる部分とそうでない部分があって当たり前と考えるべきだろう． 4 おわりに 先⼈たちの多⼤な努⼒により，⾮常に⾼度な計算技術が凝縮された第⼀原理計算プログラムが，今や⼿軽に利⽤できる時代となった．インターネットで検索すれば，様々な種類の第⼀原理計算プログラムや，そのインストール⽅法の解説が⾒つかる．本稿を読んで，第⼀原理計算をツールとして使ってみようと思われた読者がいたなら，次にすることは，計算法を知るためにテキストを開いて勉強し直すことではなく，とにかくプログラムをインストールして実際に計算を実⾏し，第⼀原理計算のルーチンを⼀度体験することである．材料開発の現場で使われるツールの１つに第⼀原理計算が加わることで，何かが変わる可能性は⼗分あると信じている． 本稿の作成にあたって，NIMS の⼩林⼀昭博⼠，新井正男博⼠，佐久間博博⼠，法政⼤学の渡邊雄⼆郎教授に貴重なご意⾒を頂きました．深く感謝致します． 参考⽂献 1) ショール，ステッケル（佐々⽊泰造，末原茂 共訳）, “密度汎函数理論⼊⾨”, 吉岡書店 (2014). 2) Foresman, Frisch（川内進 訳），”電⼦構造論による化学の探究”,第 3 版, Gaussian, Inc. (2017). 3) Gaussian : https://gaussian.com (2023.10.20) 4) VASP: https://www.vasp.at (2023.10.20) 5) Quantum-ESPRESSO : https://www.quantum-espresso.org (2023.10.20) 6) ⼩林⼀昭，“徹底解剖 第⼀原理計算（連載）”, ⾦属, 79 ,931 (2009). 7) M. Born, J. R. Oppenheimer, Ann. Phys. 84, 457(1927). 8) ザボ，オストランド（⼤野公男, 阪井健男, 望⽉祐志 共訳）, “新しい量⼦化学（上・下）”, 東京⼤学出版会 (1987). 9) 川添良幸, 三上益弘, ⼤野かおる, “コンピュータ・シミュレーションによる物質科学 分⼦動⼒学とモンテカルロ法”, 共⽴出版 (1996). 10) R. M. Martin, “Electronic Structure, Basic Theory and Practical Methods”, Cambridge Univ. Press. (2004). 11) J. Kohanoff, “Electronic Structure Calculations for Solids and Molecules”, Cambridge Univ. Press (2006). 12) D. Marx, J Hutter, “Ab Initio Molecular Dynamics ‒ Basic Theory and Advanced Methods”, Cambridge Univ. Press (2009). 13) ⾚井久純，⽩井光雲, “密度汎関数法の発展−マテリアルデザインへの応⽤”,丸善出版 (2012). 14) 岡崎誠，“べんりな変分原理”, 共⽴出版 (1993). 15) S. Suehara, H. Yamada, Geochim. Cosmochim. Acta, 109, 62-73 (2013). 16) I. Hamada, Phys. Rev. B 89, 121103 (2014). 17) H. Sakuma, S. Suehara, J. Geophys. Res.: Solid Earth, 120, 2212-2219 (2015). 18) キッテル(宇野良清，津屋昇，新関駒⼆郎，森⽥章，⼭下次郎 共訳), “キッテル固体物理学⼊⾨”, 第 8 版，丸善(2008), 第５章. 19) 末原茂, ⾦属, 80 (2010) pp. 11-17. 20) S. Suehara, T. Aizawa, T. Sasaki, Phys. Rev. B 81, 085423 (2010). 21) D. Frenkel, B. Smit, “Understanding MOLECULAR DYNAMICS from Algorithms to Applications”, 2nd Ed., Academic Press (2002). 22) 岡崎進，吉井範⾏，“コンピュータ・シミュレーションの基礎”，第２版，化学同⼈ (2011). 23) D. Sheppard, P. Xiao, W. Chemelewski, D. Jonson, J. Chem. Phys. 136, 074103 (2012). 24) C. Wang, N. Sakai, Y. Ebina, S. Suehara, T. Kikuchi, D. Tang, R. Ma and T. Sasaki, J. Mater. Chem. A, 10, 25481 (2022). 25) S. Suehara, P. Thomas, A. P. Mirgorodsky, T. Merle-Mejean, J. C. Champarnaud-Mesjard, T. Aizawa, S. Hishita, S. Todoroki, T. Konishi, and S. Inoue, Phys. Rev. B 70, 205121 (2004). 26) T. Aizawa, S. Suehara, S. Hishita, S. Otani, M. Arai, Phys. Rev. B 71, 165405 (2005). 27) P. C. Muller, Ch. Ertural, J. Hempelmann, R. Dronskowski, J. Phys. Chem. C, 125, 7959 (2021). 28) R. S. Mulliken, Int. J. Chem. Phys. 23,1833 (1955). 29) P.-O. Löwdin, Advances in Quantum Chemistry, 5, 185 (1970). 30) R. F. Bader, Acc. Chem. Res. 18, 9, (1985). 筆者紹介 末原 茂（すえはら しげる） 国⽴研究開発法⼈物質・材料研究機構（NIMS）電⼦・光機能材料研究センター資源循環材料グループ主幹研究員 博⼠（⼯学） 1990年早稲⽥⼤学理⼯学部材料⼯学科卒，同年科学技術庁無機材質研究所（現NIMS）に⼊所後，SPCTS(現ヨーロピアンセラミックスセンター; IRCER), フランス国⽴科学研究センターCNRS/リモージュ⼤学客員研究員(2002-2003)，現在に⾄る． 連絡先 〒305-0044茨城県つくば市並⽊ 1丁⽬ 1番地（勤務先） 表１ ⾦雲⺟ K2Mg6(OH)4[(Si6Al2)O20] の格⼦定数，エントロピー，⽐熱15)計算値(PBE)計算値（vdW⼒補正） 実験値格⼦定数a (Å) 5.3428 5.3128___ 5.3158b (Å) 9.2537 9.2012___ 9.2036c (Å) 10.5006 10.3164___ 10.3100α(°) 90.00 90.00___ 90.00β(°) 99.80 100.13___ 99.89γ(°) 90.00 90.00___ 90.00エントロピー(J・mol-1・K-1) 311.19 315.9⽐熱(J・mol-1・K-1) 352.88 355.1* ⼀次元調和遷移状態理論(HTST)による図５のA状態からB状態へ のホッピング頻度は下式で表される1)．𝑘!→# 𝑠$% = 𝜈 exp −𝛥𝐸𝑘#𝑇ν の範囲: 1012~1013 (0.1-1ps 原⼦の振動周期)ΔE (eV) kA→B (s-1) 待ち時間(1/kA→B)10 < 10-154 > 10145 年1 2x10-4 〜 2x10-5 < ~10時間 0.1 2x1010 〜 2x1011 < ~0.1 ns0.01 7x1011 〜 7x1012 < ~1 ps表２ HTST*による室温での反応頻度（kBT~26meV）表３ TeO4構造単位モデルの計算結果23)基底関数系 Energy（Ha)α(a.u.)γ(a.u.)Te-O結合⻑(Å)（⻑/短）SBKJC -74.120 50.3 3500 2.04/2.04SBKJC(sp) -74.138 59.3 11328 2.04/2.04SBKJC(d ) -74.227 49.8 3231 2.03/1.96SBKJC(d,sp) -74.240 57.4 11257 2.03/1.97SBKJC(2d,sp) -74.265 59.4 11333 2.02/1.96SBKJC(3d,sp) -74.280 60.2 11858 2.03/1.96※ 数値を⽐較しやすくするために原⼦単系(Ha, a.u.)を⽤いた．Ha（Energy) ： 4.3597×10-18 J ; 27.2114 eVa.u. (α） ： 1.6488×10-41 C2・m2・J-1 ; 1.8621×10-30 m3a.u. (γ） ： 6.2354×10-65 C4・m4・J-3 ; 7.0423×10-54 m5・V-2図１ 物質は核と電⼦の混合物第⼀原理計算内部では，物質は核（⽩抜きの丸記号）と電⼦（⿊点記号）の混合物として扱われる．その構造モデルは，(a)単体で真空中に浮かべるか，(b)点線のような境界で作られるユニットセル内に配置される．(b) ユニットセル内に構造モデル（固体，液体など）を置く（a) 真空中に構造モデル（分⼦，クラスターなど）を置く図２ ⼀電⼦軌道の基底関数展開近似のイメージ通常，基底関数には，s, p, d, ...などの原⼦軌道関数の系列（a）か，平⾯波関数（sin関数,cos関数）(b)が選ばれる．各々の基底関数の係数（c1,c2,...）の値を変化させれば，⼀電⼦軌道関数𝜙!は，展開された範囲内ではあるが，様々な分布を表現することが可能である．プログラムは基礎⽅程式を満たすこれらの係数を変分法で求めている．𝜙! 〜 c1 + c2 + c3 + ...𝜙! 〜 c1 + c2 + c3 + ... (b) 平⾯波の線型結合（a) 原⼦軌道型関数の線型結合図３ MD計算によるB/ZrB2(0001)表⾯構造の探索ホウ素の被覆率（monolayer; ML）が異なるZrB2(0001)⾦属表⾯モデルのMD計算前後の表⾯構造を(a)-(d)に⽰した20）．⾒やすくするために⿊線であらわされるユニットセルを2x2倍して表⽰している．MD前（初期構造） MD後（安定構造候補）図4 MD計算における第⼀原理計算のタイミング原⼦の動く軌跡はdt単位の原⼦配置の並びで表現され，dt毎に第⼀原理計算が⾏われる．dt後の原⼦配置は核にかかる⼒や速度ベクトルをもとに⽣成される．dtは少なくとも系の⼀番速い原⼦の振動を再現できる値にする必要がある．核の移動速度が上昇しそうな場合はdtを短めに設定する（⾼温MDや素早く構造が変化するケース）時間振幅（核の平衡位置からのずれ）核の速度：その時刻の原⼦配置での第⼀原理計算dt（振動周期の1/8で設定した場合）図５ 活性化エネルギーで隔てられた２つの状態A x BExΔEA状態 遷移状態 B状態エネルギー反応座標（構造の違い）化学反応・拡散の⽅向短時間MDによる軌跡図６ NEB計算によるイオン透過シミュレーション24)（a）TaO3ナノメッシュモデルとLi+イオンの拡散経路(MEP). （b）Li+, Na+, K+のMEPエネルギー.5 115Ta O H Li+Side viewTop view(a) (b) Energy /eV