分割統治型密度汎関数強束縛分子動力学 (DC-DFTB-MD)法の最近の展開

Abstract

The authors' group has developed the program named DC-DFTB-K for the on-the-fly quantum mechanical molecular dynamics (MD) simulation of huge systems using the density-functional tight binding (DFTB) method. The combination with the divide-and-conquer (DC) method enables linear-scaling calculation of DFTB energy and its derivatives. Due to the massively parallel implementation, the program can treat systems containing one million atoms on the K computer. In this paper, the recent extension of DC-DFTB-MD technique is outlined together with the illustrative application to chemical reaction dynamics in lithium-ion device.

1 はじめに

大規模かつ複雑な分子系で起こる化学プロセスの制御は，ナノサイズ機能性分子の合理設計や高効率に動作する触媒の開発などにおいて重要な課題である．特に，反応の動的過程を原子・分子レベルで追跡することは，反応経路の探索に有効であると考えられる．反応過程を微視的視点から観測し，精密に解析するための効果的なアプローチとして，量子化学的手法に基づく分子動力学(MD)シミュレーションが挙げられる．たとえば，二次電池材料の開発 [1,2,3]やカーボンナノチューブおよびグラフェンの核生成・成長過程 [4,5]に対して，密度汎関数理論(DFT)やDFTより定式化された密度汎関数強束縛(DFTB)法 [6,7,8]を使用した数百原子を含む系のMD計算などが報告されている．しかし，従来のDFTおよびDFTB計算では，計算コストが系の大きさに対して3乗(もしくはそれ以上)に比例して増加するため，数千~数万原子からなる巨大分子・分子集合体の取り扱いは困難であった．そこで，著者らの研究グループでは，ナノスケール反応系のダイナミクス実現を目指して，分割統治(DC)法 [9,10]をDFTB法に対して適用したDC-DFTB法の理論を構築した．さらに，同手法が超並列計算環境で高効率に動作するプログラムDC-DFTB-Kの開発を進めてきた．本稿では，DC-DFTB法の概要と最近の開発内容，さらに応用例を紹介する．なお，開発した並列計算アルゴリズムや達成された性能の詳細は，今後順を追って論文報告を行う予定である [11]．

2 理論とプログラムの開発

DFTB法のエネルギーは次のように表される．   

$\begin{array}{c}{E}_{\text{DFTB}}=\sum _{\mu \nu }^{\text{AO}}{D}_{\mu \nu }{H}_{\mu \nu }^0+\sum _{A>B}^{\text{atom}}{V}_{AB}^{\text{rep}}\\ +\frac{1}{2}\sum _{AB}^{\text{atom}}{\gamma }_{AB}\Delta q_A\Delta q_B+\frac{1}{3}\sum _{AB}^{\text{atom}}{\Gamma }_{AB}\Delta q_A^2\Delta q_B\end{array}$(1)

ここで，μ, νは原子軌道(AO) ϕの指標，D_μνは密度行列，$H_{\mu \nu }^0=⟨{\phi }_{\mu }|{\hat{H}}^0|{\phi }_{\nu }⟩$はハミルトニアンの電荷ゆらぎに依存しない部分，$V_{AB}^{\text{rep}}$は二体の短距離反発ポテンシャル，Δq_AはMulliken電子密度 [12]から求められる原子Aの誘起電荷，γ_ABおよびΓ_ABは電荷ゆらぎによるクーロン相互作用を適切に記述する関数である．右辺第2, 3, 4項までを含む計算はそれぞれDFTB1 (NCC-DFTB) [6], DFTB2 (SCC-DFTB) [7], DFTB3 [8]と一般に呼ばれる．DFTB2およびDFTB3では，Δq_Aを自己無撞着に決定する(SCC)ため，化学反応シミュレーションにおいて重要な化学結合の生成・開裂の取り扱いにより適したモデルとなっている．式(1)のすべての項は，パラメータから算出できる．積分計算を必要としないため，計算コストはDFTよりも数百倍程度低い．また，パラメータは原子A, Bの種類と原子間距離によって変わり，その値はDFT計算から求めるため，他の半経験的手法と比較して高精度な結果を得ることができる．

