QM/MM 法を組み合わせた Red Moon 法による SEI 膜形成シミュレーション

Abstract

Red Moon (RM) 法 (混合モンテカルロ (MC)/分子動力学反応法) は複合化学反応系のための実用的な原子レベルのシミュレーションを実現する有用な手法である．RM法では，化学反応はMC法に基づき確率的に取り扱う． 本稿では，最近開発された，RM法における量子力学/分子力学法を用いた新しいエネルギー評価法及び，応用計算であるリチウムイオン電池の還元生成物同士における二量化反応の分岐に着目した固体電解液相間膜形成シミュレーションについて解説する．

1 はじめに

実験研究において化学反応を理解するための指針として原子レベルの計算機シミュレーションは非常に有効であるため，電解液や生物系のシステムにおける溶液中の化学反応過程を追跡するために，第一原理分子動力学 (AIMD) 法が広く用いられている [1,2,3]．しかし，この方法は計算コストが非常に高くシミュレーション時間，系のサイズが大きく制限される．より大きな系を取り扱うために，密度汎関数強束縛法 [4,5]を用いた分子動力学 (MD) 法や反応力場を用いたMD法 [6,7,8]が開発されている．これらの方法はAIMD計算より，大きく計算コストが削減できるが，多数の化学反応からなる複合化学反応系を取り扱うことは難しい．そこで，大規模な複合化学反応系を取り扱うため，Red Moon (RM) 法 (混合モンテカルロ (MC)/ MD反応法) を開発してきた [9,10]．これまで，N，N-ジメチルホルムアミド溶媒中の2-クロロブタンのラセミ化反応 [9]，芳香族ポリアミド膜の形成 [11]，二次電池の固体電解液相間 (SEI) 膜形成 [12,13,14,15,16]，オレフィン重合 [17]に適用されている．SEI膜は二次電池の初回充電サイクル時に負極上に電解液の還元生成物として生じる．SEI膜は後続の充放電サイクルにおいて電解液の更なる還元分解を防ぐと共に，電極-電解液間でSEI膜を通したLiイオンの交換が可能であるため，二次電池の安定性に大きく寄与する [18,19]．RM法によるシミュレーションでは，これまでの実験研究では得ることが難しかった原子レベルの構造を理論的に提供することができた．更に，二次電池におけるSEI膜の溶媒分子 [12]や添加剤 [13,16]，塩濃度 [14]への依存性に関する微視的機構の理解に有用であることを示した．

RM法は複合化学反応過程をMC法とMD法を組み合わせたサイクル (MC/MDサイクル) を繰り返すことによりシミュレートできる．RM法においては，MC法で化学反応などを確率的に取り扱い，MD法で短時間スケールの分子運動を取り扱う．前者の処理において，それぞれの化学反応の遷移確率は反応前後のポテンシャルエネルギー差を用いて評価される．しかし，通常，RM法におけるMC処理のポテンシャルエネルギーは分子力学 (MM) 法により評価するが，この評価方法にはまだ精度向上の余地がある．精度向上を行う方法として，系全体を量子力学的に取り扱うことが考えられるが，計算コストが非常に高く現実的ではない．そこで，計算コストを削減する現実的な取り扱い手法として，反応に関わる箇所を量子力学的に，それ以外の箇所を分子力学的に取り扱う量子力学 (QM)/MM法 [20,21,22]を適用する．RM法におけるMCスキームへのQM/MM法の導入は計算コストと計算精度のバランスが良いと考えられる．

本稿では，最近開発されたRM法を拡張した手法であるQM/MM-RM法およびその応用研究について解説する．

2 計算手法

2.1 Red Moon法の概要

Red Moon (RM) 法の名称は，Rare Event-Driving MethOdology Of Necessityの頭字語に由来する [9,10]．この方法では，短時間スケールの分子運動をMD法，長時間スケールで起こる化学反応をMC法により取り扱う．このMD法とMC法による一連の処理をMC/MDサイクルと呼び，MC/MDサイクルを繰り返すことにより，従来の手法では再現できなかった長時間スケールの複合化学反応を取り扱うことができる [9,10]．

