《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）

平成 27 年度 修士論文

短距離古典分子動力学法における加速化手法の系統的構築：
近接リスト-セル分割ハイブリッド法とリスト半径・セルサイズの最適化

新潟大学大学院自然科学研究科 博士前期課程 2 年
数理物質科学専攻

物理学コース

F14A001B
淡路

大輔

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）

概要
短距離古典分子動力学 (Molecular Dynamics: MD) 法における計算時間短縮の手法として，相互作用する近傍粒子を効率よく探索する方法である Verlet の近接リスト法 (Verlet Neighbor List: VNL 法) とセル分割法 (Cell Linked List: CLL 法) は従来からよく用いられている．本論文では両者を組み合わせ
たハイブリッドな手法 (VNL-CLL 法) を提案し，その最適化の方法について述べている．VNL 法はカットオフ半径より少し大きなリスト半径を設けておき，その範囲内に存在する近傍粒子をリストに登録，各MD ステップごとにそのリストに含まれる近傍粒子との相互作用を計算する方法である．そこでは粒子配置が時間発展によって変化するため数ステップ毎にリストを更新しなければならない．一般にはこの更新
頻度は 10 ステップ程度に固定することが多いが，リストを動的に更新する手法を採用することにより計
算時間の短縮を図ることができる．この動的更新型の VNL 法では，計算システム (粒子数，温度，密度，
カットオフ半径) に対して計算コストを最適とするリスト半径がある事が従来までに分かっているが，さらなる高速化のためには各パラメータ依存性の検討が必要である．一方，CLL 法では計算の対称としている系全体を一辺がカットオフ半径よりもわずかに大きい立方体セルに分割する．これにより粒子間相互作用は，同じセルに属する粒子間，もしくは隣接セルに属する粒子間に限定される．隣接するセルの数は限られており (2 次元で 8 個，3 次元で 26 個)，相互作用の数を大幅に減らすことが出来る．さらにセル分割数を増加させ，より細分化したセルを用いることによって参照体積を減らす．これにより従来型の CLL 法と比べて，計算の加速化を図ることができる．この手法は排他的セル分割法と呼ばれ，相互作
用 (カットオフ半径) が非常に短距離的である粉体のような系に対して大きな効果がある．VNL 法において通常は全粒子探索によりリストを再構築するが，これには O(N 2 ) の計算コストを要する．しかしながら，リスト更新に CLL 法を用いて近傍粒子探索をすることによりその計算コストを O(N 2 ) から O(N )まで低減することができる．この手法を VNL-CLL 法と呼ぶ．本研究では VNL-CLL 法において Verletの近接リスト更新の際に排他的セル分割法を用い，計算時間のセルサイズ (あるいはセル分割数) 依存性を系統的に明らかにした. また，近接リストに対する動的更新の手法を導入し，VNL 法の結果に基づき，計算時間のリスト半径依存性を計算システム (粒子数，温度，密度，カットオフ半径) を変化させ，網羅的に検討した．その結果，計算時間を最小化する最適なリスト半径，セルサイズの組み合わせが計算システ
ムに依存すること見出し，さらなる加速化のためにリスト半径，セルサイズの最適化の手法を構築した．

–i–《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）

目次
i

概要
第1章

序論

1

1.1

分子間相互作用 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11.2

運動方程式の数値積分 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31.3

周期境界条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

1.4

古典統計力学の基礎理論 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

1.5

本論文の目的と構成 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

ベルレの近接リスト法

92.1

概要 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

92.2

SKIN(リスト半径) の最適化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

セル分割法

17

3.1

概要 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

17

3.2

排他的セル分割法. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

20

3.3

セルサイズの最適化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23

VNL-CLL 法

264.1

概要 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

4.2

リスト半径・セルサイズの最適化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

29

結論

34

第2章

第3章

第4章

第5章
謝辞

36

参考文献

37

– ii –

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）

第1章

序論
分子動力学法 (Molecular Dynamics: MD 法) はモンテカルロ法 (Monte Carlo: MC 法) と並ぶ分子シ
ミュレーションの手法の一つである．MC 法はハミルトニアンから，系の状態を実現するボルツマン分布に従う粒子配置を乱数を用いて求め，平衡状態における熱力学量のカノニカル平均を計算する手法である．これに対し，MD 法は運動方程式を数値積分して粒子運動の時間発展を求める．そのため熱平衡状態における熱力学量ばかりではなく，原子・分子の動的挙動についても有用な情報を得る事が可能である．
さらに，外場があるような非平衡状態における粒子系についてもシミュレーションを行うことができる利点がある．本章ではまず，MD 法の概要について紹介する [1–4]．

1.1 分子間相互作用
MD 法と MC 法，どちらを用いて分子シミュレーションを行うにしても，粒子間に働く相互作用を設定する．希ガス原子間に働く力は 2 体相互作用 (対ポテンシャル， pair potential) によってよく記述され，希ガス以外に単純な金属やイオン系に対しても用いられる．もし，構造をもつ分子を考えれば 3 体
力，4 体力が無視できなくなる．
質量が m である単成分球状 N 粒子系のポテンシャルエネルギーを UN ({ri }) と書く事にする．ri は i番目の粒子の座標ベクトルであり，{ri } はすべての粒子の座標を与えて初めて決まる関数である事を示

3

φ(r) / ε

2
1
0
-1
0

1

2

3

r/σ
図 1.1

Lennard-Jones(LJ) ポテンシャル：青線が LJ ポテンシャル，赤線は斥力項，緑線は引力項を表す．

–1–

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）

す．ここで，粒子間距離 r にのみ依存する古典的 2 体相互作用ポテンシャル (pair potential) を φ(r) とすれば，全ポテンシャルエネルギーはN

UN ({ri }) =N

1 φ(|ri − rj |)
2 i=1

(1.1)

i=j

である．代表例としてクーロン力や分散力などが挙げられる．係数 1/2 は (i，j) ペアを 2 重に数える事を防ぐためのものである．ポテンシャルエネルギー UN から，i 番目粒子に働く力をFi = −
=−

∂UN
∂ri
N
1  ∂φ(rij )
2

j=i

(1.2)

∂ri


N 
1  ∂φ(rij ) ∂φ(rji )
=
+
2
∂ri
∂ri
j=i

=−

N

dφ(rij ) ∂rij
j=i

=−

drij

∂ri

N

dφ(rij ) ri − rj
j=i

drij

rij

と与えることができる．ここで，(ri − rj )/rij は，j 粒子から i 粒子に向かう単位ベクトルを表している．

2 体ポテンシャルの代表例として Lennard-Jones(LJ) ポテンシャルを図 1.1 に示す．LJ ポテンシャル
は，アルゴンやキセノンなどの希ガスのファンデルワールス力に対するポテンシャルとして用いられ σ 12  σ 6φ(r) = 4
−
r
r

(1.3)

と与えられる． はポテンシャル深さであり，σ は粒子直径である．第 1 項は斥力項，第 2 項は引力項であるが，r < σ ではポテンシャルの勾配が急峻で強い斥力が働き，r −1 で減衰するクーロンポテンシャルに比べて非常に短距離的な相互作用である．粒子間に働く力は，

F = −∇φ(r) = −24  σ 13  σ 7 r
dφ(r) r
=
2
−
dr r
σ
r
r
r

(1.4)

となる．

LJ ポテンシャルなどの短距離相互作用では、一般に粒子間距離が粒子直径の 3 倍程度になるとその相互作用は急激に小さくなる。そこで、シミュレーションの都合上カットオフする半径を設定し、それより遠方にある粒子との間に働く力をゼロにする。これを相互作用のカットオフという。これにより、カット
オフ半径内に含まれない遠方の粒子との相互作用を計算しなくても良くなるが、そもそもカットオフ半径内に含まれるか否かを判定するために粒子間距離の計算は必須である。つまり依然として粒子間距離 (あるいは粒子ペア) の計算に O(N 2 ) の計算コストを要してしまう。この計算コストを O(N 2 ) から O(N )に低減するアルゴリズムとして，近接粒子の探索手法がいくつか提案されている．本研究においてはこの近接粒子探索手法に焦点を当てている．

–2–

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）

1.2 運動方程式の数値積分
相互作用ポテンシャル φ(r) を LJ ポテンシャルによって与えられた粒子系などで，ニュートンの運動
方程式

d2 r
F (t)
=
2
dt
m

(1.5)

を数値的に解く方法が MD 法である．実用的な数値計算においては，比較的簡単な数値積分法が使われることが多い．その中のひとつが Verlet の差分式 [1–5] と呼ばれる，式 (1.5) を時間刻み Δt に対して 2
次精度を持つ差分式である．時刻 t + Δt と t − Δt における粒子の位置 r(t ± Δt) をテイラー展開すると，r(t + Δt) = r(t) + Δtv(t) +(Δt)2 F (t)+ O((Δt)3 )
2
m

(1.6)

r(t − Δt) = r(t) − Δtv(t) +(Δt)2 F (t)+ O((Δt)3 )
2
m

(1.7)

を得る．両辺の和と差をとるとr(t + Δt) + r(t − Δt) = 2r(t) + (Δt)2F (t)+ O((Δt)4 )mr(t + Δt) − r(t − Δt) = 2Δtv(t) + O((Δt)3 )

(1.8)
(1.9)

これより，時刻 t + Δt における位置と t における速度はr(t + Δt) = 2r(t) − r(t − Δt) + (Δt)2v(t) =F (t)+ O((Δt)4 )m

1{r(t + Δt) − r(t − Δt)} + O((Δt)2 )2Δt(1.10)(1.11)

で与えられる．これが Verlet の差分式である．時刻 t + Δt における位置を求めるには，２つの時刻 t と

t − Δt での位置が必要である．初期条件を位置と速度で与えると，t = Δt における位置 r(Δt) は式 1.7から求まる．これと r(0) から r(2Δt) を計算し，式 (1.11) より，速度 v(Δt) が得られる．この差分方程式は時間反転 (Δt → −Δt) に対して対称であり，速度の符号を変えると系は描いた軌跡を逆向きに辿って行く．式 (1.11) の代わりに速度を



Δt
v t+
2
で与えると，式 (1.10) は



Δt
v t+
2





=

1
[r(t + Δt) − r(t)]
Δt


Δt
=v t−
2


+ Δt

F (t)
m

(1.12)(1.13)と書き換えられる．これと，式 (1.12) を r(t + Δt) について解いたdtr(t + Δt) = r(t) + Δtv t +2

(1.14)

を組み合わせたのが，よく使われる Verlet の蛙跳び (Verlet’s leap frog) 法である [1–5]．式 (1.13) で最初に速度を修正し, 次に式 (1.14) で位置の更新を行う．

–3–

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）

図 1.2

周期境界条件の概念図：中央の灰色のセルが基本セル，周辺のセルがイメージセルを表している．基本セルは一辺 L = V 1/3 で決まる．V は計算系の全体積である．

1.3 周期境界条件
シミュレーションをおこなう系の物性を調べる際，物質の表面の性質なのか，表面から十分離れた内部(バルク) の状態なのか，あるいは 2 つの相の共存状態なのかなど，対象とする系のおかれている状況に応じて境界条件を設定しなければならない．
現実系のような 1023 個の粒子を取り扱うことは，現在の最新鋭のコンピュータをもってしてもシリア
ル CPU で 104 ∼ 106 ，並列計算により ∼ 109 個程度が限界であり，事実上不可能である．104 ∼ 106 個程度の粒子数のシミュレーションでは界面による影響は顕著であり，限られた粒子数でバルクの性質を抽出することがシミュレーションにおいては必須となる．そこで，バルクのシミュレーションに用いる周期
境界条件 (Periodic Boundary Condition) について説明する．図 1.2 のように，計算システムを x, y, z
方向にそれぞれに繰り返し並べ，以下の条件を考慮することが周期境界条件である．中央のセルが一辺がL = V 1/3 の基本セル，周辺のセルをイメージセルと呼ぶ．ここで V は全体積である．周期境界条件では以下の点が考慮されている．