DC法は，全系を重なりのない部分系に分割し，部分系周辺の環境効果としてバッファ領域を加えた局在化領域に対する計算を行うことにより，対角化にかかる計算コストを削減する方法である．DC法における全系の密度行列$D_{\mu \nu }^{\text{DC}}$は，部分系αの密度行列$D_{\mu \nu }^{\alpha }$を用いて構築される．   

$\begin{array}{c}{D}_{\mu \nu }\approx D_{\mu \nu }^{\text{DC}}=\sum _{\alpha }^{\text{subsystem}}{p}_{\mu \nu }^{\alpha }{D}_{\mu \nu }^{\alpha }\\ =2\sum _{\alpha }^{\text{subsystem}}{p}_{\mu \nu }^{\alpha }\sum _i^{\text{MO(}\alpha \text{)}}{f}_{\beta }\left({\epsilon }_{\text{F}}-{\epsilon }_i^{\alpha }\right)c_{\mu i}^{\alpha }{c}_{\nu i}^{\alpha \ast }\end{array}$(2)

ここで，$p_{\mu \nu }^{\alpha }$は分割行列，$f_{\beta }\left(x\right)={\left[1+\exp \left(-\beta x\right)\right]}^{-1}$はFermi関数，ε_Fは電子数保存の制約条件から決定される共通のFermi準位である．部分系の軌道エネルギー${\epsilon }_i^{\alpha }$および軌道係数$c_{\mu i}^{\alpha }$は，部分系αの局在化領域について一般化固有値方程式を解くことにより求められる．ε_Fを用いる仕組みによりあらかじめ部分系の電子数を指定する必要がないため，DC法ではシミュレーション中に電荷やスピンの分布が変化し得る場合であっても，MD計算の取り扱いが可能という特長を持つ．

大規模DFTB計算のボトルネックは，ハミルトニアン行列の対角化計算である．この問題を解消したDC-DFTB法は，計算時間の大幅な高速化と，系の大きさに対する線形スケーリングをエネルギー勾配計算も含めて達成している [11]．また，各局在化領域の計算は独立であることから，部分系をMPIにより並列化する処理が可能となる．開発中のDC-DFTB-Kプログラムでは，部分系を計算負荷が均等となるようノードに分散し，各局在化領域内の作業をOpenMPで並列化したハイブリッド並列計算に対応する．著者らは，「京」コンピュータ上で100万原子を超える水溶液モデルの電子状態計算に成功している [11]．

DC-DFTB-Kプログラムで現在利用できるMD計算関連の主要な機能をTable 1に示す．上述したDFTB計算法の他に，スピン分極を考慮した計算(sDFTB) [13]にも対応している．また，弱い相互作用が問題となる系についても取り扱いができるよう種々の経験的分散力補正法を実装した．運動方程式の数値積分には速度Verlet法 [14]を採用し，NVEおよびNVTアンサンブルのMD計算が実行できる．MD計算では，RATTLE法 [15]を用いた距離拘束の動力学，ソフトポテンシャルによる分子集合体の形状の束縛，Lagrange補間を利用した [16] Mulliken電荷の時系列予測によるSCC収束性の改善などの手法がオプションとして選択可能となっている．この他に，構造最適化計算やエネルギー以外のプロパティ(Mulliken電荷，双極子モーメント，Mayer結合次数 [17,18]など)評価が，DC-DFTB法により大規模系に適用可能である．

Table 1.  Basic specification of DC-DFTB-K program.

DFTB energy and gradient	DFTB1/2/3sDFTBShell-resolved SCC-DFTB
Linear-scaling calculation	Divide-and-conquer (DC)
Dispersion correction	Slater-KirkwoodLennard-JonesDFT-D2/3
MD ensemble	NVE, NVT
NVT thermostat	Velocity scalingNosé-Hoover chainBerendsenAndersenLangevin
Additional MD option	RATTLE constraintSoft potentialSimulated annealingLagrange interpolation of Mulliken charges

3 DC-DFTB-MD法の応用：リチウムイオンデバイスの分子シミュレーション