Figure 1にMC/MDサイクル処理の概要を示す．初めに，配置状態 rにおいて反応候補を探索する (Figure 1 (a))．反応候補数$N_{\text{cand}}$ が0より大きい場合，試行する化学反応として，活性化障壁$E_{\text{a}}$を用いた指数因子$R_{\text{a}}\sim \exp (-\beta E_{\text{a}})$ に比例して反応候補を1個選択する (Figure 1 (b)，(c))(反応候補の選択方法については「3．計算条件」参照)．このとき，反応物の力場パラメータとポテンシャル関数形を生成物のものへと切り替え，起こり得る配置sを生成し，古典MDシミュレーションにより系全体を緩和する．エネルギー変化$\Delta U_{rs}$と遷移確率$W_{r\to s}$ (2.2節を参照) を計算するため，MC法に従って試行した反応の採択，棄却を判定する (Figure 1 (e))．反応候補が存在しない場合，これらの処理の代わりに短時間のMDシミュレーションを実行する (Figure (g))．最後に，もし分子混合物の組成がほとんど変化しない場合は計算を終了する．それ以外の場合はMC/MDサイクルを繰り返す (Figure 1 (h))．

Figure 1.

 Schematic flowchart of a computational scheme for the MC/MD cycle in the RM method. Reproduced from J. Chem. Phys., 149, 044113 (2018), with permission of AIP Publishing.

2.2 RM法のMCスキーム

RM法のMCスキームでは，メトロポリス法による遷移確率,   

$W_{r\to s}=\min \{1,\exp (-\beta \Delta U_{rs}\left)\right\}$(1)

は状態rとsのポテンシャルエネルギー変化$\Delta U_{rs}(=U_s-U_r)$ によって得られる．

従来のRM法では，$\Delta U_{rs}$は次式で得られる．   

$\Delta U_{rs}=\Delta U_{rs}^{\text{MM}}+\Delta U_0^{\text{reac}}$(2)

ここで，$\Delta U_{rs}^{\text{MM}}$は状態rとsのポテンシャルエネルギー差，$\Delta U_0^{\text{reac}}$は反応に伴うポテンシャルエネルギー変化である．以降，このMM法による従来のRM法をMM-RM法と呼ぶ.

一方で，QM/MM法によるエネルギー評価をRM法中のMCスキームに入れることは自然な拡張である．この拡張した手法では$\Delta U_{rs}$を次式で評価する．   

$\Delta U_{rs}=\Delta U_{rs}^{\text{QM/MM}}$(3)

ここで，$\Delta U_{rs}^{\text{QM/MM}}$は状態rとsのポテンシャルエネルギー差であり，QM/MM法によって直接得られる．今回の目的のため，反応物及び生成物である溶質分子を量子力学的に取り扱い，それ以外の分子を分子力学的に取り扱う．このようなQM/MMエネルギー評価を用いたRM法を，以降はQM/MM-RM法と呼ぶ．MM-RM法とQM/MM-RM法の違いはFigure 1 (d)における$\Delta U_{rs}$の評価方法のみである．

本研究では，$\Delta U_{rs}$ (式(1)参照) を次式で評価する.   

$\Delta U_{rs}=\Delta U_{rs}^{\text{sol}}+\Delta U_{rs}^{\text{sol-svt}}$(4)

ここで，$\Delta U_{rs}^{\text{sol}}$は溶質分子のポテンシャルエネルギー差，$\Delta U_{rs}^{\text{sol-slv}}$は溶質-溶媒間相互作用エネルギー差である．今回のRM法では，溶質分子は反応物または生成物であり，溶媒分子はそれ以外の分子を指す．計算コストを大幅に削減するため，溶媒-溶媒間相互作用を無視する．なぜなら，この値は式(4)の他の2個の項と比較して無視できるからである．よって， $\Delta U_{rs}^{\text{sol-slv}}$は次式のように$\Delta U_{rs}^{\text{IN}}$と$\Delta U_{rs}^{\text{OUT}}$に分割することができる.   