• 有限個の粒子が有限サイズの計算セルに閉じ込められていると考える．計算セルの形状は，周期的に並べる事が出来るよう，立方体，直方体，または，平行六面体とする．

• ある粒子が時間発展により運動し，位置を変え，右側のイメージセルに移ったとすれば (図 1.2 の矢印)，左隣のイメージセルから同一の粒子が補充される．

• 粒子の相互作用に関しても周期境界を考えなければならない．イメージセルを導入したとしても，実質的に同じペアの相互作用を重ねて数えることがないように，実粒子・shadow 粒子を含めて最–4–

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）

も近いペアの相互作用のみを考慮する必要がある．これは相互作用を計算セルサイズの半分のところでカットする事に相当し，minimum image convention と呼ばれている．
以上３つの条件を考慮し，周期境界条件のもとに計算を実行する．これにより，少ない粒子数のシミュ
レーションでもバルクの性質を抽出することが可能となる．他の境界条件としては，反射境界，自由境界などがある．

1.4 古典統計力学の基礎理論
MD シミュレーションから，各時刻における粒子の位置 ri (t) および速度 vi (t) が得られる．以下では古典統計力学に基づいて粒子の位置・速度から得られる代表的な熱力学量の計算方法に関して概説す
る [6–9]．

(a) エネルギー
運動方程式を数値的に解くと，全エネルギー
N

E=

N

1

N

1 
φ(|ri − rj |)
2 i=1

+

2
i=1

i=j

mvi2

(1.15)

は保存される．つまり，小正準集合 (NVE アンサンブル) を計算していることになる．ここで · · · はアンサンブル平均または時間平均を表している．

(b) 温度
ある時刻における瞬間的な温度 T (t) は，運動エネルギーからT (t) =

N1 
mvi2
3N kB i=1

(1.16)と決める事ができる．つまり，エネルギー等分配則の結果である．さらに充分長い時間をかけて計算し系が熱平衡状態に到達したとき系の温度は瞬間的な温度 T (t) の時間平均から

1
τ →∞ τ

τT = T (t)t = lim

T (t)dt(1.17)

0のように求められる．ここで · · · t は特に時間に対する平均値をとることを意味している．

(c) 圧力圧力は Virial の定理に基づいて計算される．まず，i 番目粒子が容器の壁から受ける力を Wi ，他の粒子からの力の総和を Fi とすれば

d2 r i= F i + Wi
(1.18)
dt2の運動方程式が成り立つ．この両辺の各項と i 番目粒子の位置ベクトル ri との内積をとり，時間m積分をとる．ri · Fi の項は
τ
0


τ
d2 r idriri ·dt = ri ·
−
dt2
dt 0

τ
0

dri dri·dt
dt dt(1.19)のように部分積分を行うと，右辺第 1 項の ri · vi は任意の時刻で有限の値を持つ．したがって τ
で割って τ → ∞ とすれば消える．全ての粒子に対して和をとれば，
N

1
i=1

2


m

dri
dt

2

1
=−
2
–t 5 –N
i=1

ri · (Fi + Wi )(1.20)t《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）となる．ここで左辺は運動エネルギーの長時間平均である．一方で，右辺の量は Clausiu によって名付けられた Virial である．系の運動エネルギーの長時間平均が Virial に等しいことを示すこの関係式は Virial の定理と呼ばれる．さらに，エネルギー等分配則からN


r i · Wii=1N

= −3N kB T −ri · Fi(1.21)i=1ttとすることができ，左辺は壁から受ける力の Virial である．これは圧力 P と関係づけられる．つまり，i 番目粒子が容器の壁に及ぼす力は −Wi であるが，すべての粒子が壁面の単位面積あたりに及ぼす力の平均が圧力である．壁に位置ベクトル r の面積要素 dS をとり，そこから外向き単位
法線ベクトル n とすれば，単位面積 dS が受ける平均の力は P ndS である．したがって Gauss の定理を用いればN

r i · Wi

i=1

=
S

t

P n · rdS = −3P V(1.22)である．よって圧力は

P =

=

1
N kB T+
V
3V
N kB T
1
−
V
3V

N


ri · Fii=1N(1.23)
tri · ∇i UN ({ri })i=1(1.24)tとして計算する事ができる．

(d) 動径分布関数
1 つの粒子を中心とした他の粒子の位置の統計的分布を表す関数として動径分布関数がある．原点(r = 0) に粒子があるとき半径 (r, r + dr) の球殻内にある粒子の平均個数を dn(r) とおくと，平均個数と動径分布関数とは

dn(r) = ρg(r)4πr2 dr,

ρ=

N
V

(1.25)の関係で結ばれている．ここで ρ は数密度である．この動径分布関数は液体構造の静的性質を調べる指標となる．g(r) は r → ∞ で g(r) = 1 となるように規格化され，また排除体積により r → 0
で g(r) = 0 となる．また，この動径分布関数を用いて系のエネルギー，圧力を評価することが可能である．内部エネルギーを

E = E id + E ex

(1.26)のように理想項と過剰項に分ける．ここで理想項はE id =
過剰項は

E ex = UN ({ri }) =1
ZN

3
N kB T
2

UN ({ri }) exp(−UN ({ri })/kB T )dr N(1.27)(1.28)である．またZN =exp(−UN ({ri })/kB T )dr N(1.29)

は配置積分と呼ばれる．内部エネルギーの過剰項は動径分布関数 g(r) を用いて∞E ex
φ(r)g(r)r2 dr= 2πρ
N– 6 –0(1.30)

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）

としてあらわさる．圧力は式 (1.24) より，式 (1.28) と同様の考察から

P
2πρ
=1−
ρkB T
3kB T

∞

φ (r)g(r)r3 dr(1.31)

0とあらわすことができる．

1.5 本論文の目的と構成
分子シミュレーションにおいて最も時間がかかるのは，粒子間相互作用 (ポテンシャルエネルギーや力の計算) の部分である．LJ ポテンシャルのような 2 体相互作用ですら，粒子間距離の計算に O(N 2 ) の計
算コストを要してしまう．これは位置の更新，速度の更新などの O(N ) の計算コストに対し顕著である．
大きなシステムサイズの系に対して，純粋に O(N 2 ) の計算を行うことは現代のコンピュータパフォーマ
ンスを持ってしても現実的でない．本論文の目的は，計算効率を向上するための近接粒子探索リストアル
ゴリズムを提案すると共に，それぞれの手法を最適化する方法を構築する事である．
本論文では 2 つの粒子探索手法を組み合わせた VNL-CLL 法を提案している．その要素はリストの動的更新の手法を用いたベルレの近接リスト法 (Verlet Neighbor List: VNL 法) とセル分割法 (Cell Linked

List: CLL 法) からなる．さらにこの CLL 法にはセルサイズを拡張した排他的セル分割法 (Exclusive
CLL 法) を採用した．VNL 法・CLL 法にはそれぞれ利点・欠点がある．VNL 法は基本的に計算コスト
を O(N ) に低減することができる利点があるが，数ステップに 1 度 O(N 2 ) の計算コストを要する更新を伴う欠点がある．CLL 法はリスト更新の必要がなく常に O(N ) の計算コストに低減することができる一方で，VNL 法と比較しリストに保存しておく粒子数が多い事が欠点である．本研究において新規に開発した VNL-CLL 法は VNL 法・CLL 法両手法の利点を組み合わせることによって，それぞれの欠点を克服している大きな特徴がある．
本論文の構成は以下の通りである．第 2 章では VNL 法とそれをさらに拡張したリストの動的更新の手法を紹介する．第 3 章では VNL 法と同様に良く用いられる手法である CLL 法と，より高度な排他的セル分割法に対し，アルゴリズムと最適化の方法を紹介する．また，対応するパラメータ依存性の検討，計算性能の比較を行った．第 4 章では VNL 法，CLL 法の両者の利点を組み合わせた VNL-CLL 法を提案し，その最適化手法のための計算時間に対する評価式を解析的に与える．第 5 章では結論を述べる．また，最適化された各種法 (VNL 法，リスト動的更新型 VNL 法，CLL 法，Exclusive CLL 法，VNL-CLL法) に関して計算パフォーマンスの比較を行った．
なお，本研究では CPU1 コア計算を実施し，各近接粒子探索手法の効率性について検討している．具体的には LJ ポテンシャル，周期境界条件を用いて NVE アンサンブルのもとで計算した．各手法ともに

100,000 ステップのレコードをとり計算時間を評価した．自然科学研究機構計算科学研究センターの計算機を使用した．詳細は以下の通りである．

Fujitsu PRIMERGY RX300 S7
CPU: Intel Xeon E5-2690
Memory: 8GB
Compiler: インテル (R) 64 対応インテル (R) C コンパイラー XE (インテル (R) 64 対応アプリ
ケーション用)、バージョン 15.0.2.164 ビルド 20150121 (C) 1985-2015 Intel Corporation.
また，本論文の結果は特に断らない限り，長さ σ ，質量 m，時間 τ =–7–

mσ 2 /kB T ，エネルギー ，温度

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）

/kB として無次元化してあらわす．時間刻み Δt は，0.005τ として計算を実行した．

–8–

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）

第2章

ベルレの近接リスト法
2.1 概要
カットオフ半径 rcut により，それぞれの粒子は半径 r ≤ rcut に存在する粒子としか相互作用をしない．
通常の固体や液体であれば，その数は
rcut

4πρ

drr2 g(r) ∼043πρrcut3(2.1)

で評価される．ここで密度 ρ = 1.0 の系を考える．LJ ポテンシャルで典型的に用いるカットオフ半径

rcut = 2.5 とすれば相互作用する粒子数は高々 200 個程度であり，系全体の総粒子数 (1000 粒子以上と考える) と比べると極めて少ない．そこで，各粒子についてあらかじめ相互作用をする粒子ペア，あるいは
粒子リストを設けておけば計算コストを抑えることができると期待される．実際には図 2.1 で表される
ように，カットオフ半径の少し外側に SKIN = rlist − rcut を設けておき，その中に含まれる粒子のインデックスを配列に記録しておく．この手法をベルレの近接リスト法 (Verlet Neighbor List: VNL 法) と呼ぶ [10, 11]．もちろんそのリストの作成の計算コストは O(N 2 ) であるが，その更新頻度を 10 ステップに 1 回，あるいは 100 ステップに 1 回などのように少なくできれば良い．実際には時間発展により粒子

図 2.1

ベルレの近接リスト法 (VNL) の概念図：SKIN = rlist − rcut ，rcut はカットオフ半径，rlist

はリスト半径を表す．灰色粒子は黒色粒子の近接リストに記憶されている粒子であり，毎ステップリ
ストを用いて相互作用を計算することが出来る．白色粒子はリストの外にでており，計算の必要の無い粒子である．

–9–《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）

