ポール・ピーター・エヴァルトにちなんで名付けられた計算方法
エワルド和は、 ポール・ピーター・エワルド にちなんで名付けられた、 周期系における 長距離相互作用（ 静電相互作用など）を計算する方法です。これは、 イオン結晶の静電エネルギーを計算する方法として最初に開発され、現在では 計算化学 における長距離相互作用の計算に一般的に使用されています。エワルド和は ポアソン和公式 の特殊なケースであり 、実空間での相互作用エネルギーの合計を フーリエ空間 での同等の合計に置き換えます。この方法では、長距離相互作用は、短距離寄与と、 特異点 を持たない長距離寄与の2つの部分に分割されます。短距離寄与は実空間で計算され、長距離寄与は フーリエ変換を 使用して計算されます。この方法の利点は、直接和に比べてエネルギーの 収束が 速いことです。 これは、この方法が長距離相互作用を計算する際に高い精度と妥当な速度を持つことを意味し、周期系における長距離相互作用を計算するための事実上の標準方法となっています。この方法では、総クーロン相互作用を正確に計算するために、分子系の電荷中性が必要です。無秩序な点電荷システムのエネルギーと力の計算で導入される打ち切り誤差に関する研究は、Kolafa と Perram によって提供されています。 [1]
導出 エワルド和は、相互作用ポテンシャルを 2 つの項の和として書き換えます。
ここで、 は 、和が実空間で急速に収束する短距離項を表し、 は、 和がフーリエ (逆数) 空間で急速に収束する長距離項を表します。長距離部分はすべての引数に対して有限である必要があります (最も顕著なのは r = 0) が、便利な数学的形式 (最も一般的には ガウス分布 ) を持つことができます 。この方法では、短距離部分は簡単に合計できることを前提としています。したがって、問題は長距離項の合計になります。フーリエ和を使用しているため、この方法では、調査対象のシステムが無限 周期的で あると暗黙的に想定されています(結晶の内部に対する合理的な想定)。この仮想周期システムの 1 つの繰り返し単位は、 単位セル と呼ばれます。そのようなセルの 1 つが参照用の「中心セル」として選択され、残りのセルは イメージ と呼ばれます。
φ （
r
）
=
d e ふ
φ
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つの電荷密度フィールドの
二重 積分として表すことができます
。単位セルの電荷密度フィールドは、 中心の単位セル内の
電荷の 位置の合計であり
、 総 電荷密度フィールドは、単位セルの電荷 とその周期的なイメージ
の合計です。
え
ℓ r
= ∬ d
r
d
r
′
ρ
トータル
（
r
）
ρ
あなた 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} ^{\プライム })}
ρ
あなた c
（
r
）
{\displaystyle \rho _{uc}(\mathbf {r} )}
r
け
{\displaystyle \mathbf {r} _{k}}
q
け
{\displaystyle q_{k}}
ρ
あなた c
（
r
）
=
d e ふ
∑
c h 1つの r グ e s
け
q
け
δ （
r
−
r
け
）
{\displaystyle \rho _{uc}(\mathbf {r} )\ {\stackrel {\mathrm {def} }{=}}\ \sum _{\mathrm {charges} \ k}q_{k}\delta (\mathbf {r} -\mathbf {r} _{k})}
ρ
トータル
（
r
）
{\displaystyle \rho _{\text{TOT}}(\mathbf {r} )}
q
け
{\displaystyle q_{k}}
ρ
トータル
（
r
）
=
d e ふ
∑
ん
1
、
ん
2
、
ん
3
∑
c h 1つの r グ e s
け
q
け
δ （
r
−
r
け
−
ん
1
1つの
1
−
ん
2
1つの
2
−
ん
3
1つの
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) 数学的な合計が 発散するために発生します。
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 月). 「点電荷システムの Ewald 総和公式におけるカットオフ誤差」. 分子シミュレーション . 9 (5): 351–368. doi :10.1080/08927029208049126.
^ Di Pierro, M.; Elber, R.; Leimkuhler, B. (2015)、「すべての長距離力に対するエワルド和による等圧等温アンサンブルの確率的アルゴリズム」、 化学理論と計算ジャーナル 、 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). 「粒子メッシュエワルド：大規模システムにおけるエワルド和の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。 Bibcode :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。