$\Delta U_{rs}^{\text{sol-svt}}=\Delta U_{rs}^{\text{IN}}+\Delta U_{rs}^{\text{OUT}}$(5)

ここで，$\Delta U_{rs}^{\text{IN}}$は溶質とカットオフ距離内の溶媒との相互作用エネルギー差であり，$\Delta U_{rs}^{\text{OUT}}$は溶質とカットオフ距離外の溶媒との相互作用エネルギー差である.

2.3 QM/MM-RM法におけるポテンシャルエネルギー評価方法

QM/MM法のハミルトニアンは次式で定義される．   

${\hat{H}}_{\text{QM/MM}}={\hat{H}}_{\text{QM}}+{\hat{H}}_{\text{QM-MM}}+{\hat{H}}_{\text{MM}}$(6)

${\hat{H}}_{\text{QM}}$ と${\hat{H}}_{\text{MM}}$はそれぞれQM, MM領域のハミルトニアンである．一方で，${\hat{H}}_{\text{QM-MM}}$はQM-MM領域間の相互作用であるカップリング項である．周期境界条件下のQM/MM法で，系を3個の領域に分割した．すなわち，(i) 溶質領域であるQM領域 (Figure 2青色領域)，(ii) カットオフ距離内の溶媒領域 (IN領域 (Figure 2破線内側の白色領域))，(iii) カットオフ距離外の溶媒領域 (OUT領域 (Figure 2破線外側の白色領域)) である．よって，QM/MM法による総エネルギーは次式のように記述される.   

$U_{\text{QM/MM}}=⟨\Psi |{\hat{H}}_{\text{QM/MM}}|\Psi ⟩$

(7a)   

$=⟨\Psi |{\hat{H}}_{\text{QM}}+{\hat{H}}_{\text{QM-MM}}|\Psi ⟩+U^{\text{vdW}}+U^{\text{MM-sub}}$

, (7b)   

$=U^{\text{QM-sub}}+U^{\text{es-IN}}+U^{\text{es-OUT}}+U^{\text{vdW}}+U^{\text{MM-sub}}$

, (7c)

Figure 2.

 Schematic illustration of potential energy terms. The central blue region shows the QM region and that surrounded by a blue dotted closed curve shows the cutoff region, while the black solid lines show periodic boundaries. Reproduced from J. Chem. Phys., 149, 044113 (2018), with permission of AIP Publishing.

ここで，$\Psi$ はQM領域を記述する波動関数,$U^{\text{QM-sub}}$ はQM領域のポテンシャルエネルギー，$U^{\text{es-IN}}$はIN領域におけるQM原子とMM原子のQM-MM静電相互作用エネルギー，$U^{\text{es-OUT}}$はOUT領域におけるQM原子とMM原子のQM-MM静電相互作用エネルギー，$U^{\text{MM-sub}}$ はMM領域のポテンシャルエネルギーである．$U^{\text{vdW}}$はQM領域，QM領域における原子間のレナード-ジョーンズ (LJ) ポテンシャル関数で記述される．Figure 2にQM/MM法における各ポテンシャルエネルギー項の略図を示す (式(7c)参照)．MM-RM法と同様に，状態rとsのMM相互作用エネルギー差は無視する．また，今回はカットオフ法を用いるため，OUT領域とQM原子のQM-MM静電相互作用$U^{\text{es-OUT}}$も無視される．結果として，QM/MM-RM法では，式(4), (5)の各項は次式で得られる.   

$\Delta U_{rs}^{\text{sol}}=\Delta U_{rs}^{\text{QM-sub}}$(8)

  

$\Delta U_{rs}^{\text{IN}}=\Delta U_{rs}^{\text{es-IN}}+\Delta U_{rs}^{\text{vdW}}$(9)

  

$\Delta U_{rs}^{\text{OUT}}=0$(10)

最後に，QM/MM-RM法では，状態rと状態sのポテンシャルエネルギー変化は次式で計算される.   

$\Delta U_{rs}=\Delta U_{rs}^{\text{QM-sub}}+\Delta U_{rs}^{\text{es-IN}}+\Delta U_{rs}^{\text{vdW}}+\Delta U_{rs}^{\text{es-OUT}}$(11)