は次第に移動していくので，ある程度計算が進んだ後に，そのリストを更新する必要がある．つまり粒子がリストの範囲から出てしまう前にリスト更新を行わなければならない．リストの更新を除けば，相互作用の計算コストは O(N 2 ) から O(N ) に低減することが可能である．固体の様に粒子がほとんど移動しな
いような系に対してはリスト更新の頻度を少なくする事が可能であり VNL 法は非常に効果的である．逆に言えば高温で粒子が激しく運動しているような液体系に対しては頻繁にリストを更新しなければならず，あまり効果的ではない．

VNL 法における近接リストの作成に関する該当部分を Algorithm1 に示す．まず始めにリストを構成する時点の粒子の配置を保存する．これは後にリストの更新のタイミングを決定するために必要となる．具体的には，毎ステップ，各粒子の変位を計算しそれが SKIN/2 を超えてしまうとリストを更新するということを繰り返す (詳細に関しては後に説明する)．次に，カットオフ半径 rcut よりも少し大きい半径 rlist = rcut + SKIN を設けておき，その中に含まれる粒子をリストに保存していく．point[i] は i 番目の粒子と相互作用する j 粒子のインデックスがリスト配列の何番目から読み取ればよいかを示している．
したがって，i 番目の粒子の VNL に属している相手方 j 粒子は list[i] の i = point[i] から point[i+1] - 1
までに収められていることになる．そして，リストを使用して力を計算する際には Algorithm2 の様にすれば良い．最後に相手方 j 粒子のインデックスを配列 list[ ] から参照し相互作用を計算する．
リストの更新間隔は 10 ステップなどに固定しても，次の更新までに粒子がリストから飛び出した場合には相互作用する粒子の誤参照が起こってしまい計算が破綻する．そこで毎ステップ，リストを構築した
時点からの粒子の変位を計算し，それがある閾値を超えてしまったらリストを更新すれば，これにより相互作用する粒子の誤参照を防ぐ事ができ，かつその更新間隔を最大限に取る事が可能となる．つまり，リストの更新間隔を事前に固定する必要がなくなる．この行程をリストの動的更新と呼ぶ [12–14]．リスト
動的更新の構造は Algorithm3 に示した．
ここで，閾値に関して考察する．図 2.2(a) に記すように黒い 2 つの粒子の運動に注目する．これら 2つの粒子は互いの近接リストに含まれ合っている．ここで，これら 2 つの粒子が互いに対角方向に速度 v

Algorithm 1 Make the Verlet Neighbor List
1: =====Save current conﬁguration.=====
2: for i=0 to NumParticle do
3:
rx0[i] ← rx[i]4:ry0[i] ← ry[i]5:
rz0[i] ← rz[i]6: end for
7: =====Make the Verlet Neighbor List.=====
8: (int)nlist ← 0
9: for i=0 to NumParticle - 1 do
10:
point[i] ← nlist + 1
11:
for j=i+1 to NumParticle do
12:
calculate rij
13:
if rij < rlist then14:
nlist ← nlist + 1
15:
list[nlist] ← j
16:
if rij < rcut then17:
calculate force
18:
end if19:
end if
20:
end for
21: end for
22: point[N-1] ← nlist + 1– 10 –

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）

Algorithm 2 Use the Verlet Neighbor List
1: for i=0 to NumParticle - 1 do
2:jbegin ← point[i]3:jend ← point[i+1] - 1
4:nlist ← point[i]5:if jbegin <= jend then6:for jnab = jbegin to jend do
7:j ← list[jnab ]8:calculate rij
9:
if rij < rcut then10:
calculate force
11:
end if
12:
end for13:
end if
14: end for

で等速直線運動したと仮定すると，n ステップ後の位置は右図の様になる．図 2.2(b) では各粒子は互いの近接リストの範囲から出てしまっている．この時各粒子がそれぞれ移動した距離は

v × ndt = SKIN/2(2.2)と計算する事が出来る．つまり，粒子変位が SKIN/2 になったときリストを更新すれば良い．この仮定は実際の多粒子系の MD 計算においてはほぼ起こりえない．なぜなら粒子は周囲の粒子と互いに衝突し合い，対角方向へ等速直線運動など起こりえないからである．したがって，粒子の変位が閾値 SKIN/2 を超えてしまったらリストを更新するという設定が妥当である事を意味している．
以上の事から，VNL 法の構造は Algorithm4 のようになる．毎ステップ，リストを更新するか否かチェックし，更新するならばリストを再構築すると共に力を計算する．更新しないならリストを用いて力を計算する事ができる．

(b) t = t0 + ndt における粒子配置．

(a) t = t0 における粒子配置．

Verlet Neighbor List criterion：(a)t = t0 において，黒色粒子はそれぞれ近接リストに粒子を記録してあり，互いにリストに含まれ合っている．(b)n ステップ経過した後，黒色粒子は互いの近図 2.2

接リストから飛び出す．

– 11 –

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）

Algorithm 3 Check update of the Verlet Neighbor List
1: dispmx ← 0
2: for i=0 to NumParticle do3:
dispmx ← max ( | rx[i] - rx0[i] | , dispmx)
4:
dispmx ← max ( | ry[i] - ry0[i] | , dispmx)
5:
dispmx ← max√ ( | rz[i] - rz0[i] | , dispmx)
6:
dispmx ← 3.0 dispmx
7: end for
8: if dispmx > SKIN / 2 then9:Make the Verlet Neighbor List
10: end if

Algorithm 4 the Verlet Neighbor List
1: ==== MD loop start =====
2: for step = 0 to nstep do
3:
move r(t) → r(t+dt) and v(t) → v(t+dt/2)
4:
Check update of the VNL
5:
if Need update of the VNL then6:
Make the VNL and calculate force O(N 2 )
7:
else
8:
Use the VNL and calculate force O(N )9:
end if10:
move v(t+dt/2) → v(t+dt)
11: end for

2.2 SKIN(リスト半径) の最適化
式 (2.2) より，リスト更新の条件はmax |r(t0 + ndt) − r(t0 )| >
i

SKIN
2

(2.3)

と書き換えられる．ここで n は更新間隔，dt はシミュレーションの時間刻みを示している．t = t0 の粒子位置からの変位を全粒子に関して毎ステップ計算しておき，そのうちのどれか１粒子でも変位がSKIN/2 を越えてしまうとリストを更新することになる．前述のように，更新間隔 n は粒子の運動に従い，つまり動的に変化する．よって，更新間隔 n は温度 T ，密度 ρ，粒子数 N ，カットオフ半径 rcut に依存することになる．
いま n + 1 ステップの計算を実行することを考える．このとき，n ステップの間はリストを使用 (O(N ))し，n + 1 ステップ目にリストを更新する為に全粒子探索 (O(N 2 )) をするので，その間，相互作用の計算に要する時間は解析的に



N −1
n
1
+
4πρ
TVNL = N τf
2
n+1
n+1

rlist


drr g(r) − 1
2

(2.4)

0

と，あらわすことができる．ここで τf は各粒子ペアの計算に要する単位時間量である．第 1 項はリスト更新に要する時間，第 2 項はリスト使用に要する時間を表している．

VNL 法に対してパラメータ (リスト半径) を最適化するために，リスト半径 rlist をシステムを設定しているパラメータ (N, T, ρ, rcut ) から見積もる近似を考える．式 (2.3) から，SKIN は更新間隔 n と粒子変位から見積もる事ができるので，SKIN = 2nx– 12 –(2.5)

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）とすればよい．ここで，x は最も長い距離移動し SKIN/2 を超える粒子の最大変位をあらわしている．
ここで x をどのように評価するかが重要となる．まず始めに考えられることは，拡散係数 D から長さスケール lD を見積もることである．つまり，x ∼ lD とし，ここで lD =√Dδt という近似である．しかしながらこの lD は，拡散過程から特徴付けられる長さスケールであり，リスト半径 rlist = 2.5 ∼ 5.0 に比べると非常に長く不適切である．そこで，短い長さスケールを特徴付けるための平均自由行程 lth を考える．系の温度 T から，エネルギー等分配則 12 mv 2 = 32 kB T (3 次元) を用いて平均速度 v を評価する
事で見積もる事ができる．これにより，より適切な近似 x ∼ lth が可能となり，lth = vδt とすればよい．しかしながらまだ十分ではない．平均速度 v から粒子の “平均変位” を求めるのではなく，全粒子の中から最も長い距離移動した粒子の “最大変位” を求めなればならない．そこで粒子の “最大変位” を

3kB T
m

x ∼ cδt(2.6)

とし，定数 c は Maxwell 分布の中で最も速い粒子を抽出するため c = 3.0 とする．図 2.3 は粒子数 1,000，

10,000，100,000 の３つのシステムに対する速度分布の結果を示している．
次により簡単化のため，粒子の分布を平均密度で近似する．つまり，g(r) = 1 とし，rlist

drr2 g(r) ∼

0

3
rlist
3

(2.7)

となる近似である．これにより，1 粒子が n + 1 ステップに要する計算時間は



4π 3
1
n
1
TVNL
=
+
ρr − 1fVNL =
2τf N (N − 1)
n + 1 n + 1 N − 1 3 list


4π
1
n
1
3
=
+
ρ(rcut + 2nx) − 1
n+1 n+1N −1 3

(2.8)

で評価する事が可能となる [15]．右辺第 1 項は n + 1 ステップ目に一度リストを更新する事を反映しており，第 2 項は n ステップの間，半径 rlist の球形の VNL を使用する事に対応している．リストの動的更新の手法を採用することで，式 (2.2) から明らかなように更新間隔 n はリスト半径 rlist あるいは SKIN に
10

N = 1,000
N = 10,000
N = 100,000

9
8
7

p(vx)

6
5
4
3
2
1
0

図 2.3

-4

-2

0

vx

2

4

粒子速度の x 成分の分布：T = 0.772，ρ = 0.8 として N = 1, 000, 10, 000, 100, 000 の 10

ステップ間の平均速度分布の結果をあらわしている. 緑線はエネルギー等分配則から見積もった速度

3kB T /m，紫線は 3



3kB T /m をあらわす．

– 13 –

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）

依存する．SKIN が大きければ，それだけ更新間隔を長く取る事ができるが，一方で，そもそも SKIN が大きい事でその中に含まれる粒子数
rlist

4πρ

4
3
πρrlist3

drr2 g(r) ∼

0

(2.9)

も増えてしまう．したがって，O(N 2 ) のリスト更新はほとんどしなくてもよくなるが，リストを使用して毎ステップ相互作用を計算する際に参照する粒子数が多くなるので，O(N ) (正確には O(N ×
リストに含まれている隣接粒子)) のリスト使用に計算時間がかかってしまう．逆に SKIN が小さければ頻繁にリスト更新を行ってしまい，O(N 2 ) の計算コストが大半を占めてしまう．極端な例である SKIN = 0

time [sec / N]

100

N = 1,000
N = 10,000
N = 100,000

10

1

0.1

0

1

2

3
4
SKIN = rlist - rcut

56(a) システムサイズ (粒子数) 依存性
1

1

T = 0.5
T = 1.0
T = 2.0

time [sec / N]time [sec / N]

ρ = 0.5
ρ = 0.8
ρ = 1.0

0.1 0

0.5

1
SKIN

(b) 密度依存性
図 2.4

1.5

2

0.1 0

time [sec / N]

1

0.5

1
SKIN

(c) 温度依存性

1.5

2

0.1 0

0.5

rcut = 2.5
rcut = 3.5
rcut = 4.5
1
1.5
SKIN

(d) カットオフ半径依存性