リチウムイオン二次電池やキャパシタは，携帯電子機器用の小型蓄電池や電気自動車用の大型な蓄電池など様々な用途で日常生活に密着している．その性能・信頼性向上のためには，グラファイト負極表面で溶媒分子が還元的に分解して生成する電極界面被膜(SEI)の物性の解明が不可欠である．しかし，電極界面の詳細な反応機構のみならず，SEI膜の構造や組成といった基礎的な性質であっても明らかでない点が多く存在する．そこで，著者らは溶媒・極板・添加剤を含むデバイス全体を取り込んだ化学反応シミュレーションにDC-DFTB-MD法を適用し，SEI膜の構造予測に向けた理論的検討を行ってきた．これまでに，電気二重層の形成に伴う負極界面近傍の電位変化 [19]を報告した．ここでは，負極表面における溶媒の結合生成および分解挙動に関する解析結果を報告する．

本研究で使用したモデル構造をFigure 1左に示す．負極極板は，電解質溶媒が接触する端を-OH, -O, -COOHの置換基で修飾し，それ以外を水素原子でキャップしたグラフェン(炭素原子数は114) 7枚で構成されている．電解質溶媒は，リチウムイオン4原子，EC 24分子，VC 2分子から形成されており，系の全原子数は1364である．時間刻み1.0 fs，温度300 K，電解質溶媒の散逸を防ぐため密度に応じた半球状のポテンシャルを課した条件を用いて，NVTアンサンブルにおけるDFTB2レベルでのシミュレーションを実施した．

Figure 1.

 Left panel: Investigated simulation model of electrolyte decomposition reaction at graphite anode. Right panel: Chemical structure of ethylene carbonate (EC) and vinylene carbonate (VC).

トラジェクトリを解析した結果，大きく分けて(i)溶媒1分子の分解反応，(ii)極板と溶媒間での結合形成反応が確認できた．まず(i)について，代表的なスナップショットと結合長の時間変化をFigure 2に示す．Figure 2上は，C1-O3結合が100 fsで切断した例である．この時，C2-O2結合は一度切断した後，弱い相互作用を持ったまま移動する様子が観察された．Figure 2下は，C1-O2とC1-O3結合が50 fsでほぼ同時に切断し，COが脱離した例である．同様に，2箇所の結合が切断される機構でC₂H₄が発生する分解反応も確認された．COやC₂H₄分子は実験にて検出されており [20,21]，DC-DFTB-MD法は電気分解反応の理論的検討に妥当な手段であることを示唆している．

Figure 2.

 Upper panel: Snapshot of 1 C–O bond cleavage reaction (left) and time course changes of C–O distance (right). Bottom panel: Snapshot of 2 C–O bond cleavages reaction (left) and time course changes of C–O distance (right). Red and green dashed lines are the C–O bond length (1.43 Å) and van der Waals distance (3.22 Å), respectively.

次に，(ii)に関する反応のスナップショットをFigure 3に示す．Figure 3左にある通り，リチウムイオンが溶媒和構造を崩し，溶媒由来の酸素原子と極板に修飾された置換基由来の酸素原子両方と結合している様子が観察された．また，Figure 3右はグラファイト極板と溶媒分子の模式的な反応過程を示している．はじめに，溶媒分子は静電的引力により極板側に水素を向けて接近する．続いて，グラファイト末端にて溶媒分子の水素引き抜き反応が生じる．最後に，水素を引き抜かれた溶媒の炭素原子が極板の炭素・酸素原子と結合を形成する．極板と溶媒間での結合形成反応が生じたトラジェクトリのうち80%以上がこの反応を経由していた．従って，この過程はSEI膜生成反応の第一段階として支配的に進行しているものと考えられる．

Figure 3.

 Left panel: Snapshot of bond formation between Li⁺ and oxygen functional groups of anode. Right panel: Schematic illustration on the bond formation between graphite and electrolyte.

4 まとめ

