Enabling Ab-Initio Molecular Dynamics under Bias: The CP2K+SMEAGOL Interface for Integrating Density Functional Theory and Non-Equilibrium Green Functions

技術要約：CP2K+SMEAGOLインターフェースによるバイアス下でのアブイニシオ分子動力学の実現

問題提起
密度汎関数理論（DFT）と非平衡グリーン関数法（NEGF）を組み合わせた手法は、外部バイアス電位下での量子輸送をモデル化するための標準的なアプローチであるが、従来の計算手法は静的な系に限定されている。既存の実装は通常、固定された原子位置における輸送特性を計算するカノニカルな枠組み内で動作する。これにより、理論的理解とマクロな実験結果との間に大きな隔たりが生じている。特に、変換プロセスが本質的にグランドカノニカルである電気化学（EC）システムにおいて、この問題は顕著である。電気化学的な変換が起こる自然な環境を正確に記述するためには、バイアス下での分子動力学（MD）と、計算セルに入り込み、あるいは出ていく電子の明示的な開放境界の記述が必要となる。さらに、既存のDFT-NEGF実装の多くは、構造最適化やダイナミクスに必要な電流誘起力を計算する能力を欠いているか、あるいは凝縮相シミュレーションに最適化された最先端のパッケージが持つ機能を備えていない特定のDFTコード（例：SIESTA）に制限されている。

手法
著者らは、普及している混合ガウス/平面波電子構造パッケージであるCP2Kと、NEGFコードであるSMEAGOLとの間の新しいインターフェースを提示している。CP2K+SMEAGLEと名付けられたこのインターフェースは、元のSMEAGOLコードを可能な限り再利用しつつ、あらゆるDFTソフトウェアとリンク可能なスタンドアロンのライブラリとして設計されている。

主な技術的実装は以下の通りである：

• 電流誘起力： このインターフェースは、バイアス下でイオンに作用する非保存的な力の計算を実装している。これらの力は、原子位置に関するフル・ハミルトニアンの勾配を利用して、原子運動量の時間微分から導出される。極めて重要な点として、本実装には密度行列（$\Omega_{ij}$）の微分から生じる力項が含まれており、これは「電子の風（electron wind）」成分（電子からイオンへの運動量転移）を考慮したものである。この項は他の実装では無視されることが多いが、定性的な正当性を確保するために不可欠であることが示されている。
• 密度行列の構築： バイアス窓内の束縛状態を正しく占有するために、インターフェースは密度行列の計算を更新する。単一のコンター評価や左右の寄与の単純な平均を用いるのではなく、Brandbygeら（Brandbyge et al.）に従い、非平衡寄与によって決定される重みを用いた、重み付きダブルコンター・スキームを採用して積分を合計する。
• シミュレーション設定： 本手法は、バイアス下でのボーン・オッペンハイマー分子動力学を可能にする。系は、半無限電極に取り付けられた中央の「拡張分子」領域でモデル化される。リードの自己エネルギーは、中央領域の電荷密度とは独立して静的であるものとして扱われる。これは、中央領域が電荷の変化を遮蔽するのに十分な大きさである場合に有効な仮定である。

主な結果と検証
著者らは、いくつかの系に対してSIESTA+SMEAGOLを用いた参照計算と比較することで、CP2K+SMEAGOLインターフェースの検証を行った。

1. ゼロバイアス力（無限金ワイヤ）： インターフェースは、ゼロバイアスにおける原子力を正常に再現し、CP2K単独の計算と一致した。また、エネルギー密度行列項（$\Omega$）を無視すると、定性的に誤った力が算出されることも実証している。
2. 平行平板キャパシタ： インターフェースは、キャパシタにおける期待される電位降下を再現し、リード内では一定の電位を示し、真空中で線形に降下することを示した。また、高バイアスにおいて対称な電荷密度と正しいハートリーポテンシャルを生じさせるためには、重み付きダブルコンター・スキームが必要であることを確認した。一方、シングルコンター評価では、非対称で誤った結果が得られた。
3. Au-H$_2$-Au 接合の構造最適化： 印加バイアス下（最大1.5 V）において、インターフェースはH-H結合の伸長および原子変位を正しく予測し、これは従来のSIESTA+SMEAGOLの結果と一致している。これらの変化は、電場ではなく電流に起因するものである。
4. 溶媒和された金ワイヤの分子動力学： 著者らは、大規模な凝縮相系（166個の水分子によって溶媒和された単原子金ワイヤ）に対する初のDFT-NEGF分子動力学シミュレーションを実施した。
   • 安定性： 開放系であるにもかかわらず、ダイナミクスは安定しており、長期的なエネルギードリフトは最小限であった。
   • 電子移動： 1 Vのバイアス下での200 fsのトラジェトリにおいて、電子移動の効果が観察された。具体的には、局所的な電子密度の最小値が、中央のワイヤ原子から右端の原子へとシフトしており、これは右側の電極の正電荷を遮蔽する挙動と一致している。
   • 透過： 溶媒和されたワイヤの透過スペクトルは、理想的な無限ワイヤと定性的に類似していたが、その差異は水分子によって開かれた追加のトンネルチャネルに起因していた。

パフォーマンス
主要な計算上のボトルネックは、グリーン関数の積分としての電子密度の評価であり、1回のSCFステップあたり、単一のハミルトニアン対角化に代わって、96個のエネルギーポイント計算（実軸64点、複素軸32点）を必要とする。

• 並列化： コードは、実空間計算とk点のためにMPIを利用し、DOループとLAPACKルーチンのためにOpenMPを利用している。性能は32 OpenMPスレッドまで良好にスケールする。
• コスト： 溶媒和された金ワイヤ系において、CP2K+SMEAGOLの計算は、標準的なCP2Kの計算と比較して、1回のSCFステップあたり約3.2倍遅かった（CP2Kは256コアに対し、本手法は8倍の2048コアを使用）。著者らは、これはコストがオーダー単位で増加することを意味しており、さらなる加速化（例：ScaLAPACKや機械学習力場）なしでは、長時間のMDトラジェトリの実行が困難であることを指摘している。

意義と主張
本論文は、CP2K+SMEAGOLが、現在利用可能なCP2K向けの最も機能豊かで完全なDFT-NEGF実装であると主張している。これは、$\Gamma$点輸送（力を含まない）に限定されているネイティブのCP2K実装を凌駕するものである。

本研究の主要な意義は、DFT-NEGFを用いて、現実的な動作条件下での大規模な凝縮相系における分子動力学シミュレーションを実行できることを実証した点にある。著者らの知る限り、これはこのような計算が行われた初めての事例である。このインターフェースは、印加電位、電流の流れ、および原子ダイナミクスの相互作用が極めて重要となる電気触媒反応や電気化学システムの研究を可能にし、微視的な理論モデルとマクロな実験条件との間の溝を埋めるものである。