計算時間の SKIN サイズ依存性：点はシミュレーションの結果を，実線は式 (2.8) をそれぞれ表している．(a) システムサイズ依存性，(b) 密度依存性，(c) 温度依存性，(d) カットオフ半径依存性を検証した．計算時間を最小化する最適リスト半径 (SKIN サイズ) はシステムサイズ，密度，温度，
カットオフ半径に依存する．SKIN サイズを小さくとるとリスト更新 (O(N 2 )) の頻度が増え，計算コ
ストを要する．SKIN サイズを大きくすると参照体積が増え，隣接粒子数のループに計算コストを要する．

– 14 –

2

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）の時は毎ステップ毎にリストを更新する事になる．また，SKIN → L の極限ではシステム全体の粒子がリストに含まれることになり，リストを使用したとしても O(N 2 ) の計算コストを要してしまう．このよ
うに，SKIN と更新頻度は互いにトレードオフの関係にあり，計算時間を最適とする SKIN の特定が必要とされる．
式 (2.8) の結果とシミュレーションによる結果を図 2.4(a)∼(d) に示す．シミュレーションは LJ ポテン
シャルで NVE アンサンブルを 100,000 ステップ行った．縦軸は 1 粒子あたりに要する計算時間，横軸はSKIN である．各パラメータに対する依存性は以下のようにまとめられる．
(a) 粒子数依存性
温度 T = 0.772，密度 ρ = 0.8，カットオフ半径 rcut = 2.5 とし，システムサイズ (粒子数)N =1, 000, 10, 000, 100, 000 の 3 つの系に対しシミュレーションを行った．シミュレーションの結果，計算時間を最適化する SKIN サイズはそれぞれのシステムサイズに対して，SKIN = 0.5σ(N =1, 000), 1.1σ(N = 10, 000), 6.3σ(N = 100, 000) である．n + 1 ステップに一度行う更新のステッ
プでは VNL 法は O(N 2 ) の計算コストを伴う．更新までは O(N ) を維持することができるため，粒子数が増えるとその効果は顕著になる．従って，粒子数が増えるに伴い，なるべくリスト更新をしない方が (更新間隔 n が大きい方が) 計算コストを落とすことができると考えられる．つまり，粒子数の増加に対して更新間隔 n を延ばすために，最適 SKIN は大きくなる．

(b) 密度依存性
粒子数 N = 1, 000，温度 T = 0.772，カットオフ半径 rcut = 2.5 とし，密度 ρ = 0.5, 0.8, 1.0の 3 つの系に対しシミュレーションを行った．計算時間を最適とする SKIN サイズはそれぞれ，SKIN = 0.8σ(ρ = 0.5), 0.6σ(ρ = 0.8), 0.5σ(ρ = 1.0) である．密度はリストに含まれる粒子数に関係する. 当然，低密度であればリストに記憶する相互作用ペアの数は小さくなり，計算コストも下がる. また，低密度になるにつれ粒子が移動する平均自由行程が長くなる．これにより，その分だけ粒子はリストの範囲から出やすくなり，頻繁に O(N 2 ) の更新を行ってしまう．したがって，更新間隔を長くした方が計算効率は向上する．そのため，低密度なほど最適 SKIN は大きい．

(c) 温度依存性
粒子数 N = 1, 000，密度 ρ = 0.8，カットオフ半径 rcut = 2.5 とし，温度 T = 0.500, 1.000, 2.000の 3 つの系に対しシミュレーションを行った．計算時間を最適とする SKIN サイズはそれぞれ，SKIN = 0.4σ(T = 0.5), 0.5σ(T = 1.0), 0.7σ(T = 2.0) である．温度は粒子の速度に依存する．シ
ステムの温度が高いほど粒子の速度は増加し，それに伴い粒子がリストの範囲から出る頻度も増加
する．そのため温度が高い系では，なるべく更新間隔 n を大きくし，リストをより長く使い続けた方が計算効率が良くなるため，SKIN を大きく取った方が良い．

(d) カットオフ半径依存性
粒子数 N = 1, 000，温度 T = 1.000，密度 ρ = 0.8 とし，カットオフ半径 rcut = 2.5, 3.5, 4.5の 3 つの系に対しシミュレーションを行った．計算時間を最適とする SKIN サイズはそれぞれ，SKIN = 0.6σ(rcut = 2.5), 0.4σ(rcut = 3.5), 0.3σ(rcut = 4.5) である．カットオフ半径が大きいほど，カットオフ球体積内に含まれる相互作用を計算するペアの数は多くなる．つまり，リストに含まれる粒子ペアも多くなる ( O(N × リストに含まれている隣接粒子数) )．極端な例をあげれば，
カットオフ半径がシステムの一辺の長さと等しい場合，リスト使用による計算コストが O(N 2 ) となってしまう．つまり，カットオフ半径が大きくなるにつれ，リスト更新による計算コストの寄与は相対的には小さくなる．逆に，カットオフ半径が小さければシステム全体に比べてリスト体積は– 15 –

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）

非常に小さいので，なるべくリスト更新をしないほうが効率が良い．そのため，カットオフ半径が小さいほど最適 SKIN は大きい．
図 2.4(a)∼(d) において，SKIN が大きくなるにつれ，式 (2.8) の結果とシミュレーションの結果に差が生じていることがわかる．これは式 (2.6) において c = 3.0 と仮定したことに起因する．この c = 3.0 は系に含まれている粒子の最大変位を見積もるための指標として採用した値であった．粒子数が増加するに伴い，あるステップにおける Maxwell 分布において 3

3kB T /m を上回る速度を持つ粒子が存在する確率は，粒子数が多い方が大きい (例えば 100,000 粒子のうち 3

3kB T /m を上回る速度をもつ粒子を見出す確率が 1/100,000 であるとすると，100,000 粒子系では毎ステップ 10,000 粒子系では 10 ステップに 1
度，そのような速度を持つ粒子を観測することになる)．c が小さければ最大速度をもつ粒子が実際よりも高い確率で観測されると仮定する事になる．つまり，実際よりも更新頻度が上がると仮定してしまうことになり，したがって，SKIN が大きい領域で評価式はシミュレーション結果より大きな値をとる．しかしながら，式 (2.8) において最も興味があるのは，計算時間を最小化する最適な SKIN の特定であるから，

c = 3.0 としても N ∈ [1000, 100, 000] のシステムサイズに対しては十分に最適値を特定する事ができる．
なお，VNL 法においてはリスト更新の際の全粒子探索の計算コスト O(N 2 ) がボトルネックとなるため，巨大なシステムサイズ (N > 100, 000) に対しては用いられることはない．

– 16 –

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）

第3章

セル分割法
3.1 概要高速化のためのもうひとつの有益な方法として，計算系を一辺が Lcell の小さな立方体セルに分割する事を考える．長さ Lcell を rcut ≤ Lcell を満たすように選ぶことで，粒子間の相互作用は同じセルに属す
る粒子間，もしくは隣接セルに属する粒子間に限定される．図 3.1 に示すように，隣接するセルの数は限られている (2 次元なら 8 つ，3 次元なら 26 個) ことから，大きな系であれば相互作用の数を大幅に減らすことができる．この方法をセル分割法 (Cell Linked List: CLL 法) と呼ぶ [10, 11]．それぞれのセルに属する粒子の数は有限 (高々 100 粒子程度) であるから，相互作用に要する計算時間は O(N ) に抑える事ができる．分割するシステムが直方体であれば一般に分割セルも直方体となるが，それぞれのセルの一辺の長さが rcut よりも長くすればよい．

3 次元系における CLL 法では，参照すべき隣接セルを特定するには Algorithm5 のようにおこなう．
隣接セルの特定作業をマッピングと呼ぶ．ここで用いている関数 icell(ix, iy, iz) の詳細は Algorithm6に示している． ここで 3 × 3 のセルの中心のセルインデックスを imap とする (図 3.2)．さらにこ

図 3.1

セル分割法 (CLL) の概念図：Lcell はセルの一辺の長さ，rcut はカットオフ半径をそれぞれ表している．黒色粒子は中央の灰色セルに属しており，相互作用の計算の為に参照する粒子は隣接している白色セルに含まれている．

– 17 –

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）

Algorithm 5 Create the maps of the cell list (Cell Index)
1: int ix, iy, iz
2: int imap, mapsize, maps[mapsize]
3: int M (システムボックス一辺の分割数)
4: mapsize = 13 * M * M * M (13 は周囲の参照セル数，M*M*M はシステムの全セル数)
5: ===== Initialize the array of the map =====
6: for imap = 0 to mapsize do
7:
map[imap] ← 0
8: end for
9: ===== create the maps of the cell list =====
10: for iz=0 to M do
11:
for iy=0 to M do
12:
for ix=0 to M do
13:
imap ← ( icellNo(ix, iy, iz) ) * 13;14:map[ imap + 1 ] ← icellNo( ix+1, iy , iz );15:map[ imap + 2 ] ← icellNo( ix+1, iy+1, iz );16:
map[ imap + 3 ] ← icellNo( ix , iy+1, iz );17:
map[ imap + 4 ] ← icellNo( ix - 1, iy+1, iz );18:
map[ imap + 5 ] ← icellNo( ix+1, iy , iz - 1);
19:map[ imap + 6 ] ← icellNo( ix+1, iy+1, iz - 1);
20:map[ imap + 7 ] ← icellNo( ix , iy+1, iz - 1);
21:map[ imap + 8 ] ← icellNo( ix - 1, iy+1, iz - 1);
22:map[ imap + 9 ] ← icellNo( ix+1, iy , iz+1);23:map[ imap + 10 ] ← icellNo( ix+1, iy+1, iz+1);24:map[ imap + 11 ] ← icellNo( ix , iy+1, iz+1);25:map[ imap + 12 ] ← icellNo( ix - 1, iy+1, iz+1);26:map[ imap + 13 ] ← icellNo( ix , iy , iz+1);
27:
end for28:
end for
29: end for

Algorithm 6 function icellNo(ix, iy, iz): Return the Index of the Cell
1: int ix, iy, iz
2: int icellNo
3: int M (システムボックス一辺の分割数)
4: icellNo ← mod(ix + M, M)+ mod(iy + M, M) * M+ mod(iz + M, M) * M * M(mod(a, b) は a を b で割った余りを返す)

5: return(icellNo)の imap から周囲に隣接した参照セルのインデックスを imap+1 から imap+13 まで割り振る．周囲に隣接したセルは 26 個あるが，粒子間相互作用は作用・反作用の法則から図 3.2 に示すように，

Ncell = (3 × 3 × 3 − 1)/2 = 13 個のみ参照できればよい．ここで Ncell は参照セル数と呼ぶ．つまり，配列 map[] に隣接セルのインデックスを埋め込んで行く．このマッピングにはもちろん計算コストを要するが，空間を分割するだけであるから粒子がどのセルに属するかは考える必要はないため，MD を開始
する前にあらかじめ行っておけばよい．
次に，各ステップにおいて計算される粒子の位置とそれに対応したセルインデックス imap をリンク
させ，隣接したセルから粒子ペアを特定する (Algorithm7)．ここで便宜上，システムの一辺を [−L/2,

L/2] から [0, 1] に変数変換して i 粒子の位置を格納し，そこから対応するセルインデックス icell を決めていることに注意する．セルインデックス icell に含まれる i 粒子を配列 head[icell](はじめは 0 に初期化されている) から list[i] として保存する．その次に i 粒子のインデックスを head[icell] に上書き保存して行く．これを繰り返すことで CLL 法におけるリスト構造を実現する．すでに述べたように，この作業に– 18 –

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）

