Computation method named after Paul Peter Ewald
エワルド和は、 ポール・ピーター・エワルド にちなんで名付けられた、 周期系における 長距離相互作用（ 静電相互作用など）を計算する手法です。これはもともと イオン結晶の静電エネルギーを計算する手法として開発され、現在では 計算化学 における長距離相互作用の計算に広く用いられています 。エワルド和は ポアソン和公式 の特殊なケースであり、実空間における相互作用エネルギーの総和を フーリエ空間 における同等の総和に置き換えたものです。この手法では、長距離相互作用は短距離寄与と 特異点を 持たない長距離寄与の2つの部分に分けられます。短距離寄与は実空間で計算され、長距離寄与は フーリエ変換 を用いて計算されます。この手法の利点は、 直接和を計算する場合に比べてエネルギーが急速に 収束すること です。これは、この方法が長距離相互作用の計算において高い精度と妥当な速度を有することを意味し、周期系における長距離相互作用の計算における 事実上の標準 法となっています。この方法では、全クーロン相互作用を正確に計算するために、分子系の電荷が中性であることが必要です。無秩序な点電荷系のエネルギーと力の計算で導入される打ち切り誤差に関する研究は、KolafaとPerramによって提供されています [1] 。
導出 エワルド和は、相互作用ポテンシャルを2つの項の和として書き換えます。 ここで、 は実空間で和が急速に収束する短距離項、 は フーリエ（逆関数）空間で和が急速に収束する長距離項を表します。長距離項はすべての引数に対して有限（特に r = 0）ですが、任意の便利な数学的形式（最も典型的には ガウス分布） を取ることができます。この手法では、短距離項は容易に和を求めることができると仮定しているため、問題は長距離項の和を求めることになります。フーリエ和を用いるため、この手法では、対象とする系が無限 周期的 であると暗黙的に仮定しています（結晶内部では妥当な仮定です）。この仮想的な周期系の1つの繰り返し単位は 単位セル と呼ばれます。このようなセルの1つが参照用の「中心セル」として選択され、残りのセルは イメージ と呼ばれます。 φ ( r ) = d e f φ s r ( r ) + φ ℓ r ( r ) , {\displaystyle \varphi (\mathbf {r} )\ {\stackrel {\mathrm {def} }{=}}\ \varphi _{sr}(\mathbf {r} )+\varphi _{\ell r}(\mathbf {r} ),} φ s r ( r ) {\displaystyle \varphi _{sr}(\mathbf {r} )} φ ℓ r ( r ) {\displaystyle \varphi _{\ell r}(\mathbf {r} )}
長距離相互作用エネルギーは、中心単位胞の電荷と格子全体の電荷との間の相互作用エネルギーの和である。したがって、これは 単位胞と結晶格子の場を表す2つの電荷密度場にわたる 二 重積分として表すことができる
。ここで、単位胞の電荷密度場は 中心単位胞内の 電荷の 位置にわたる和で あり、 全 電荷密度場は 単位胞の電荷 とその周期像にわたる
同じ和である。 E ℓ r = ∬ d r d r ′ ρ TOT ( r ) ρ u c ( r ′ ) φ ℓ r ( r − r ′ ) {\displaystyle E_{\ell r}=\iint d\mathbf {r} \,d\mathbf {r} ^{\prime }\,\rho _{\text{TOT}}(\mathbf {r} )\rho _{uc}(\mathbf {r} ^{\prime })\ \varphi _{\ell r}(\mathbf {r} -\mathbf {r} ^{\prime })} ρ u c ( r ) {\displaystyle \rho _{uc}(\mathbf {r} )} r k {\displaystyle \mathbf {r} _{k}} q k {\displaystyle q_{k}} ρ u c ( r ) = d e f ∑ c h a r g e s k q k δ ( r − r k ) {\displaystyle \rho _{uc}(\mathbf {r} )\ {\stackrel {\mathrm {def} }{=}}\ \sum _{\mathrm {charges} \ k}q_{k}\delta (\mathbf {r} -\mathbf {r} _{k})} ρ TOT ( r ) {\displaystyle \rho _{\text{TOT}}(\mathbf {r} )} q k {\displaystyle q_{k}} ρ TOT ( r ) = d e f ∑ n 1 , n 2 , n 3 ∑ c h a r g e s k q k δ ( r − r k − n 1 a 1 − n 2 a 2 − n 3 a 3 ) {\displaystyle \rho _{\text{TOT}}(\mathbf {r} )\ {\stackrel {\mathrm {def} }{=}}\ \sum _{n_{1},n_{2},n_{3}}\sum _{\mathrm {charges} \ k}q_{k}\delta (\mathbf {r} -\mathbf {r} _{k}-n_{1}\mathbf {a} _{1}-n_{2}\mathbf {a} _{2}-n_{3}\mathbf {a} _{3})}
ここで、 は ディラックのデルタ関数 、、 は 格子ベクトル、であり 、 すべて の 整数にわたる。全場は 格子関数 との 畳み込み として 表される。 δ ( x ) {\displaystyle \delta (\mathbf {x} )} a 1 {\displaystyle \mathbf {a} _{1}} a 2 {\displaystyle \mathbf {a} _{2}} a 3 {\displaystyle \mathbf {a} _{3}} n 1 {\displaystyle n_{1}} n 2 {\displaystyle n_{2}} n 3 {\displaystyle n_{3}} ρ TOT ( r ) {\displaystyle \rho _{\text{TOT}}(\mathbf {r} )} ρ u c ( r ) {\displaystyle \rho _{uc}(\mathbf {r} )} L ( r ) {\displaystyle L(\mathbf {r} )} L ( r ) = d e f ∑ n 1 , n 2 , n 3 δ ( r − n 1 a 1 − n 2 a 2 − n 3 a 3 ) {\displaystyle L(\mathbf {r} )\ {\stackrel {\mathrm {def} }{=}}\ \sum _{n_{1},n_{2},n_{3}}\delta (\mathbf {r} -n_{1}\mathbf {a} _{1}-n_{2}\mathbf {a} _{2}-n_{3}\mathbf {a} _{3})}
これは畳み込み なので 、 の フーリエ変換 は積になります。 ここで、格子関数のフーリエ変換は、 逆空間ベクトルが定義されているデルタ関数 （および巡回置換）の別の和です。ここで、 は中心単位胞の体積です（幾何学的に 平行六面体 である場合、多くの場合はそうなりますが、必ずしもそうとは限りません）。 と はどちら も 実関数で偶関数であることに注意してください。 ρ TOT ( r ) {\displaystyle \rho _{\text{TOT}}(\mathbf {r} )} ρ ~ TOT ( k ) = L ~ ( k ) ρ ~ u c ( k ) {\displaystyle {\tilde {\rho }}_{\text{TOT}}(\mathbf {k} )={\tilde {L}}(\mathbf {k} ){\tilde {\rho }}_{uc}(\mathbf {k} )} L ~ ( k ) = ( 2 π ) 3 Ω ∑ m 1 , m 2 , m 3 δ ( k − m 1 b 1 − m 2 b 2 − m 3 b 3 ) {\displaystyle {\tilde {L}}(\mathbf {k} )={\frac {\left(2\pi \right)^{3}}{\Omega }}\sum _{m_{1},m_{2},m_{3}}\delta (\mathbf {k} -m_{1}\mathbf {b} _{1}-m_{2}\mathbf {b} _{2}-m_{3}\mathbf {b} _{3})} b 1 = d e f 2 π a 2 × a 3 Ω {\displaystyle \mathbf {b} _{1}\ {\stackrel {\mathrm {def} }{=}}\ 2\pi {\frac {\mathbf {a} _{2}\times \mathbf {a} _{3}}{\Omega }}} Ω = d e f a 1 ⋅ ( a 2 × a 3 ) {\displaystyle \Omega \ {\stackrel {\mathrm {def} }{=}}\ \mathbf {a} _{1}\cdot \left(\mathbf {a} _{2}\times \mathbf {a} _{3}\right)} L ( r ) {\displaystyle L(\mathbf {r} )} L ~ ( k ) {\displaystyle {\tilde {L}}(\mathbf {k} )}
簡潔にするために、有効な単一粒子ポテンシャルを定義する。 v ( r ) = d e f ∫ d r ′ ρ u c ( r ′ ) φ ℓ r ( r − r ′ ) {\displaystyle v(\mathbf {r} )\ {\stackrel {\mathrm {def} }{=}}\ \int d\mathbf {r} ^{\prime }\,\rho _{uc}(\mathbf {r} ^{\prime })\ \varphi _{\ell r}(\mathbf {r} -\mathbf {r} ^{\prime })}
これも畳み込みなので、同じ方程式のフーリエ変換は、 フーリエ変換が定義される
積となる。 V ~ ( k ) = d e f ρ ~ u c ( k ) Φ ~ ( k ) {\displaystyle {\tilde {V}}(\mathbf {k} )\ {\stackrel {\mathrm {def} }{=}}\ {\tilde {\rho }}_{uc}(\mathbf {k} ){\tilde {\Phi }}(\mathbf {k} )} V ~ ( k ) = ∫ d r v ( r ) e − i k ⋅ r {\displaystyle {\tilde {V}}(\mathbf {k} )=\int d\mathbf {r} \ v(\mathbf {r} )\ e^{-i\mathbf {k} \cdot \mathbf {r} }}
エネルギーは 単一の 場の積分
として表すことができる。 E ℓ r = ∫ d r ρ TOT ( r ) v ( r ) {\displaystyle E_{\ell r}=\int d\mathbf {r} \ \rho _{\text{TOT}}(\mathbf {r} )\ v(\mathbf {r} )}
プランシュレルの定理 を用いると 、エネルギーはフーリエ空間でも合計できる。 E ℓ r = ∫ d k ( 2 π ) 3 ρ ~ TOT ∗ ( k ) V ~ ( k ) = ∫ d k ( 2 π ) 3 L ~ ∗ ( k ) | ρ ~ u c ( k ) | 2 Φ ~ ( k ) = 1 Ω ∑ m 1 , m 2 , m 3 | ρ ~ u c ( k ) | 2 Φ ~ ( k ) {\displaystyle E_{\ell r}=\int {\frac {d\mathbf {k} }{\left(2\pi \right)^{3}}}\ {\tilde {\rho }}_{\text{TOT}}^{*}(\mathbf {k} ){\tilde {V}}(\mathbf {k} )=\int {\frac {d\mathbf {k} }{\left(2\pi \right)^{3}}}{\tilde {L}}^{*}(\mathbf {k} )\left|{\tilde {\rho }}_{uc}(\mathbf {k} )\right|^{2}{\tilde {\Phi }}(\mathbf {k} )={\frac {1}{\Omega }}\sum _{m_{1},m_{2},m_{3}}\left|{\tilde {\rho }}_{uc}(\mathbf {k} )\right|^{2}{\tilde {\Phi }}(\mathbf {k} )}
最終的な合計の どこに。 k = m 1 b 1 + m 2 b 2 + m 3 b 3 {\displaystyle \mathbf {k} =m_{1}\mathbf {b} _{1}+m_{2}\mathbf {b} _{2}+m_{3}\mathbf {b} _{3}}
これが本質的な結果です。 を計算してしまえば、 上の和／積分は 簡単で、すぐに収束するはずです。収束しない最も一般的な理由は、単位胞の定義が不十分なことです。無限和を避けるためには、単位胞は電荷が中性でなければなりません。 ρ ~ u c ( k ) {\displaystyle {\tilde {\rho }}_{uc}(\mathbf {k} )} k {\displaystyle \mathbf {k} }
粒子メッシュエヴァルト（PME）法 エワルド和法は、コンピュータ が登場するはるか以前、 理論物理学 における手法として開発されました 。しかし、エワルド法は1970年代以降、粒子系、特に 重力 や 静電気 のように 逆二乗 則 によって相互作用する粒子系の コンピュータシミュレーションにおいて広く利用されてきました。最近では、PMEは レナード・ジョーンズポテンシャル の部分 を計算する際にも使用され、 打ち切りによるアーティファクトを排除しています。 [2] プラズマ 、 銀河 、 分子 のシミュレーションなどに応用されています 。 r − 6 {\displaystyle r^{-6}}
粒子メッシュ 法では 、標準的なエワルド和と同様に、一般的な相互作用ポテンシャルは2つの項に分離されます 。粒子メッシュエワルド和の基本的な考え方は、点粒子間の相互作用エネルギーの直接和を
、 実空間における短距離ポテンシャルの 直接和 （これは 粒子メッシュエワルドの 粒子 部分です ）とフーリエ空間における長距離部分の和の2つの和に置き換えることです。 φ ( r ) = d e f φ s r ( r ) + φ ℓ r ( r ) {\displaystyle \varphi (\mathbf {r} )\ {\stackrel {\mathrm {def} }{=}}\ \varphi _{sr}(\mathbf {r} )+\varphi _{\ell r}(\mathbf {r} )} E TOT = ∑ i , j φ ( r j − r i ) = E s r + E ℓ r {\displaystyle E_{\text{TOT}}=\sum _{i,j}\varphi (\mathbf {r} _{j}-\mathbf {r} _{i})=E_{sr}+E_{\ell r}} E s r {\displaystyle E_{sr}} E s r = ∑ i , j φ s r ( r j − r i ) {\displaystyle E_{sr}=\sum _{i,j}\varphi _{sr}(\mathbf {r} _{j}-\mathbf {r} _{i})} E ℓ r = ∑ k Φ ~ ℓ r ( k ) | ρ ~ ( k ) | 2 {\displaystyle E_{\ell r}=\sum _{\mathbf {k} }{\tilde {\Phi }}_{\ell r}(\mathbf {k} )\left|{\tilde {\rho }}(\mathbf {k} )\right|^{2}}
ここで 、およびは ポテンシャル と 電荷密度 の フーリエ変換を 表します （これは エワルド 部分です）。両方の和はそれぞれの空間（実空間とフーリエ空間）で急速に収束するため、精度の低下をほとんど伴わず、必要な計算時間を大幅に短縮しながら切り捨てることができます。 電荷密度場のフーリエ変換を効率的に評価するには、 高速フーリエ変換を 使用します。この変換では、密度場を空間内の離散格子（これは メッシュ 部分です）上で評価する必要があります。 Φ ~ ℓ r {\displaystyle {\tilde {\Phi }}_{\ell r}} ρ ~ ( k ) {\displaystyle {\tilde {\rho }}(\mathbf {k} )} ρ ~ ( k ) {\displaystyle {\tilde {\rho }}(\mathbf {k} )}
エワルド和に暗黙的に含まれる周期性仮定のため、PME法を物理系に適用するには周期対称性を課す必要があります。したがって、この方法は空間的に無限に広がる系としてシミュレートできる系に最適です。 分子動力学 シミュレーションでは、通常、電荷中性の単位セルを意図的に構築し、これを無限に「タイル状に並べることで」像を形成させます。しかし、この近似の効果を適切に考慮するために、これらの像は元のシミュレーションセルに再統合されます。この全体的な効果は 周期境界条件 と呼ばれます。これを最も明確に視覚化するには、 単位立方体 を想像してみてください。上面は下面と、右面は左面と、前面は背面と実質的に接触しています。したがって、単位セルのサイズは、「接触」している2つの面間の不適切な運動相関を回避するのに十分な大きさでありながら、計算上実行可能な程度に小さいサイズになるように慎重に選択する必要があります。短距離相互作用と長距離相互作用のカットオフの定義によっても、アーティファクトが発生する可能性があります。
密度場をメッシュに制限することで、PME法は密度が「滑らかに」変化する系、あるいは連続的なポテンシャル関数を持つ系に対してより効率的になります。局所的な系や密度が大きく変動する系は、 GreengardとRokhlinの 高速多重極法を用いるとより効率的に扱えます。
双極子項 極性結晶（すなわち、単位胞内に 正味の双極子を持つ結晶）の静電エネルギーは 条件付き収束 、すなわち和の順序に依存します。例えば、中心の単位胞と、常に増加する立方体上に配置された単位胞との双極子間相互作用が、相互作用エネルギーを球状に和した場合とは異なる値に収束する場合を例に挙げます。大まかに言えば、この条件付き収束は、(1) 半径の殻上で相互作用する双極子の数が のように増加する 、(2) 単一の双極子間相互作用の強度が のように低下する 、(3) 数学的な和 が発散する、という3つの理由から生じます。 p u c {\displaystyle \mathbf {p} _{uc}} R {\displaystyle R} R 2 {\textstyle R^{2}} 1 / R 3 {\textstyle 1/{R^{3}}} ∑ n = 1 ∞ 1 n {\textstyle \sum _{n=1}^{\infty }{\frac {1}{n}}}
このいくぶん意外な結果は、実際の結晶の有限なエネルギーと調和できます。なぜなら、そのような結晶は無限ではなく、特定の境界を持っているからです。より具体的には、極性結晶の境界には、表面に有効表面電荷密度があります。 ここで、 は表面法線ベクトル、 体積あたりの正味双極子モーメントを表します。 その表面電荷密度を持つ中心単位胞内の双極子の相互作用エネルギーは次のように書き表すことができます [3]。 ここで 、 と は単位胞の正味双極子モーメントと体積、 は結晶表面上の微小領域、 は 中心単位胞から微小領域へのベクトルです。この式は、エネルギーを積分することによって得られます。 ここで、 は 微小表面電荷によって生成される微小電場を表します ( クーロンの法則 )。 負の符号は の定義に由来し 、電荷に向かう方向を指し、電荷から離れる方向を指しません。 σ = P ⋅ n {\displaystyle \sigma =\mathbf {P} \cdot \mathbf {n} } n {\displaystyle \mathbf {n} } P {\displaystyle \mathbf {P} } U {\displaystyle U} U = 1 2 V u c ∫ ( p u c ⋅ r ) ( p u c ⋅ n ) r 3 d S {\displaystyle U={\frac {1}{2V_{uc}}}\int {\frac {\left(\mathbf {p} _{uc}\cdot \mathbf {r} \right)\left(\mathbf {p} _{uc}\cdot \mathbf {n} \right)}{r^{3}}}\,dS} p u c {\displaystyle \mathbf {p} _{uc}} V u c {\displaystyle V_{uc}} d S {\displaystyle dS} r {\displaystyle \mathbf {r} } d U = − p u c ⋅ d E {\displaystyle dU=-\mathbf {p} _{uc}\cdot d\mathbf {E} } d E {\displaystyle d\mathbf {E} } d q = d e f σ d S {\displaystyle dq\ {\stackrel {\mathrm {def} }{=}}\ \sigma dS} d E = d e f ( − 1 4 π ϵ ) d q r r 3 = ( − 1 4 π ϵ ) σ d S r r 3 {\displaystyle d\mathbf {E} \ {\stackrel {\mathrm {def} }{=}}\ \left({\frac {-1}{4\pi \epsilon }}\right){\frac {dq\ \mathbf {r} }{r^{3}}}=\left({\frac {-1}{4\pi \epsilon }}\right){\frac {\sigma \,dS\ \mathbf {r} }{r^{3}}}} r {\displaystyle \mathbf {r} }
歴史 エワルド和は、 イオン結晶の 静電エネルギー（および マーデルング定数）を決定するために、1921 年に ポール ピーター エワルドによって開発されました（下記の参考文献を参照）。
スケーリング 一般的に、異なるエワルド和法では 計算時間が 異なります。直接計算では （ は系内の原子数）となります。PME法では となります 。 [4] O ( N 2 ) {\displaystyle O(N^{2})} N {\displaystyle N} O ( N log N ) {\displaystyle O(N\,\log N)}
参照
参考文献 ^ Kolafa, Jiri; Perram, John W. (1992年9月). 「点電荷系におけるエワルド和公式のカットオフ誤差」. 分子シミュレーション . 9 (5): 351– 368. doi :10.1080/08927029208049126. ^ Di Pierro, M.; Elber, R.; Leimkuhler, B. (2015)「すべての長距離力に対するエワルド和を用いた等圧等温アンサンブルの確率的アルゴリズム」 Journal of Chemical Theory and Computation 、 11 (12): 5624– 5637、 doi :10.1021/acs.jctc.5b00648、 PMC 4890727 、 PMID 26616351 ^ Herce, HD; Garcia, AE; Darden, T (2007年3月28日). 「静電表面項：(I)周期系」. The Journal of Chemical Physics . 126 (12): 124106. Bibcode :2007JChPh.126l4106H. doi :10.1063/1.2714527. PMID 17411107. ^ Darden, Tom; York, Darrin; Pedersen, Lee (1993-06-15). 「粒子メッシュEwald：大規模システムにおけるEwald和を求めるN⋅log(N)法」. The Journal of Chemical Physics . 98 (12): 10089– 10092. Bibcode :1993JChPh..9810089D. doi :10.1063/1.464397. ISSN 0021-9606. エワルド、P (1921)。 「Die Berechnung optischer und elektrostatischer Gitterpotentiale」。 アンナレン・デア・フィジーク 。 369 (3): 253–287 。 書誌コード :1921AnP...369..253E。 土井 ：10.1002/andp.19213690304。 Darden, T; Perera, L; Li, L; Pedersen, L (1999). 「結晶学ツールキットからモデラーへの新しいテクニック：粒子メッシュEwaldアルゴリズムと核酸シミュレーションへの応用」 Structure . 7 (3): R55 – R60 . doi : 10.1016/S0969-2126(99)80033-1 . PMID 10368306. S2CID 40964921. Frenkel, D., Smit, B. (2001). 『分子シミュレーションの理解：アルゴリズムから応用まで 』Academic press. 