3 計算条件

電解液，反応生成物の力場パラメータとして，generalized AMBER force filed (GAFF) [23] を用いる．PF${}_6^{-}$の力場パラメータは過去に行われた理論研究を参照した [24]．すべての量子力学計算はGaussian09 [25]でB3LYP/6-31+G (d)の理論レベルを用いて実行した [26]．原子点電荷はrestrained electrostatic potential (RESP) 法を用いて得た [27]．このとき，実験的な双極子モーメントの再現の為，電解液の溶媒であるエチレンカーボネート (EC) 分子の電荷を0.8倍にスケーリングした [18]．QM/MMエネルギー評価のために，Amber9 [28]とGaussian09 [25]を利用したAmber-Gaussianインターフェースプログラム [29]を用いた．LJ相互作用のカットオフ距離は8 Åとした．QM/MM法では，QM原子とMM原子間のカットオフ距離は12 Åとした．外部のMM点電荷を伴うQM電荷密度の静電相互作用は，QM領域のハミルトニアンにMM点電荷の寄与を加えて計算するelectrostatic embedding法 [30]を用いて評価した．

すべてのMD計算にはAmber9のSANDERモジュール [28]を用いた．SHAKE法 [31]により水素原子と他の重い原子の結合を拘束し，積分の時間刻みを1 fsに設定した．温度はBerendsen thermostat [32]により298 K に維持した．MC/MDサイクルにおける緩和のMDシミュレーション時間は4 psに設定した．本研究では，分子間における原子対の距離がファンデルワールス (vdW) 半径の和よりも小さいかどうかを条件として，反応候補の探索を行った．ここで，$R_{\text{a}}$の値はすべての化学反応で1.0とした．分子の可視化ソフトウェアとして，Visual Molecular Dynamics (VMD) [33]を用いた．

Figure 3に本シミュレーションにおけるモデルと反応スキームを示す．過去の実験研究によれば，SEI層は主にエチレンジカーボネート (EDC^2-) で構成されていると一般的に信じられている [34,35,36]．一方で，いくつかの理論研究ではブチレンジカーボネート (BDC^2-) がより安定な生成物であると言われている [37,38,39]．このような実験と理論研究で意見の分かれる矛盾の起源はまだ明らかにされていない．それゆえ，本研究では，${\text{2EC}}^{·-}$の二量化反応 ($2{\text{EC}}^{·-}\to {\text{EDC}}^{\text{2-}}+{\text{C}}_{\text{2}}{\text{H}}_{\text{4}}$または$2{\text{EC}}^{·-}\to {\text{BDC}}^{\text{2-}}$)の分岐に注目する．本モデルは1.1 mol L⁻¹ LiPF₆/EC電解液と負に帯電 (−2 e) した負極を含んでいる [12]．それに応じて，二個のLiカチオンを電解液に追加することで系を中性化した．電解液はグラファイト(層間距離3.7Å)様の負極と真空領域に隣接している．負極構造はシミュレーション中で固定し，電解液分子が電極中に侵入しないように制限した．初期構造依存性を排除するため，10個の独立した初期構造を用いてシミュレーションを実行した.

Figure 3.

 Model simulation system (3.41 nm × 3.70 nm × 30.0 nm) and the reaction scheme (cyan, carbon; red, oxygen; white, hydrogen; orange, phos- phorus; green, fluorine; pink, Li⁺). The origin of the z coordinate is taken at the position of the carbon atoms on the anode surface in contact with the electrolyte. The right side numbers represent the number of ions/molecules of the simulated system.

QM/MM-RM法との比較のため，MM-RM法でも同様の計算を行った．本計算において，$\Delta U_0^{\text{reac}}$(式(2)参照) は量子力学計算により得た (Table 1).

Table 1. Reaction energy $\Delta U_0^{\text{reac}}$ in isolated system at the B3LYP/6-31+G (d) level of theory.