図 3.2

隣接セルのマッピング：中央のセル imap が参照する隣接セルのインデックスを imap+1 から imap+13 まで割り振る．作用反作用の法則から参照セルは対角方向のみを割り当てれば良い．

図 3.3

左：セルインデックス icell

図 3.4

head 配列，list 配列のリンク構造の概念
図：position2 のセルには head[2] = 8，つまり 8番の粒子が先頭に入っている．さらに list[8] = 7，
つまり 8 番粒子の次は 7 番粒子が入っている．list
配列が-1 を返すまでこれを繰り返す．

右：粒子インデックス i のリンクを表している．配列 head[icell]を用いてセルに含まれている粒子インデックス i の先頭を参照し，粒子インデックス i から配列 list[i]を用いて次に内包している粒子インデックスを参照していく．

は各ステップ O(N ) の計算コストを要する．リンク作業の具体的な様子を図 3.3，3.4 に示す．第 1 に，
セルインデックスから配列 head[ ] を用いて，セルに含まれる最も大きいインデックスを持つ粒子を参照
する．第 2 に次に小さいインデックスを持つ粒子を配列 list[ ] を用いて参照する．これを繰り返すことにより，最終的に-1 を終了判定として次のセルへとこの作業を移行する．例えば，図 3.4 では，セルイン
デックス 2 から head[2] = 8，つまり 2 のセルに含まれる最も大きなインデックスを持つ粒子は 8 だとわ
かる．ここから，list[8] = 7 であることから，8 の粒子の次は 7 を参照し，その次は 4，1，-1 といった具合である．リスト参照の具体的なアルゴリズムを Algorithm8 に示している．

Algorithm5 から Algorithm8 に基づき，全体のアルゴリズムをまとめると Algorithm9 のようになる．つまり各ステップで相互作用を計算する前に隣接セルのマッピングを行い粒子の位置とセルの住所をリンクしセル参照を行う．

– 19 –

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）

Algorithm 7 Link the Particle Index to the Cell Index
1: int icell, i, j
2: int M (システムボックス一辺の分割数)
3: double Systemsize (システムの一辺の長さ)
4: double cellsize (セルの一辺の長さ)
5: cellsize ← Systemsize / (double)M
6: if cellsize < rcut then7:

Error: Cell size is too small for cutoﬀ radius.

8: end if
9: for i=0 to NumParticle do
10:dmyrxi ← rx[i] - ﬂoor( rx[i] + 0.5)

dmyryi ← ry[i] - ﬂoor( ry[i] + 0.5)
dmyrzi ← rz[i] - ﬂoor( rz[i] + 0.5)
13:
icell ← (int)( ( dmyrxi + 0.5 ) / (double)M )+ (int)( ( dmyryi + 0.5 ) / (double)M ) * M+ (int)( ( dmyrzi + 0.5 ) / (double)M ) * M * M
14:list[i] ← head[icell]15:
head[icell] ← i
16: end for
11:12:

3.2 排他的セル分割法
CLL 法は実装が比較的簡単なことから良く用いられる手法であり，システムサイズに関係なく計算効率を O(N ) で抑えることが出来る．しかしながら，少ない粒子数の系 (∼ 1, 000 粒子数以下) では VNL

(a)
図 3.5

(b)

(c)

排他的セル分割法: (a) 第 1 隣接セル参照法 (fc = 1, Ncell = 33 個)，(b) 第 2 隣接セル参照法(fc = 2, Ncell = 53 個)，(c) 第 3 隣接セル参照法 (fc = 3, Ncell = 73 個)．括弧の fc は x, y, z 方向の隣接参照セル数，Ncell は参照セル数をあらわしている．太枠で囲った領域が参照体積 Vref ，円で囲まれた領域がカットオフ体積 Vcutoﬀ である．参照セル数 Ncell は第 1 隣接セル参照なら 33 個あるが，粒子間相互作用は作用・反作用の法則から (33 − 1)/2 = 13 個，第 2 隣接セル参照なら (53 − 1)/2 = 62個，第 3 隣接セル参照なら (73 − 1)/2 = 171 個というように，対角方向に属する半分のセルのみを参照できればよい．

– 20 –

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）

Algorithm 8 Use the cell linked List and calculate force
1: for icell=0 to Ncell do2:====== select i particle contain in each cells ======
3:i ← head[icell]
4:
while i >= 0 do
5:====== j particle contain in the current cell ======
6:j ← list[i]7:while j >= 0 do
8:
Caluculate rij
9:
if rij < rcut then10:
Calculate force
11:
end if12:
end while13:
j ← list[j]14:====== j particle contain in the neighbor cells ======
15:
jcell0 ← reference cell × icell
16:
for neighbor = 0 to reference cell do
17:
jcell ← map[ jcell0 + neighbor ]18:
j ← head[ jcell ]19:while j >= 0 do
20:
Caluculate rij
21:
if rij < rcut then22:
Calculate force
23:
end if24:
j ← list[j]25:
end while26:
end for27:i ← list[i]28:
end while
29: end for

Algorithm 9 Cell Linked List
1: Create the map of the cell list.
2: ===== MD Loop Start =====
3: for step = 0 to nstep do
4:
move r(t) → r(t+dt) and v(t) → v(t+dt/2)
5:Link the particle index to the cell index.
6:
Use the cell linked list and calculate force.
7:
move v(t+dt/2) → v(t+dt)
8: end for法に比べ計算効率が劣る．本質的に計算に必要とする領域は，カットオフ半径 rcut 内に含まれる領域
3(以下，カットオフ体積を Vcutoﬀ = 4πrcut/3 とする) であり，それに比べて参照しなければならない体積 (以下，参照体積 Vref = 33 × L3cell ) が大きいことが要因として挙げられる．具体的に 3 次元系において Lcell = rcut = 2.5 とすれば，Vref /Vcutoﬀ は 6.4 倍である．つまり，カットオフ半径に含まれないが参照しなければならない領域が 85 ％も存在していることを意味している．少ない粒子数の系に対し
3て VNL 法の方が計算効率が優れているのは，VNL の参照体積は Vref = 4πrlist/3 であり，rcut = 2.5,

rlist = rcut + 0.5 に取ったとしてもカットオフ球の外側にある参照領域はおおよそ 42 ％である．粒子数が少なければ，VNL 法のボトルネックであったリスト更新の際の O(N 2 ) はそれほど全体の計算時間に対し甚大な寄与にはならない．したがって参照体積を小さく取ることが計算効率の向上に直結している．そ– 21 –

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）

Algorithm 10 Create the maps of the exclusive cell list (Cell Index)
1: int imap, mapsize, maps[mapsize], count2: int M (システムボックス一辺の分割数)
3: int reference cell (参照セル数), adjacent cell (各方向の隣接セル数)
4: refernce cell ← (2.0∗adjacent cell+1)3 − 1 / 2.0
5: =====Initialize the array of the map=====
6: for imap = 0 to mapsize do
7:
map[imap] ← 0
8: end for
9: =====Create the array of the map=====
10: for iz=0 to M do
11:
for iy=0 to M do
12:
for ix=0 to M do
13:
imap ← icell(ix, iy, iz) * reference cell
14:
count ← 0
15:
for jx= -adjacent cell to adjacent cell do
16:
for jy= -adjacent cell to adjacent cell do
17:
for jz= -adjacent cell to adjacent cell do
18:
if jy > 0 then19:
count ++
20:
map[ imap + count] ← icell(ix + jx, iy + jy, iz + jz)
21:else if jy == 0 then22:
if jx >= 0 AND jx > 0 then23:
count ++
24:
map[ imap + count] ← icell(ix + jx, iy + jy, iz + jz)
25:else if jx >= 0 AND jz == 0 then26:
count ++
27:
map[ imap + count] ← icell(ix + jx, iy + jy, iz + jz)
28:else if jx == 0 AND jz < 0 then29:
count ++
30:
map[ imap + count] ← icell(ix + jx, iy + jy, iz + jz)
31:
else
32:
...nothing to do
33:
end if34:
end if
35:
end for36:
end for37:
end for38:=====End of the adjacent cell loop=====
39:
end for40:
end for
41: end for
42: =====End of the M loop=====

こで，セルをより細分化し，さらに外側のセルまで参照することを考える．つまり，セルサイズの条件を

Lcell >rcut
rcut
rcut
> Lcell >fc − 1fc(fc = 1)

(3.1)

(fc > 1)

(3.2)

とすることで参照体積を小さくする事が可能となる (図 3.5)．ここで，fc は x, y, z 方向の隣接参照セル数をあらわしている．fc = 1 の場合は，前節に述べた従来型の CLL 法である．このようにセルの細分化を行う方法を排他的セル分割 (Exclusive cell linked list) 法と呼ぶ [16, 17]．図 3.5 に排他的セル分割法の概念図を示した．セルをより細分化することで 1 つのセルに 1 つ以上の粒子が入ることができなくなることに注意する．

– 22 –

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）

前節で説明され，また図 3.5(a) で示されている CLL 法は，第 1 隣接のセルまで参照する方法 (以下，第一隣接セル参照) である．この方法ではセルの一辺の長さ (セルサイズ Lcell ) はカットオフ半径 rcut よりも大きいものとしていたが，第 2 隣接セル法では注目している粒子が属しているセルから 2 つ外側まで
セル参照を行う．この場合セルサイズは rcut > Lcell > rcut /2 を満たさなければならない図 3.5(b)．セルサイズはLcell =

L
M

(3.3)

のようにシステム一辺を均等に分割して設定する．ここで一辺のセル分割数は M である．例えば，

N = 1, 000，ρ = 0.8，T = 1.0，rcut = 2.5 の系においてはシステム一辺の長さは L ∼ 10.772 であり，
• Lcell /rcut ∈ [1.0772(M = 4), 1.436(M = 3)]：第 1 隣接セル参照 (fc = 1, 図 3.5(a))• Lcell /rcut ∈ [0.478(M = 9), 1.0772(M = 4)] ：第 2 隣接セル参照 (fc = 2, 図 3.5(b))• Lcell /rcut ∈ [0.331(M = 13), 0.478(M = 9)] ：第 3 隣接セル参照 (fc = 3, 図 3.5(c))のようにセル分割を行うことができる．それぞれ，セル分割はシステムの一辺の長さを等分割しなければいけないことからセルサイズは離散値を取る．分割数 M = 4 とすれば

Lcell =

10.772L=∼ 2.693M
4(3.4)

Lcell =

10.772L=∼ 2.154M5(3.5)となり，さらに M = 5 とすればとなるが，従来型の CLL 法 (fc = 1) では Lcell ≥ rcut = 2.5 でなければならないので 5 分割はできない．
しかし，第 2 隣接セル参照 (fc = 2) の場合，2 つ隣りまで参照数を増やすことで rcut > Lcell ≥ rcut /2 のようにセルサイズを設定することが出来るようになる．第 2 隣接セル参照 (fc = 2) を取れば 5 分割が可能になる．Algorithm10 でアルゴリズムとして示しているように隣接セルのマッピングを行うことで，
参照セル数を増やしセルサイズを小さくすることが可能となる．

3.3 セルサイズの最適化
排他的セル分割法を用いることで，従来型 CLL 法からさらなる計算効率の向上が期待される．セルを細かくすることで参照体積が小さくなり，計算に関係のない参照体積は減少する．しかし，その一方で参照しなければならないセル数が 33 , 53 , 73 · · · と増加する．これによりセルインデックスのループに計算時間を要するようになる．従って，計算時間に対するセルサイズには最適値が存在することになる．そこ
でセル分割法に関しても，VNL 法と同様の解析的評価式を提案する．セル分割法のアルゴリズムをまとめると，

