エバルトの方法

エバルトの方法（エバルトのほうほう、英: Ewald method）は、分子動力学法、量子化学的手法、バンド計算手法などで、単位胞（またはスーパーセル）内の原子核（またはイオン芯）同士のクーロン相互作用を効率良く計算する手法である。

実空間と逆格子空間のどちらでも発散してしまうクーロン相互作用を、実空間での計算が都合のよい部分と、逆格子空間での計算が都合のよい部分との2つに分けて別々に計算を行い、これら2つの計算結果の和が求めるべき答となる。

具体例

格子点 l {\displaystyle {\boldsymbol {l}}} におかれたイオンの作る静電ポテンシャル ϕ ( r ) {\displaystyle \phi ({\boldsymbol {r}})} を例にとって説明する。 ϕ ( r ) {\displaystyle \phi ({\boldsymbol {r}})} は係数を省略すれば次の式で書ける。

ϕ ( r ) = ∑ l 1 | r − l | {\displaystyle \phi ({\boldsymbol {r}})=\sum _{\boldsymbol {l}}{\frac {1}{\left|{\boldsymbol {r}}-{\boldsymbol {l}}\right|}}}

静電ポテンシャル ϕ ( r ) {\displaystyle \phi ({\boldsymbol {r}})} を、関数 f ( r ) {\displaystyle f({\boldsymbol {r}})} を用いて、次のように分割する。

ϕ ( r ) = ∑ l f ( r ) | r − l | + ∑ l 1 − f ( r ) | r − l | {\displaystyle \phi ({\boldsymbol {r}})=\sum _{\boldsymbol {l}}{\frac {f({\boldsymbol {r}})}{\left|{\boldsymbol {r}}-{\boldsymbol {l}}\right|}}+\sum _{\boldsymbol {l}}{\frac {1-f({\boldsymbol {r}})}{\left|{\boldsymbol {r}}-{\boldsymbol {l}}\right|}}}

ここで、 f ( r ) {\displaystyle f({\boldsymbol {r}})} として短距離で0に収束する関数をうまく選び、第1項の和が実空間で短距離で0に収束し、第２項の和が逆格子空間で短距離で0に収束すると計算上都合がよい。 実際には、 f ( r ) {\displaystyle f({\boldsymbol {r}})} として相補誤差関数 erfc がよく使われる。Gを任意の定数として f ( r ) = erfc ⁡ ( G | r − l | ) {\displaystyle f({\boldsymbol {r}})=\operatorname {erfc} (G\left|{\boldsymbol {r}}-{\boldsymbol {l}}\right|)} を、 1 − f ( r ) {\displaystyle 1-f({\boldsymbol {r}})} は誤差関数の定義式を代入すると

ϕ ( r ) = ∑ l 1 | r − l | erfc ⁡ ( G | r − l | ) + ∑ l 1 | r − l | 2 π ∫ 0 G | r − l | e − t 2 d t {\displaystyle \phi ({\boldsymbol {r}})=\sum _{\boldsymbol {l}}{\frac {1}{\left|{\boldsymbol {r}}-{\boldsymbol {l}}\right|}}\operatorname {erfc} (G\left|{\boldsymbol {r}}-{\boldsymbol {l}}\right|)+\sum _{\boldsymbol {l}}{\frac {1}{\left|{\boldsymbol {r}}-{\boldsymbol {l}}\right|}}{\frac {2}{\sqrt {\pi }}}\int _{0}^{G\left|{\boldsymbol {r}}-{\boldsymbol {l}}\right|}e^{-t^{2}}\,\mathrm {d} t}

第2項で t = | r − l | ρ {\displaystyle t=\left|{\boldsymbol {r}}-{\boldsymbol {l}}\right|\rho } とおいて変数変換する。

ϕ ( r ) = ∑ l 1 | r − l | erfc ⁡ ( G | r − l | ) + ∫ 0 G ∑ l 2 π exp ⁡ ( − | r − l | 2 ρ 2 ) d ρ {\displaystyle \phi ({\boldsymbol {r}})=\sum _{\boldsymbol {l}}{\frac {1}{\left|{\boldsymbol {r}}-{\boldsymbol {l}}\right|}}\operatorname {erfc} (G\left|{\boldsymbol {r}}-{\boldsymbol {l}}\right|)+\int _{0}^{G}\sum _{\boldsymbol {l}}{\frac {2}{\sqrt {\pi }}}\exp(-|{\boldsymbol {r}}-{\boldsymbol {l}}|^{2}\rho ^{2})\,\mathrm {d} \rho }

第２項の被積分関数は、格子の周期を持つ r {\displaystyle {\boldsymbol {r}}} の関数であるから、フーリエ級数に展開できる。単位格子の体積を v c e l l {\displaystyle v_{\mathrm {cell} }} とおくと、次のように第２項を展開できる。