Elementary reaction	$\Delta U_0^{\text{reac}}$ [kcal·mol−1]
${\text{EC + e}}^{\text{-}}\to {\text{EC}}^{·\text{-}}$	-5.2
${\text{2EC}}^{·-}\to {\text{BDC}}^{\text{2-}}$	-88.5
${\text{2EC}}^{·-}\to {\text{EDC}}^{\text{2-}}{\text{ + C}}_2{\text{H}}_4$	-16.3
${\text{EC}}^{·-}{\text{ + e}}^{\text{-}}\to {\text{CO}}_3^{2-}$	-44.2
${\text{PF}}_6^{-}{\text{ + e}}^{\text{-}}\to {\text{PF}}_6^{2-}$	116.9
${\text{PF}}_6^{2-}{\text{ + e}}^{\text{-}}\to {\text{PF}}_3+3{\text{F}}^{-}$	-55.1

4 結果と考察

今回実行したRMシミュレーションでは，${\text{EC}}^{·-}$, BDC^2-, EDC^2-, CO${}_3^{2-}$がSEI膜を形成し，この膜はEC分子が負極表面に接近することを防ぐ．Figure 4は平衡状態 (∼2,500 MC/MDサイクル時) における負極表面付近のSEI膜，ガス分子，電解液の代表的なスナップショットを示す．ここでz=0は負極表面を表している．このときのSEI膜の時間平均質量密度分布${\rho }_m$ (Figure 5) は平衡状態 (2,500 MC/MDサイクル時) を初期構造として追加のMDシミュレーションを実行して得た． 結果として，EDC^2- (紫色) の量がBDC^2- (緑色) の量に比べて多くなることがわかった．今回のシミュレーションで得られた傾向は，LIBのEDC^2-がSEI膜の主要な構成物であるという実験的観測とよく似ている [34,35,36]．一方で，従来のMM法を用いたRM法 (MM-RM法) による追加のSEI膜形成シミュレーションでは，EDC^2-の反応がほとんど棄却され，代わりに相対的にBDC^2-が採択されたためBDC^2-の量が大幅に増えた (Figure 6)．これは，MM-RM法のMCスキームにおけるEDC^2-への二量化に関する$\Delta U_{rs}$の不正確な推定に起因しており，それゆえ，今回のQM/MM-RM法はRM法のMCスキームの精度向上により，SEI膜の実験の特徴を適切に再現することができる．実際，平均値$<\Delta U_{rs}>$の値は，MM-RM法とQM./MM-RM法で，それぞれ126.8 ± 2.5 kcal mol⁻¹と−52.9 ± 3.1 kcal mol⁻¹となる.

Figure 4.

 A typical snapshot of the SEI films and LiPF₆/EC-based electrolytes in the QM/MM-RM method (blue, Li⁺ in the SEI film). In this panel, the solvent EC and gas molecules are depicted by a stick model, while the remaining molecules are depicted by a ball model (full vdW radii).

Figure 5.

 Mass density distribution ${\rho }_m$ of the SEI film components in the atomistic reaction simulation using the QM/MM-RM methods. Reproduced from J. Chem. Phys., 149, 044113 (2018), with permission of AIP Publishing.

Figure 6.

 Mass density distribution ${\rho }_m$ of the SEI film components in the atomistic reaction simulation using the MM-RM methods.

このとき，Figure 7はRMシミュレーション中のSEI膜における反応生成物の表面数密度${\rho }_{\text{n}}^{\text{s}}$ 変化を示す．EDC^2-とBDC^2-の${\rho }_{\text{n}}^{\text{s}}$を比較すると，収束時におけるEDC^2-の値 (6.0 nm^-2) はBDC^2-の値 (3.0 nm^-2) より大きくSEI膜中で豊富であることを示している．さらに，EDC^2- (紫色) への二量化はBDC^2- (緑色) に比べてすべてのSEI膜形成プロセス中で優先的に発生している．しかし，BDC^2-への二量化における${\rho }_{\text{n}}^{\text{s}}$曲線の勾配が600 MC/MDサイクルの後に大きくなっている．このMC/MDサイクル数は${\text{EC}}^{·-}$の表面数密度のピークに一致していており，${\text{EC}}^{·-}$の濃度が高い状態においてBDC^2-への二量化反応が促進されていることを示唆している.