1. システムをセルで領域分割し，隣接セルのマッピング (Algorithm9)
2. 各粒子がそれぞれどのセルに属するかリンクしリストを作成 (Algorithm7)
3. セルに関するループ (Algorithm8)
4. その中に含まれる粒子で距離の計算 (Algorithm8)
5. さらにカットオフ半径内に含まれる粒子のみ相互作用の計算 (Algorithm8)の 5 つのフェーズから構成される．この中で計算時間がセル数に依存するのは 1, 3, 4, 5 のフェーズである．フェーズ 2 は粒子のインデックスに関してループを行うだけであるから，計算時間はセル数には依存– 23 –

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）

しない．さらに，フェーズ 1 は相互作用計算を開始する直前に 1 回だけ行うため，全体の計算時間にはほとんど寄与しない．従って，計算時間のセルサイズ依存性が顕著にあらわれてくるのは 3, 4, 5 のフェーズとなる．以上の考察から，1 粒子が各ステップに要する計算時間の解析的評価式を

1
TCLL = (Ncell + 1) × τα2
1 3
ρLcell (Ncell + 1) − 1 × τβ+23
4 πrcut+× τγ
3 Ncell L3cell(3.6)のようにあらわすことができる [15]．右辺第 1 項はセルインデックスのループに要する時間，第 2 項はセル内に含まれる粒子インデックスのループに要する時間，第 3 項はカットオフ半径体積に含まれる粒子ペアに対して相互作用を計算する時間に対応している．ここで，τα ,τβ ,τγ は計算機に依存し，あらかじめ 1
ステップでの，各フェーズに要する時間を計測することによって与えられる．式 (3.6) から明らかなように，CLL 法の 1 粒子あたりの計算時間は密度，カットオフ半径に依存するが温度には依存しない．つまり各ステップにおいて粒子の “位置” とセルのインデックスとをリンクし，隣り合ったセルから粒子のイ
ンデックスを参照してくるだけの手法であるから，この手順には一切粒子の “速度” に依存する要素は含
まれない．言い換えれば，CLL 法はシステムの動的要素に依存しないリストアルゴリズムである．
式 (3.6) による計算時間のセルサイズ依存性と，シミュレーションの結果を図 3.6(a),(b) に示す．VNL
法と同様，LJ ポテンシャルで NVE アンサンブルを 100,000 ステップ計算し，粒子数あたりの計算時間として評価している．例えば図 3.6(a) の赤線にあらわされるように，この系では第 2 隣接セル参照の手法で Lcell /rcut ∼ 0.478(9 分割) のセルサイズが最も効率が良い．各パラメータに対する依存性は以下のようにまとめられる．

(a) 密度依存性システムの体積 V は密度 ρ，粒子数 N によって決まり，システムの一辺の長さは L = (N/ρ)1/3
で与えられる．このことからセルはシステムの一辺を均等に分割して作らなければならないので，取りうるセルサイズ Lcell はそれぞれの密度に応じて変わってくる．また，当然のことながら，高密度なほどカットオフ球体積に含まれる粒子数が増えるので計算コストがかさむ．これは式 (3.6)の第 2 項の寄与であるが，ただし，セル内に含まれる粒子ペア数のカウントに要する時間に対応しているので，カットオフ半径 rcut = 2.5 ではどの密度においても第 2 隣接セル参照で最も細分化したセルサイズが最適であることがわかる．

(b) カットオフ半径依存性
カットオフ半径 rcut = 4.5 では第 3 隣接セル参照 (fc = 3) が最適であり，rcut = 3.5, 2.5 では第2 隣接セル参照 (fc = 2) が最適，rcut = 1.5 では第 1 隣接セル参照 (fc = 1) が最適であることがわかる．これは式 (3.6) の第 3 項に着目すると理解できる．カットオフ半径 rcut が大きければ，対応して大きなセルでシステムを分割しなければならない．つまり，セル分割には式 (3.1) の条件を満たさなければならない．粒子数 N = 1, 000，密度 ρ = 0.8 の系では一辺が Lcell = 10.772 となり，一辺を 3 等分することを考えると，セルサイズは Lcell = L/M = 10.772/3 = 3.5906 となり，カットオフ半径が rcut = 4.5 であったとすると上記の条件を満たさない．つまり，カットオフ半径

rcut = 4.5 の系ではそもそも第 1 隣接セル参照はできない．これらのことから，カットオフ半径が大きくなるにしたがい，セル分割数，隣接参照セル数を増やした方が計算効率が良くなる．

– 24 –

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）

100
10
time [sec / N]time [sec / N]

10

10.2

0.4

0.6

ρ = 0.80
ρ = 0.60
ρ = 0.30
0.8
1

0.1
0.2

1.2

Lcell / rcut

0.4

0.6

rcut = 1.50
rcut = 2.50
rcut = 3.50
rcut = 4.50
0.8
1

1.2

Lcell / rcut(a) 密度依存性
図 3.61(b) カットオフ半径依存性

計算時間のセルサイズ依存性：点はシミュレーションの結果を，実線は式 (3.6) をそれぞれ表している．(a) 密度依存性，(b) カットオフ半径依存性を検証した．計算時間を最小化する最適リスト
半径 (SKIN) サイズはシステムサイズ，密度に依存する．セルサイズを小さいくとるとセル数が増え，
セルインデックスのループに計算コストを要する．セルサイズを大きくすると参照体積 (粒子数) が増
え，粒子数のループに計算コストを要する．図 (a) において，いずれの密度に対しても第 2 隣接セル参照法が最適である．図 (b) において，rcut = 1.5 では第 3 隣接セル参照 (fc = 3) が最適，rcut = 2.5，

rcut = 3.5 では第 2 隣接セル参照 (fc = 2) が最適，rcut = 4.5 では第 1 隣接セル参照 (fc = 1) が最適である．

– 25 –

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）

第4章

VNL-CLL 法第 2 章で紹介したベルレの近接リスト法 (Verlet Neighbor List: VNL 法)，第 3 章で紹介したセル分割法 (Cell Linked List: CLL 法) および，排他的セル分割法の長所・短所を以下にまとめる．

• VNL 法– リスト更新に O(N 2 ) の計算コストを要するため，10,000 粒子数以上の系に対してはそれがボトルネックとなり現実的でない．一方で，少ない粒子数の系に対しては高いパフォーマンスを発揮する．

– 最も速く移動する粒子によって更新間隔 n が決定されてしまう．これは気液共存系などの場合に適さないが，運動性の乏しい結晶固体系などには n は長くなるメリットがある．

– リスト半径 rlist を最適に選ぶためには試行錯誤が必要であったが，式 (2.8) によって決定することが可能．

– CLL 法に比べ参照体積が小さく，リスト上に保っておく粒子数も少なく (メモリ量が少ない)，
その分だけ計算速度は速い．

– 並列計算との相性が悪い．
• CLL 法– 全体として O(N ) の計算量であるため，10,000 粒子数以上の系に対しても高いパフォーマンスを実現することが可能．

– 粒子の運動性に影響しないため，気液共存系などの不均一系の場合に対してもデメリットはない．

– 排他的セル分割法では，最適化なセルサイズは自明ではなかったが，式 (3.6) によって決定することが可能．

– VNL 法と比べ参照体積は大きく，リスト上に保っておく粒子数は多い．さらにリスト配列も複数ある (メモリ量が多い)．その分だけ計算速度は遅い．

– 並列計算との相性が良い．
これらより，参照体積の点で VNL 法では計算効率が良いが，数ステップに一度リスト更新 O(N 2 ) を必要とするためそれがボトルネックとなっている．そこで，このリスト更新によるボトルネックを CLL法を用いて補うことが出来ればさらなる高速化が期待される．そこで VNL-CLL 法を提案する．

– 26 –

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）

4.1 概要VNL-CLL 法は VNL 法の動的更新を基本としている．つまり，更新するまでの間 n ステップはベルレの近接リストを用いて相互作用を計算する．つまり，各ステップに近接リストを更新するか否かチェッ
クを行い，動的リスト更新の手法を採用する．n + 1 ステップ目に全粒子探索 (O(N 2 )) を行うことで，

VNL を更新するのではなく，その代わりに CLL(O(N )) を用いて粒子探索を行う．これにより全体でO(N ) の計算を実現することが可能となる．VNL を基本的に用いており，さらに参照体積を最小化する
ことができることから高速化が見込める．VNL-CLL 法の詳細なアルゴリズムは Algorithm11 であり，Algorithm12, Algorithm13 から構成されている．

図 4.1

VNL-CLL 法の概念図：青色円は VNL 法における近接リスト半径をもった円，CLL 法において緑四角は参照する隣接セル領域を表している．セルの中に Verlet の近接リストを内包させるため，セルサイズの取りうる大きさの下限が rlist で決まっていることに注意する．

Algorithm 11 VNL-CLL
1: Create the map of the cell list.
2: Make the VNL and calculate force.
3: ===== MD loop start =====
4: for step = 0 to nstep do
5:
move r(t) → r(t+dt) and v(t) → v(t+dt/2)
6:
Check update of the VNL.
7:
if Need update of the VNL then8:=====Make the VNL using CLL=====
9:
Link the particle index to the cell index.
10:
Use the CLL and make the VNL.
11:
Calculate force.
12:
else
13:
=====Use the VNL=====
14:
Use the verlet neighbor list and calculate force.
15:
end if16:
move v(t+dt/2) → v(t+dt)
17: end for– 27 –

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）

Algorithm 12 VNL-CLL: Make the VNL using CLL
1: for icell = 0 to Ncell do2:=====select i particle contain in each cell=====
3:i ← head[icell]
4:
while i >= 0 do
5:=====j particle contain in the current cell=====
6:j ← cell list[i]7:while j >= 0 do
8:
calculate rij
9:
if rij < rlist then10:neighbor list[i][point[i]] ← j
11:point[i] ++12:
if rij < rcut then13:
calculate force
14:
end if15:
end if16:
j ← cell list[j]17:
end while18:=====j particle contain in the neighbor cell=====
19:
jcell0 ← reference cell * icell
20:
while neighbor = 1 to reference cell do
21:
jcell ← map[ jcell0 + neighbor ]22:
j ← head[jcell]23:
while j >= 0 do
24:
calculate rij
25:
if rij < rlist then26:
neighbor list[i][point[i]] ← j
27:point[i] ++28:
if rij < rcut then29:
calculate force
30:
end if31:
end if32:
j ← cell list[j]33:
end while34:
end while35:
i ← cell list[i]36:
end while
37: end for
38: for i=0 to NumParticle do39:
end point[i] ← point[i]40: end for

第 1 に，CLL 法と同様に相互作用を計算する前に隣接セルインデックスのマッピングを行う．第 2に，初期配置をもとに Verlet の近接リストを作成し相互作用を計算する．位置の更新，速度の更新をした後，VNL を更新するか否かチェックする．更新するのであれば CLL 法を用いて VNL を再構成する(Algorithm12)．更新しないのであれば，そのまま VNL を用いて力を計算する (Algorithm13)．ここで，VNL-CLL 法ではセルサイズの設定条件が CLL 法とは異なることに注意する．CLL 法ではカット
オフ半径より大きな一辺を持ったセルでシステムを均等に分割していたが，VNL-CLL 法では参照セルの内部にカットオフ半径 rcut だけでなく VNL 法のリスト半径 rlist も含まれる．したがって，セルサイズ