ϕ ( r ) = ∑ l 1 | r − l | erfc ⁡ ( G | r − l | ) + ∫ 0 G ∑ g 2 π v c e l l 1 ρ 3 exp ⁡ ( − | g | 2 / 4 ρ 2 ) exp ⁡ ( i g ⋅ r ) d ρ = ∑ l 1 | r − l | erfc ⁡ ( G | r − l | ) + 4 π v c e l l ∑ g exp ⁡ ( − | g | 2 / 4 G 2 ) | g | 2 exp ⁡ ( i g ⋅ r ) {\displaystyle {\begin{aligned}\phi ({\boldsymbol {r}})&=\sum _{\boldsymbol {l}}{\frac {1}{\left|{\boldsymbol {r}}-{\boldsymbol {l}}\right|}}\operatorname {erfc} (G\left|{\boldsymbol {r}}-{\boldsymbol {l}}\right|)+\int _{0}^{G}\sum _{\boldsymbol {g}}{\frac {2\pi }{v_{\mathrm {cell} }}}{\frac {1}{\rho ^{3}}}\exp(-|{\boldsymbol {g}}|^{2}/4\rho ^{2})\exp(i{\boldsymbol {g}}\cdot {\boldsymbol {r}})\,\mathrm {d} \rho \\&=\sum _{\boldsymbol {l}}{\frac {1}{\left|{\boldsymbol {r}}-{\boldsymbol {l}}\right|}}\operatorname {erfc} (G\left|{\boldsymbol {r}}-{\boldsymbol {l}}\right|)+{\frac {4\pi }{v_{\mathrm {cell} }}}\sum _{\boldsymbol {g}}{\frac {\exp(-\left|{\boldsymbol {g}}\right|^{2}/4G^{2})}{|{\boldsymbol {g}}|^{2}}}\exp(i{\boldsymbol {g}}\cdot {\boldsymbol {r}})\end{aligned}}}

よって、 erfc ⁡ ( x ) {\displaystyle \operatorname {erfc} (x)} や e − x {\displaystyle e^{-x}} といった、 x {\displaystyle x} について速やかに0に収束する関数が現れる形に変形することができた。 後は、各項がそれぞれ l {\displaystyle {\boldsymbol {l}}} と g {\displaystyle {\boldsymbol {g}}} に対して速く収束するように適当なGの値を選べば効率よく計算できる。

粒子・メッシュ・エバルト (PME) 法

エバルト法は、計算機の出現よりずっと以前に、理論物理学における手法として開発された。しかしながら、エバルト法は1970年代以降、粒子系のコンピュータシミュレーション、特に重力や静電気学といった逆2乗力を介して相互作用する粒子系において広範に使用されている。最近、PME法は打ち切りによるアーティファクトを除去するためにレナード-ジョーンズ・ポテンシャルの r − 6 {\displaystyle r^{-6}} 部分の計算にも使用されている[1]。PME法はプラズマ、銀河、分子のシミュレーションに応用されている。

粒子・メッシュ法では、標準のエバルト和と同じく、包括的相互作用ポテンシャルが2つの項へと分離される。

φ ( r )   = d e f   φ s r ( r ) + φ l r ( r ) {\displaystyle \varphi ({\boldsymbol {r}})\ {\stackrel {\mathrm {def} }{=}}\ \varphi _{\mathrm {sr} }({\boldsymbol {r}})+\varphi _{\mathrm {lr} }({\boldsymbol {r}})} .

粒子・メッシュ・エバルト和の基本的考えは、点粒子間の相互作用エネルギーの直接和

E TOT = ∑ i , j φ ( r j − r i ) = E s r + E l r {\displaystyle E_{\text{TOT}}=\sum _{i,j}\varphi ({\boldsymbol {r}}_{j}-{\boldsymbol {r}}_{i})=E_{\mathrm {sr} }+E_{\mathrm {lr} }}

を実空間における短距離ポテンシャルの直接和 E s r {\displaystyle E_{\mathrm {sr} }} （すなわち粒子・メッシュ・エバルトの粒子部分）

E s r = ∑ i , j φ s r ( r j − r i ) {\displaystyle E_{\mathrm {sr} }=\sum _{i,j}\varphi _{\mathrm {sr} }({\boldsymbol {r}}_{j}-{\boldsymbol {r}}_{i})}

と長距離部分のフーリエ空間における和

E l r = ∑ k Φ ~ l r ( k ) | ρ ~ ( k ) | 2 {\displaystyle E_{\mathrm {lr} }=\sum _{\boldsymbol {k}}{\tilde {\Phi }}_{\mathrm {lr} }({\boldsymbol {k}})\left|{\tilde {\rho }}({\boldsymbol {k}})\right|^{2}}

へと置き換えることである。 Φ ~ ℓ r {\displaystyle {\tilde {\Phi }}_{\ell r}} および ρ ~ ( k ) {\displaystyle {\tilde {\rho }}({\boldsymbol {k}})} はポテンシャルおよび電荷密度のフーリエ変換を表す（すなわちエバルト部分）。どちらの和もそれぞれの空間（実空間ならびにフーリエ空間）において素早く収束するため、精度の損失がほとんどなく打ち切ることができ、必要な計算時間を大きく改善することができる。電荷密度場のフーリエ変換 ρ ~ ( k ) {\displaystyle {\tilde {\rho }}({\boldsymbol {k}})} を効率的に評価するため、高速フーリエ変換が用いられる。高速フーリエ変換では、密度場は空間中の離散格子（すなわちメッシュ部分）上で評価される必要がある。

エバルト和はポテンシャルの周期性を仮定する。PME法の物理系への適用にあたり，ポテンシャルの周期的な対称性が必要となる。そのため，この手法は空間的に無限に広がる系(バルク固体など)のシミュレーションに適している。分子動力学シミュレーションでは、これは電荷が中性の単位セルを無限に並べることによって通常達成される。しかしながら、この近似の効果を適切に説明するために、無限に続くセルは元のシミュレーションセルへと再取り込みされる。これは周期的境界条件と呼ばれる。

密度場のメッシュへの制限は、密度の変動が「滑らか」な系（連続ポテンシャル関数を持つ系）についてPME法をより効率的にしている。局在した系、または密度の揺らぎが大きな系については、GreengardとRokhlinの高速多重極法（英語版）を用いてより効率的に扱うことができる。

脚注

  1. ^ Di Pierro, M.; Elber, R.; Leimkuhler, B. (2015), “A Stochastic Algorithm for the Isobaric-Isothermal Ensemble with Ewald Summations for all Long Range Forces.”, Journal of Chemical Theory and Computation, doi:10.1021/acs.jctc.5b00648 .

関連項目