Figure 7.

 10 Change in the surface number densities ${\rho }_{\text{n}}^{\text{s}}$ of the reaction products during the SEI film formation processes with respect to the MC/MD cycles in the atomistic reaction simulations. Reproduced from J. Chem. Phys., 149, 044113 (2018), with permission of AIP Publishing.

EDC^2-への優先的な二量化の微視的機構を明らかにするため，RMシミュレーション中の化学反応条件を満たした二量化反応の候補対の表面数密度${\rho }_{\text{cand}}$ (Figure 8) を調べた．結果として，EDC^2- (紫色) への二量化の${\rho }_{\text{cand}}$は，BDC^2- (緑色) への二量化反応ものより大きかった．この結果は，多くの${\text{EC}}^{·-}$はEDC^2-への二量化の反応条件 (O_aとC_bとの接触) を満たすように配向することを明確に意味する．過去の理論研究では，BDC^2-への二量化がEDC^2-の反応より発熱的であると報告されている [37]．反応間の熱力学的な違いだけではなく，QM/MM-RM法は瞬間的な配向や拡散のような動力学的効果を考慮することで，より現実的に二量化反応を取り扱うことができたと結論付けられる．

Figure 8.

 Change in the surface number densities ${\rho }_{\text{cand}}$ of reaction candidate pairs of dimerization reactions in the atomistic reaction simulations. Atom names of the EC radical anion are shown with the molecule drawn in the Figure. Reproduced from J. Chem. Phys., 149, 044113 (2018), with permission of AIP Publishing.

5 おわりに

本稿では，RM法におけるMCスキームの遷移確率をより正確に評価して適用対象を増やすことを目的とした，QM/MM法を利用したRM法，QM/MM-RM法を解説した．実用的な応用への例として，本手法をSEI膜の主要な構成物として実験で観測された二量体の違いに注目して，EC電解液を用いたLIBのSEI膜形成の解析へ適用した．実際にシミュレーションによって得られたSEI膜は主要な構成物が還元生成物の二量体であるEDC^2-であり，実験研究と同じ傾向を再現することができた．

MM-RM法においても，力場パラメータやポテンシャル関数形を修正することで計算精度の向上は可能である．しかし，特定の対象系の実験的特性を再現するためには力場パラメータを慎重に最適化する必要があり，通常このような変更は高コストである．一方，QM/MM-RM法では，瞬間的な溶質や溶媒の構成に応じて量子力学的にポテンシャルエネルギーを得ることができるため，力場選択の依存性を大きく減らすことができる．また，今回のQM/MM-RM法では，反応候補探索にはMM-MD法を使っているため，QM/MM法の適用は，その候補に対する反応前後の反応エネルギーを計算する(2回の一点計算)のみであるため，計算コストは少なく，本系におけるQM/MM-RM法の計算時間は，MM-RM法の高々1.5倍程度であった．このような結果から，QM/MM法によるエネルギー評価を用いたRM法，QM/MM-RM法は様々な複合化学反応が関与する材料設計及び機能開発に重要な貢献をすると期待している．

今回のシミュレーションでは，負極にグラファイト様の炭素電極モデルを用いた．一方，実験研究において，炭素電極の焼成温度を変えることで大きく特性が異なることが知られており [40]，この傾向は電極構造に起因していると考えられる．そこで，現在，反応力場を用いたMDシミュレーションにより不均一な表面構造を持つ炭素電極モデルを作成し，電極-電解液界面の解析を進めている (Figure 9)．このような不均一な表面構造を持つ電極と電解液の界面におけるSEI膜形成を取り扱うには，十分大きなサイズの電極表面を用いた大規模なRMシミュレーションを実施する必要があり，スーパーコンピュータのような並列計算機の活用が不可欠である.

Figure 9.

 Electrode-electrolyte interface model using complex carbon electrode.

謝辞