Lcell の取りうる条件は以下のようになる．
Lcell >rlist
rlist
rlist
> Lcell >fc − 1fc(rlist = rcut + SKIN)– 28 –(fc = 1)

(4.1)

(fc > 1)(4.2)

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）

Algorithm 13 VNL-CLL: Use the VNL
1: for i=0 to NumParticle do2:
jbegin ← 0
3:jend ← end point[i] - 1
4:
if jbegin <= jend then5:for jnab = jbegin to jend do
6:j ← neighbor list[i][jnab]
7:
if j > -1 then8:
calculate rij
9:
if rij < rcut then10:
calculate force
11:
end if12:
end if13:
end for14:
end if
15: end forとなる．セルサイズはリスト半径によって下限が決められていることに注意が必要である．

4.2 リスト半径・セルサイズの最適化
VNL 法，CLL 法の結果に基づいて VNL-CLL 法に対する 1 粒子が 1 ステップに要する計算時間に対する解析的評価式は


4
1
1 TCLLn3(τCLL /τVNL ) +πρ(rcut + SKIN) − 1
fVNL−CLL =n+1 N
n+1N −1 341n
1 TCLL CLL VNL
3(Vref /Vref ) +πρ(rcut + 2n x) − 1=
n+1 N
n+1N −1 3(4.3)のようになる．VNL 法を基にベルレの近接リストを更新する際に CLL 法を用いるのであるから，右辺第1 項のリスト更新に対応する項に CLL 法に要する時間 TCLL が含まれる．τCLL , τVNL はそれぞれ，1 ス
テップに CLL，VNL を使用し計算する時間であるが，ともに参照体積に比例する．実際には，半径 rlistの球である VNL の参照体積と，CLL 法での参照体積は一辺 2rlist の立方体との体積比であるから，

τCLL
V CLL
(2rlist )3
= ref
=
3 /3
VNL
τVNL
4πrlist
Vref

1.91

(4.4)

として考えれば良い．したがって，VNL-CLL 法の計算時間は温度 T ，密度 ρ，粒子数 N ，カットオフ半
径 rcut に依存する．計算時間の SKIN 依存性，セルサイズ依存性として図 4.2(b) に，シミュレーションの結果を，図 4.2(a) に式 (4.3) の結果を示す．カラーバーは粒子数あたりの計算時間を表しており，横軸
はセルサイズ，縦軸は SKIN である．セルサイズ Lcell が大きいと，その参照体積が大きく，計算コストが大きい．セルサイズが小さすぎるとセル数が増え計算コストがかかってしまう．また，セルサイズはシ
ステムの一辺を均等に分割しなければならないため離散値を取る．SKIN は大きすぎると，VNL の更新
頻度を低減することができるが，参照体積が大きく計算コストが大きい．また，SKIN の上限はセルサイズによって決まっている．SKIN が小さいと更新頻度が上がる．更新には CLL 法を用いるため O(N 2 )ではないにせよ，参照体積は CLL 法の方が大きいため計算コストを要する．したがって，計算コストを最適とする両者の組み合わせが存在する．各パラメータに対する依存性は以下のようにまとめられる．

– 29 –

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）

0.5
3.4

0.45
0.4

rcut + SKIN

3.2

0.35
3

0.3
0.25

2.8

0.2

2.6

0.15
0.25

0.5

1

2

Lcell / rcut(a) シミュレーションの結果．

0.5
3.4

0.45
0.4

rcut + SKIN

3.2

0.35

3

0.3

2.8

0.25
0.2

2.6

0.15
0.25

0.5

1

2

Lcell / rcut(b) 式 (4.3) の結果．
図 4.2

計算時間に対するリスト半径 (SKIN サイズ)，セルサイズ依存性. N = 1, 000，T = 1.0,

ρ = 0.5, rcut = 2.5 で 1 粒子あたりの計算時間を計測した．この系において最適リスト半径 rlist =rcut + SKIN，最適セルサイズの組み合わせはそれぞれ rlist = 3.1, Lcell /rcut = 0.42(fc = 3)((a) シ
ミュレーション)，rlist = 3.15, Lcell /rcut = 0.42(fc = 3)((b) 式 4.3) である．セルサイズが大きいと参照体積が増え，計算コストを要する．セルサイズが小さいとセル数が増え，セルインデックスのループに計算コストを要する．SKIN サイズが大きいと参照体積が増え，計算コストを要する．SKIN– 30 –

サイズが小さいと更新頻度が増え，リストの再構成に計算コストを要する．

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）(a) 粒子数依存性 (図 4.3)
温度 T = 1.0，密度 ρ = 0.8，カットオフ半径 rcut = 2.5 とし，システムサイズ (粒子数)N =500, 1, 000, 10, 000, 100, 000 の 4 つの系に対し計算時間の SKIN，セルサイズ依存性をプロットした．粒子数 N が大きいとシステムの体積 Lbox も大きくなる．セルサイズはシステムの一辺を均等に分割するため，それだけシステムが大きい方が分割するセルの間隔を細かくとることができ
る．n + 1 ステップに 1 度行う更新では VNL-CLL 法は CLL 法をベースとするため O(N 2 ) ではなく，O(N ) に計算コストを低減するが，それでもなお VNL 法の方が計算効率は良い．VNL のでは，粒子数が増えると計算コストがリスト更新のため O(N 2 ) のボトルネックを残していたが，CLL 法とハイブリッド化することで更新のボトルネックは CLL の計算コストに置き換わる．
(b) 温度依存性 (図 4.4)
粒子数 N = 1, 000，密度 ρ = 0.8，カットオフ半径 rcut = 2.5 とし，温度 T = 0.2, 0.5, 1.0, 2.0 の4 つの系に対して計算時間の SKIN，セルサイズ依存性をプロットした．システムの温度が高いほ
ど粒子速度は平均的に増加することから，それに伴い粒子がリストの範囲から出る頻度も増加す
る．そのため温度が高い系では，なるべく更新間隔 n を大きくしリストをより長く使い続けた方が計算効率が良くなる．したがって SKIN は大きく取った方が良い．VNL 法のみの場合は低温にするほど，ただ SKIN を長く取ればよいだけであったが，VNL-CLL 法の場合，SKIN はセルサイズ
によって上限が決まっているため，それほど長く取る事はできない．セルサイズを大きく設定し，SKIN を長く取ろうとしても，そもそもセルサイズが大きいと参照体積が大きくなり，計算コストがかかってしまう．これは VNL-CLL 法にすることで更新によるボトルネックが支配的でなくなった結果であると言える．

(c) 密度依存性 (図 4.5)
粒子数 N = 1, 000，温度 T = 1.0，カットオフ半径 rcut = 2.5 とし，密度 ρ = 0.3, 0.5, 0.8, 1.0 の4 つの系に対して計算時間の SKIN，セルサイズ依存性をプロットした．密度によりシステムの体
積 Lbox = V /ρ は変化する．したがって，低密度なほどシステムの一辺の長さは大きくなり，

その分だけセル分割を細かくする事ができる．また，密度はリストに含まれる粒子数に関係する．明らかに，低密度であればリストに記憶する相互作用ペアの数は小さくなり，計算コストも下がる．
また，低密度になるにつれ粒子が移動する平均自由行程が長くなる．これにより，その分だけ粒子はリストの範囲から出やすくなり頻繁に O(N 2 ) の更新を行ってしまう．したがって，更新間隔を長くした方が計算効率は向上する．そのため，低密度なほど最適 SKIN は大きくなるが，温度依存性と同様に SKIN はセルサイズで上限が決まる．

(d) カットオフ半径依存性 (図 4.6)
粒子数 N = 1, 000，温度 T = 1.0，密度 ρ = 0.8 とし，カットオフ半径 rcut = 1.5, 2.0, 2.5, 3.0 の4 つの系に対し計算時間の SKIN，セルサイズ依存性をプロットした．カットオフ半径が大きいほど，カットオフ球体積内に含まれる相互作用を計算するペアの数は多くなる．対応して，SKIN，セルサイズもその分だけ大きく取らなければならない．つまり，リストに含まれる粒子ペアも多くなる．カットオフ半径 rcut が大きくなるにつれ，リスト更新による計算コストの寄与は相対的に小さくなる．

– 31 –

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）

1

1
3.4

0.9
0.8

3.2

rcut + SKIN

rcut + SKIN

3.4

0.7
0.6

3

0.5
0.4

2.8

0.3
0.2

2.6

0.9
0.8

3.2

0.7
0.6

3

0.5
0.4

2.8

0.3
0.2

2.6

0.1
0.25

0.5

Lcell / rcut

1

0.1

2

0.25

N = 500

0.5

Lcell / rcut

1

2

N = 1, 000
1

1
3.4

0.9
0.8

3.2

rcut + SKIN

rcut + SKIN

3.4

0.7
0.6

3

0.5
0.4

2.8

0.3
0.2

2.6

0.9
0.8

3.2

0.7
0.6

3

0.5
0.4

2.8

0.3
0.2

2.6

0.1
0.25

0.5

Lcell / rcut

1

0.1

2

0.25

N = 10, 000
図 4.3

0.5

Lcell / rcut

1

2

N = 100, 000

式 (4.3) の結果：計算時間の SKIN，セルサイズ依存性．T = 1.0，ρ = 0.8，rcut = 2.5 としてそれぞれ粒子数を変化させている．

0.5

0.5
3.4

0.45
0.4

3.2

rcut + SKIN

rcut + SKIN

3.4

0.35
3

0.3
0.25

2.8

0.2
0.15

2.6

0.45
0.4

3.2

0.35
3

0.3
0.25

2.8

0.2
0.15

2.6

0.1
0.25

0.5

Lcell / rcut

1

0.1

2

0.25

T = 0.2

0.5

Lcell / rcut

1

2

T = 0.5
0.5

0.5
3.4

0.45
0.4

3.2

rcut + SKIN

rcut + SKIN

3.4

0.35
3

0.3
0.25

2.8

0.2
0.15

2.6

0.45
0.4

3.2

0.35
3

0.3
0.25

2.8

0.2
0.15

2.6

0.1
0.25

0.5

Lcell / rcut

1

0.1

2

0.25

T = 1.0
図 4.4

0.5

Lcell / rcut

1

2

T = 2.0

式 (4.3) の結果：計算時間の SKIN，セルサイズ依存性．N = 1, 000，ρ = 0.8，rcut = 2.5 としてそれぞれ温度を変化させている．

– 32 –

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）

0.5

0.5
3.4

0.45
0.4

3.2

rcut + SKIN

rcut + SKIN

3.4

0.35
3

0.3
0.25

2.8

0.2
0.15

2.6

0.45
0.4

3.2

0.35
3

0.3
0.25

2.8

0.2
0.15

2.6

0.1
0.25

0.5

Lcell / rcut

1

0.1

2

0.25

ρ = 0.3

0.5

Lcell / rcut

12

ρ = 0.5
0.5

0.5
3.4

0.45
0.4

3.2

rcut + SKIN

rcut + SKIN

3.4

0.35
3

0.3
0.25

2.8

0.2
0.15

2.6

0.45
0.4

3.2

0.35
3

0.3
0.25

2.8

0.2
0.15

2.6

0.1
0.25

0.5

Lcell / rcut

1

0.1

2

0.25

ρ = 0.8
図 4.5