著者らの研究グループで取り組んできたDC-DFTB法の開発について，最近追加した機能も含めて述べた．また，DC-DFTB-Kプログラムの利用により，通常のPCクラスタを用いても千原子程度の化学反応ダイナミクスがルーチンワークで実行できることを，リチウムイオンデバイス負極表面の電気分解反応を例に示した．我々は，この他に水中のプロトン移動に対する温度依存性の検討や，アミン溶液中の二酸化炭素化学吸収シミュレーションにDC-DFTB-MD法を応用している．後者については，1万原子を超えるモデルのトラジェクトリを得ることに成功しており，実在に近い系の化学反応ダイナミクスが実用的となりつつある．一方，特に周期系では，長距離相互作用計算の効率化がシミュレーションのさらなる高速化に不可欠となることが判明している．今後は，この課題解決に向けたプログラム設計を重点的に進めていく予定である．また，広範囲な領域に対してDC-DFTB法を利用可能とするためのパラメータの拡充を行う必要もある．

DC-DFTB-Kプログラムの使用には，事前相談が必要となる．将来的な公開に向けた準備が計画段階にある．

Acknowledgment

本研究は，文部科学省のHPCI戦略プログラム，および計算物質科学イニシアティブ(CMSI)の支援を受け行われた．プログラム開発の一部は，理化学研究所のスーパーコンピュータ「京」を利用して行った(課題番号：hp140164)．

参考文献

• 1   K. Leung, J. L. Budzien, Phys. Chem. Chem. Phys., 12, 6583 (2010).
• 2   K. Ushirogata, K. Sodeyama, Y. Okuno, Y. Tateyama, J. Am. Chem. Soc., 135, 11967 (2013).
• 3   K. Sodeyama, Y. Yamada, K. Aikawa, A. Yamada, Y. Tateyama, J. Phys. Chem. C, 118, 14091 (2014).
• 4   A. J. Page, Y. Ohta, S. Irle, K. Morokuma, Acc. Chem. Res., 43, 1375 (2010).
• 5   A. J. Page, F. Ding, S. Irle, K. Morokuma, Rep. Prog. Phys., 78, 036501 (2015).
• 6   D. Porezag, T. Frauenheim, T. Köhler, G. Seifert, R. Kaschner, Phys. Rev. B, 51, 12947 (1995).
• 7   M. Elstner, D. Porezag, G. Jungnickel, J. Elsner, M. Haugk, T. Frauenheim, S. Suhai, G. Seifert, Phys. Rev. B, 58, 7260 (1998).
• 8   M. Gaus, Q. Cui, M. Elstner, J. Chem. Theory Comput., 7, 931 (2011).
• 9   W. Yang, T.-S. Lee, J. Chem. Phys., 103, 5674 (1995).
• 10   M. Kobayashi, H. Nakai, in Linear-Scaling Techniques in Computational Chemistry and Physics: Methods and Applications (2011), pp. 97–127.
• 11   H. Nishizawa, Y. Nishimura, M. Kobayashi, S. Irle, H. Nakai, in preparation.
• 12   R. S. Mulliken, J. Chem. Phys., 23, 1833 (1955).
• 13   C. Köhler, G. Seifert, U. Gerstmann, M. Elstner, H. Overhof, T. Frauenheim, Phys. Chem. Chem. Phys., 3, 5109 (2001).
• 14   W. C. Swope, H. C. Andersen, P. H. Berens, K. R. Wilson, J. Chem. Phys., 76, 637 (1982).
• 15   H. C. Andersen, J. Comput. Phys., 52, 24 (1983).
• 16   T. Atsumi, H. Nakai, J. Chem. Phys., 128, 094101 (2008).
• 17   I. Mayer, Chem. Phys. Lett., 97, 270 (1983).
• 18   I. Mayer, J. Comput. Chem., 28, 204 (2007).
• 19   M. Okoshi, H. Nakai, Electrochemistry, 82, 1098 (2014).
• 20   H. Yoshida, T. Fukunaga, T. Hazama, M. Terasaki, M. Mizutani, M. Yamachi, J. Power Sources, 68, 311 (1997).
• 21   M. Onuki, S. Kinoshita, Y. Sakata, M. Yanagidate, Y. Otake, M. Ue, M. Deguchi, J. Electrochem. Soc., 155, A794 (2008).