本研究は， 文部科学省・元素戦略研究拠点「実験と理論計算科学のインタープレイによる触媒・電池の元素戦略研究拠点 (ESICB) 」，文部科学省ポスト「京」重点課題5「エネルギーの高効率な創出，変換・貯蔵，利用の新規基盤技術の開発」，科学技術振興機構，CREST「マクロ化学現象シミュレーションに向けた計算分子技術の構築」の支援を受けて実施された．また，計算の一部は自然科学研究機構 岡崎共通研究施設 計算科学研究センター及び名古屋大学情報基盤センターの計算機を利用して行った．

参考文献

• [1]   R. Car, M. Parrinello, Phys. Rev. Lett., 55, 2471 (1985). doi:10.1103/PhysRevLett.55.247110032153
• [2]   M. E. Tuckerman, J. Phys. Condens. Matter, 14, R1297 (2002). doi:10.1088/0953-8984/14/50/202
• [3]   M. C. Payne, M. P. Teter, D. C. Allan, T. A. Arias, J. D. Joannopoulos, Rev. Mod. Phys., 64, 1045 (1992). doi:10.1103/RevModPhys.64.1045
• [4]   D. Porezag, T. Frauenheim, T. Köhler, G. Seifert, R. Kaschner, Phys. Rev. B, 51, 12947 (1995). doi:10.1103/PhysRevB.51.12947
• [5]   M. Elstner, D. Porezag, G. Jungnickel, J. Elsner, M. Haugk, T. Frauenheim, Matter Mater. Phys., 58, 7260 (1998). doi:10.1103/PhysRevB.58.7260
• [6]   J. Tersoff, Phys. Rev. Lett., 56, 632 (1986). doi:10.1103/PhysRevLett.56.63210033244
• [7]   D. W. Brenner, Phys. Rev. B, 42, 9458 (1990). doi:10.1103/PhysRevB.42.9458
• [8]   A. C. T. van Duin, S. Dasgupta, F. Lorant, W. A. Goddard, J. Phys. Chem. A, 105, 9396 (2001). doi:10.1021/jp004368u
• [9]   M. Nagaoka, Y. Suzuki, T. Okamoto, N. Takenaka, Chem. Phys. Lett., 583, 80 (2013). doi:10.1016/j.cplett.2013.08.017
• [10]   Y. Suzuki, M. Nagaoka, J. Chem. Phys., 146, 204102 (2017). doi:10.1063/1.498339628571340
• [11]   Y. Suzuki, Y. Koyano, M. Nagaoka, J. Phys. Chem. B, 119, 6776 (2015). doi:10.1021/jp512333h25973839
• [12]   N. Takenaka, Y. Suzuki, H. Sakai, M. Nagaoka, J. Phys. Chem. C, 118, 10874 (2014). doi:10.1021/jp5018696
• [13]   N. Takenaka, H. Sakai, Y. Suzuki, P. Uppula, M. Nagaoka, J. Phys. Chem. C, 119, 18046 (2015). doi:10.1021/acs.jpcc.5b04206
• [14]   N. Takenaka, T. Fujie, A. Bouibes, Y. Yamada, A. Yamada, M. Nagaoka, J. Phys. Chem. C, 122, 2564 (2018). doi:10.1021/acs.jpcc.7b11650
• [15]   T. Fujie, N. Takenaka, Y. Suzuki, M. Nagaoka, J. Chem. Phys., 149, 044113 (2018). doi:10.1063/1.503477130068180
• [16]   A. Bouibes, N. Takenaka, T. Fujie, K. Kubota, S. Komaba, M. Nagaoka, ACS Appl. Mater. Interfaces, 10, 28525 (2018). doi:10.1021/acsami.8b0753030070476
• [17]   K. Matsumoto, M. Takayanagi, Y. Suzuki, N. Koga, M. Nagaoka, J. Comput. Chem., 40, 421 (2019). doi:10.1002/jcc.2570730351517
• [18]   K. Xu, Chem. Rev., 104, 4303 (2004). doi:10.1021/cr030203g15669157
• [19]   K. Xu, Chem. Rev., 114, 11503 (2014). doi:10.1021/cr500003w25351820
• [20]   A. Warshel, M. Levitt, J. Mol. Biol., 103, 227 (1976). doi:10.1016/0022-2836(76)90311-9985660
• [21]   U. C. Singh, P. A. Kollman, J. Comput. Chem., 7, 718 (1986). doi:10.1002/jcc.540070604
• [22]   M. J. Field, P. A. Bash, M. Karplus, J. Comput. Chem., 11, 700 (1990). doi:10.1002/jcc.540110605
• [23]   J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman, D. A. Case, J. Comput. Chem., 25, 1157 (2004). doi:10.1002/jcc.2003515116359
• [24]   Z. Liu, S. Huang, W. Wang, J. Phys. Chem. B, 108, 12978 (2004). doi:10.1021/jp048369o
• [25]   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, M.C. X. Li, A. Marenich, J. Bloino, B.G. Janesko, R. Gomperts, B. Mennucci, H.P. Hratchian, J. V. Ortiz, A.F. Izmaylov, et al., Gaussian 09, Revision A.02 (2009).
• [26]   J. M. Vollmer, L. A. Curtiss, D. R. Vissers, K. Amine, J. Electrochem. Soc., 151, A178 (2004). doi:10.1149/1.1633765
• [27]   C. I. Bayly, P. Cieplak, W. Cornell, P. A. Kollman, J. Phys. Chem., 97, 10269 (1993). doi:10.1021/j100142a004
• [28]   D.A. Case, T.A. Darden, T.E. Cheatham, C.L. Simmerling, J. Wang, R.E. Duke, R. Luo, J. Merz, K. M., D.A. Pearlman, M. Crowley, R.C. Walker, W. Zhang, B. Wang, S. Hayik, A.E. Roitberg, G. Seabra, K.F. Wong, F. Paesani, X. Wu, et al., Amber9 (2006).
• [29]   T. Okamoto, K. Yamada, Y. Koyano, T. Asada, N. Koga, M. Nagaoka, J. Comput. Chem., 32, 932 (2011). doi:10.1002/jcc.2167820949515
• [30]   D. Bakowies, W. Thiel, J. Phys. Chem., 100, 10580 (1996). doi:10.1021/jp9536514
• [31]   J. P. Ryckaert, G. Ciccotti, H. J. C. Berendsen, J. Comput. Phys., 23, 327 (1977). doi:10.1016/0021-9991(77)90098-5
• [32]   H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, J. Chem. Phys., 81, 3684 (1984). doi:10.1063/1.448118
• [33]   W. Humphrey, A. Dalke, K. Schulten, J. Mol. Graph., 14, 33, 27 (1996). doi:10.1016/0263-7855(96)00018-58744570
• [34]   D. Aurbach, Y. Ein‐Ely, A. Zaban, J. Electrochem. Soc., 141, L1 (1994). doi:10.1149/1.2054718
• [35]   D. Aurbach, Y. Ein-Eli, O. Chusid, J. Electrochem. Soc., 141, 603 (1994). doi:10.1149/1.2054777
• [36]   M. Nie, D. Chalasani, D. P. Abraham, Y. Chen, A. Bose, B. L. Lucht, J. Phys. Chem. C, 117, 1257 (2013). doi:10.1021/jp3118055
• [37]   Y. Wang, S. Nakamura, M. Ue, P. B. Balbuena, J. Am. Chem. Soc., 123, 11708 (2001). doi:10.1021/ja016452911716728
• [38]   D. Bedrov, G. D. Smith, A. C. T. van Duin, J. Phys. Chem. A, 116, 2978 (2012). doi:10.1021/jp210345b22352952
• [39]   M. M. Islam, G. Kolesov, T. Verstraelen, E. Kaxiras, A. C. T. van Duin, J. Chem. Theory Comput., 12, 3463 (2016). doi:10.1021/acs.jctc.6b0043227399177
• [40]   J.R. Dahn, T. Zheng, Y. Liu, and J.S. Xue, Science (80-.)., 270, 590–593 (1995). doi:10.1126/science.590