0.5

Lcell / rcut

12

ρ = 1.0

式 (4.3) の結果：計算時間の SKIN，セルサイズ依存性．N = 1, 000，T = 1.0，rcut = 2.5

としてそれぞれ密度を変化させている．

3

0.5

0.5

0.45

0.45

0.4

2.2

rcut + SKIN

rcut + SKIN

2.4

0.35
2

0.3
0.25

1.8

0.2

2.8

0.4
0.35

2.6

0.3
2.4

0.25
0.2

2.2

0.15

1.6

0.15

0.1
0.25

0.5

Lcell / rcut

1

2

2

0.1
0.25

rcut = 1.5

2

4

0.5

0.5

0.45

0.45

0.4

3.2

rcut + SKIN

rcut + SKIN

Lcell / rcut

1

rcut = 2.0

3.4

0.35
3

0.3
0.25

2.8

0.2

3.8

0.4
0.35

3.6

0.3
3.4

0.25
0.2

3.2

0.15

2.6

0.15

0.1
0.25

0.5

Lcell / rcut

1

3

2

0.1
0.25

rcut = 2.5
図 4.6

0.5

0.5

Lcell / rcut

1

2

rcut = 3.0

式 (4.3) の結果：計算時間の SKIN，セルサイズ依存性．N = 1, 000，T = 1.0，ρ = 0.8 としてそれぞれ粒子数を変化させている．

– 33 –

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）

第5章

結論
本論文では，LJ ポテンシャルなどの短距離相互作用をモデルポテンシャルとする短距離古典分子動力学法の計算高速化アルゴリズムに対して最適化技法を構築した．Verlet の近接リスト法 (Verlet Neighbor

List: VNL 法) では，動的リスト更新とそれに伴う計算時間のリスト半径 (SKIN) 依存性，リスト半径の最適化に関して述べた．セル分割法 (Cell Linked List: CLL 法) では，より高度な手法である排他的セル分割法 (Exclusive CLL 法) に対して，計算時間のセルサイズ依存性，セルサイズの最適化を行った．これらに基づき，本論文では新たに VNL 法，CLL 法の両者を組み合わせた手法である VNL-CLL 法の最適化手法を構築した．この最適化手法により，計算システムを決めるパラメータ (粒子数 N , 温度 T , 密度 ρ, カットオフ半径 rcut ) の組み合わせに対し，計算効率を最大化する最適リスト半径 (SKIN)，最適セ
ルサイズの組み合わせが一意に決定される．SKIN= 0.5 に固定した VNL 法、SKIN を最適化した VNL
法、CLL 法、SKIN= 0.5 に固定した VNL-CLL 法、SKIN、セルサイズ共に最適化した VNL-CLL 法に関して計算効率の比較を行った。図 5.1 にそれぞれの手法における計算パフォーマンスについて結果をまとめた．パフォーマンスは単位時間あたりの粒子数 × ステップ数として定義した。よって，この値が高いほど計算性能が良いことになる．
粒子数が小さい系 (N < 1000) に対しては、全粒子探索 O(N 2 ) の効果が少ないため CLL 法に比べ、SKIN を最適化した VNL 法が圧倒的に優れている。粒子数が大きい系に対しては VNL 法は O(N 2 ) のため急激に計算効率が低下してしまう。一方で CLL 法は O(N ) を維持する。

VNL 法、CLL 法の両者の利点を上手く組み合わせることが出来る VNL-CLL 法では、今回、SKIN とセルサイズの最適化を行う事によってさらなる計算効率の向上に成功した。CLL 法と比較して 2 ∼ 3 倍の計算性能を実現する。粒子数の小さな系に対しては VNL 法の恩恵を顕著に受けることができ大きな系に対してもシステムサイズによらず O(N ) を維持することが可能である。
本研究では，シングルコアプロセッサにおける VNL-CLL 法の計算効率最適化を行った．これは並列化によるさらなる大規模計算を見据えた上での最適化パラメータを提供している．並列計算においては一
般に，VNL 法に比べ CLL 法が並列化効率が良い事が知られており，領域 (セル) ごとにスレッドを割り当てることができるためである [18–20]．したがって，本研究における VNL-CLL 法に関しても，同様に高い並列化効率が期待される．並列計算における高速化技法としては，CLL 法で用いる配列のインデックスを並び替える手法 (Pairwise Cell List) [21] や，セルの領域内で粒子のインデックスをソートする手法 (空間ソーティング) [22–26]，さらに GPU を用いた並列計算 [27, 28] など様々な視点からのアプロー
チが盛んに行われている．今後，それらの要素を複合的に取り入れた VNL-CLL 法を実装し，パラメー
タ依存する計算時間の最適化が必須の課題となる．

– 34 –

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）

6

Performance [atom × steps / sec × 10 ]

12

10

8

6

4

2

0

図 5.1

VNL-CLL
SKINopt
VNL-CLL SKIN = 0.5
CLL
VNL
SKINopt
VNL SKIN = 0.5

1000

10000
Number of Particle

100000

各種法におけるパフォーマンス：(赤) 最適化 VNL-CLL 法，(緑)SKIN=0.5 とした VNL-CLL法，(青)CLL 法，(紫) 最適化 VNL 法，(水色)SKIN=0.5 とした VNL 法に対する計算性能を表している．

– 35 –

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）

謝辞
本研究成果は自然科学研究機構計算科学研究センターの計算機を利用して得られたものです．
本論文をまとめるにあたり、多くの方々の御協力を頂きました。皆様に心より感謝申し上げます。
指導教官である金鋼准教授には、本研究を進めるにあたり非常に多くの助言を頂きました。本研究がこの様に一応の決着を迎えられたのは、ひとえに金先生のご尽力のおかげです。この場を借りて深く御礼申し上げます。
東京大学物性研究所の渡辺宙志先生，名古屋工業大学の礒部雅晴先生，東北大学金属材料研究所の芝隼人先生には研究会において，丁寧かつ有益な助言を頂きました．
新潟大学理学部物理学科物性理論研究室の先生方、大野義章教授、吉森明教授、奥西巧一准教授、柳瀬陽一准教授にも研究内容へのご意見や、研究室での活動など様々なご支援、ご協力を頂いております。
物性理論研究室の先輩方にも本研究を進めるにあたり、多くのご助言を頂きました。特に山田武見さん、関孝一さん、中村康晴さんにはプログラミングに関する応用の議論や数値計算の高速化、プレゼン
テーションの流れなど、細かな部分にまで数々の助言を頂きました。ありがとうございます。
同研究室である、樋口沙希さん、池田光佑くんには研究内容の議論から日々のゼミまで、多くの知識を得る事ができました。
同期、後輩のみなさんには、研究意外の場面でも大変お世話になりました。皆様のおかげで博士前期課
程をとても楽しく過ごす事ができました。
最後に、私を支えてくれた家族と友人に感謝します。

– 36 –

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）

参考文献[1] 岡崎 進，「コンピュータ・シミュレーションの基礎」，化学同人 (2000).
[2] 上田 顕，「計算物理入門 分子シミュレーションを中心に」，サイエンス社 (2001)．
[3] 上田 顕，「分子シミュレーション–古典系から量子系手法まで–，裳華房 (2003)．
[4] 田中 實，山本 良一，「計算物理学と計算化学」，海文堂 (1988)．
[5] 佐藤 明，「HOW TO 分子シミュレーション –分子動力学法，モンテカルロ法，ブラウン動力学法，散逸粒子動力学法–」，共立出版 (2004).

[6] J. P. Hansen, I. R. Mcdonald, “Theory of Simple Liquids, Fourth Edition”, Elsevier (2013).
[7] D. J. Evans, Gary. Morriss, “Statistical Mechanics of Nonequilibrium Liquids”, Cambridge University Press (2008).
[8] M. E. Tuckerman, “Statistical Mechanics: Theory and Molecular Simulation”, Oxford University
Press (2010).
[9] D. Frenkel, B. Smit, “Understanding Molecular Simulation”, Academic Press (2002).
[10] M. P. Allen, D. J. Tildesley, “Computer Simulation of Liquids”, Oxford University Press (2009).
[11] D. C. Rapaport, “The Art of Molecular Dynamics Simulation”, Cambridge University Press(1995).
[12] A. A. Chialvo and P. G. Debenedetti, “On the use of the Verlet neighbor list in molecular
dynamics”, Computer Physics Communications 60 (1990) 215-224.
[13] A. A. Chialvo and P. G. Debenedetti, “On the performance of an automated Verlet neighborlist alogritm for large systems on a vector processer”, Computer Physics Communications 64(1991) 15-18.
[14] A. A. Chialvo and P. G. Debenedetti, “An automated Verlet neighbor list algorithm with a multiple time-step approach for the simulation of large systems”, Computer Physics Communications
70 (1992) 467-477.
[15] G. Sutmann, V. Stegailov, “Optimization of neighbor list techniques in liquid matter simulations”, Journal of Molecular Liquids 125 (2006) 197-203.
[16] M. Isobe, “Simple and eﬃcient algorithm for large scale molecular dynamics simulation in harddisk sistem” International Journal of Modern Physics C 10 (1999) 1281-1293.
[17] W. Mattson, B. M. Rice, “Near-neighbor calculations using a modiﬁed cell-linked method”,
Computer Physics Communications 119 (1999) 135-148.
[18] 牛島 省，「OpenMP によるプログラミングと数値計算表」，丸善株式会社 (2006)．
[19] 菅原 清文，「C/C++ プログラマーのための OpenMP 並列プログラミング」，株式会社シナノ(2009)．
– 37 –

《修士論文》

物性研究・電子版 Vol. 5, No. 2, 052601 （2016年5月号）[20] 金田 康正，「並列数値処理 –高速化と性能向上のために–」，コロナ社 (2010)．
[21] P. Gonnet, “Pairwise Verlet Lists: Combining Cell Lists and Verlet Lists to Improve Memory
Locality and Parallelism”, Journal of Computational Chemistry 33 (2012) 76-81.
[22] Z. Yao, J. S. Wang, G. R. Liu and M. Cheng, “Improved neighbor list algorithm in molocularsimulations using cell decomposition and data sorting method”, Computer Physics Communications 161 (2004) 27-35.
[23] S. Meloni, M. Rosati, “Eﬃcient particle labeling in atomistic simulations”, Journal of Computational Chemistry 126 (2007) 121102[24] U. Welliing, G. Germano, “Eﬃciency of linked cell algorithms”, Computer Physics Communications 182 (2011) 611-615.
[25] H. Watanabe, M. Suzuki and N. Ito, “Eﬃcient Implementations of Molecular Dynamics Simulations for Lennard-Jones Systems”, Progress of Theoretical Physics, 126 (2011) 203-235.
[26] H. Watanabe, M. Suzuki and N. Ito, “Huge-scale Molecular Dynamics Simulation of Multibubble
Nuclei”, Computer Physics Communicaitons, 184 (2013) 2775-2784.
[27] J. A. Anderson, C. D. Lorenz, A. Travesset, “General purpose molecular dynamics simulations
fully implemented on graphics processing units”, Journal of Computational Physics, 227 (2008)
5342-5359.
[28] A. J. Proctor, C. A. Stevens and S. S. Cho, “GPU-Optimized Hybrid Neighbor/Cell List Algorithm for Coarse-Grained MD simulations of Protein and RNA Folding and Assembly”, Proceedings of the ACM Conference on Bioinformatics, Computational Biology and Biomedicine, (2013)633-640.

– 38 –