s - 松本充弘の個人的ページ

by user

on 28 марта 2017

Category: Documents

>> Downloads: 1

views

Report

Comments

Description

Download s - 松本充弘の個人的ページ

Transcript

s - 松本充弘の個人的ページ

シミュレーションの全体像・概論：松本充弘
Oct.14, 15:10-16:40
シミュレーションの全体像・概論
松本充弘
京都大学大学院工学研究科機械理工学専攻
現在の専門分野
ミクロスケールの熱流体現象
•
•
•
•
•
分子シミュレーションスクール/TCCIウィンターカレッジ201４
微小気泡，微小液滴，液膜，・・・
界面やイオン系の統計熱力学
固体中のエネルギー輸送
マルチスケールシミュレーション
・・・
岡崎)
1
シミュレーションの全体像・概論：松本充弘
１．本講の目的
２．分子シミュレーションとは
３．モンテカルロ法 vs. 分子動力学法
3.1 統計力学のエッセンス
3.2 モンテカルロ法
3.3 分子動力学法
3.4 その他、さまざまな分子シミュレーション
４．例題：単原子分子系の分子シミュレーション
4.1 分子間の相互作用
4.2 境界条件
4.3 分子動力学シミュレーション
4.4 モンテカルロシミュレーション
4.5 データ解析入門
５．計算効率を上げるには
5.1 物理学的な視点から
5.2 数値計算法の視点から
5.3 計算機の視点から
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
2
シミュレーションの全体像・概論：松本充弘
１．本講の目的
分子シミュレーションの基本事項を説明し，本スクールを受講するにあたって
の共通基盤をつくる，あるいは知識を再確認する．
２．分子シミュレーションとは
狭義には
原子・分子集合体を対象とする数値シミュレーション
対象系のスケール：例えば，1-10 nm, 0.1-10 ns
広義には
粒子集合体を対象とする数値シミュレーション
対象系のスケール：例えば，1 mm-1 m, 1 ms-1 s
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
3
シミュレーションの全体像・概論：松本充弘
本スクールのロードマップ： 2014年度は・・・
注：私が独断で作成した参考図です．
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
4
シミュレーションの全体像・概論：松本充弘
分子シミュレーションの教科書をいくつか・・・
例えば amazon のホームページで
和書 “分子シミュレーション” や洋書 “molecular simulation”
等でキーワード検索してみてください．
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
5
シミュレーションの全体像・概論：松本充弘
１．本講の目的
２．分子シミュレーションとは
３．モンテカルロ法 vs. 分子動力学法
3.1 統計力学のエッセンス
3.2 モンテカルロ法
3.3 分子動力学法
3.4 その他、さまざまな分子シミュレーション
４．例題：単原子分子系の分子シミュレーション
4.1 分子間の相互作用
4.2 境界条件
4.3 分子動力学シミュレーション
4.4 モンテカルロシミュレーション
4.5 データ解析入門
５．計算効率を上げるには
5.1 物理学的な視点から
5.2 数値計算法の視点から
5.3 計算機の視点から
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
6
シミュレーションの全体像・概論：松本充弘
３．モンテカルロ法 vs. 分子動力学法
☞ 詳細は，松林先生の講義を．
3.1 統計力学のエッセンス
☞ 自由エネルギーに関しては，
岡崎先生の講義を．
微視的状態を区別する指標 s
その出現確率
∝ P(s)
物理量 A(s) の統計平均
A
 A(s) P(s)

 P( s)
s
s
主な統計集団
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
我々の課題：適切な P(s) をもつ
ような集団を，いかにして計算
機の中に作り出すか
7
シミュレーションの全体像・概論：松本充弘
３．モンテカルロ法 vs. 分子動力学法
3.2 モンテカルロ法
Monte Carlo法，ＭＣ法
確率過程の考え方を利用して，ボルツマン分布 P(E) ∝ exp(-E/kT) に従う出現確率で
微視的状態を作り出す方法．
２状態モデルを例とすると
遷移行列
ベクトルで略記すると
平衡状態のベクトルが満たすべき条件
平衡状態がボルツマン分布になるための十分条件
すなわち
これが Metropolis MC 法（1953）の原理である．
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
8
シミュレーションの全体像・概論：松本充弘
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
double frand( ) /* Simple Random Number Generator of [0, 1] */
{
return rand( )/(RAND_MAX+1.0);
}
int main( ) //-------------------------------------------------{
int num_step,step,state;
int hist1,hist2;
double eng1=0.0, eng2=1.0;
double temperature=1.0;
while(printf("Number of Steps?")!=0 && scanf("%d",&num_step)==1) {
hist1=hist2=0;
state=1;
/* Initial State */
Number of Steps? 10
10 0.8000 0.2000
Number of Steps? 100
100 0.7300 0.2700
Number of Steps? 1000
1000 0.7340 0.2660
Number of Steps? 10000
10000 0.7289 0.2711
Number of Steps? 100000
100000 0.7305 0.2695
Number of Steps? 1000000
1000000 0.7313 0.2687
Number of Steps? 10000000
10000000 0.7311 0.2689
for (step=0;step<num_step;step++) {
if (state==1) {
if (frand( ) < exp(-(eng2-eng1)/temperature)) {
state=2;
hist2++;
}
else {
hist1++;
}
}
else {
state=1;
hist1++;
}
}
printf("%10d %8.4f %8.4f\n", num_step,hist1/(double)num_step, hist2/(double)num_step);
厳密な比は
1:exp(-1) ～ 0.73106 : 0.26894
}
return 0;
}
２状態モデルのメトロポリスＭＣプログラムの例
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
9
シミュレーションの全体像・概論：松本充弘
３．モンテカルロ法 vs. 分子動力学法
3.2 モンテカルロ法
Metropolisモンテカルロ法
N状態の場合に容易に拡張できて，
これにより，十分なステップが経
過した後は，適切な確率で各事
象が出現することになるから，
いろいろな物理量の平均
 A(s) P(s)
A 
 P( s)
を求めることができる．
s
s
補足
1) 出現頻度が少ない(=エネルギーが高い)状態をいかにうまくサンプ
リングするかが鍵となる．
2) 最適化アルゴリズムとしてよく用いられる Simulated Annealing（焼き
なまし）法，遺伝アルゴリズム法なども，モンテカルロ法の一種と考
えることができる．
3) ボルツマン分布を人工的にずらすことにより，効率や精度を格段に
向上させる種々の拡張アンサンブル法が考案されている．
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
☞奥村先生，志賀先生，
岡本先生の講義を．
10
シミュレーションの全体像・概論：松本充弘
３．モンテカルロ法 vs. 分子動力学法
3.3 分子動力学法
分子動力学法，Molecular Dynamics法，MD法
原子・分子の運動方程式を数値積分することで位相空間中での軌跡を追跡する方法．
やがては熱平衡状態（場合によっては定常状態）に至ると考えて，データ解析を行う．
簡単な例：外場中の１粒子系
⇒
⇒
int main( )
// 1次元調和振動子モデル
{
double x=INIT_X;
double p=INIT_P;
double t=0.0;
double dt, force,omega;
double exact;
// exact solution x=cos(omega*t)
FILE *fout;
fout=fopen("naive.dat","w");
printf("Input dt: ");
scanf("%lf",&dt);
omega=sqrt(SPRING/MASS);
while (t<MAX_TIME) {
t += dt;
force = -SPRING * x;
x += dt*p/MASS;
p += dt*force;
exact=cos(omega*t);
fprintf(fout,"%8.5f %8.5f %8.5f\n",t,x,exact);
}
素朴な前進差分法は使い物にならない！
分子シミュレーションスクール/TCCIウィンターカレッジ201４
fclose(fout);
0;
}
return
岡崎)
11
シミュレーションの全体像・概論：松本充弘
３．モンテカルロ法 vs. 分子動力学法
3.3 分子動力学法
もう少し実用的な数値積分法：前進差分ではなく中心差分をとろう ⇒Leapfrog法
⇒
Leapfrog法による調和振動子計算のサンプルコード
int main( )
{
double x=INIT_X;
double p=INIT_P;
double t=0.0;
double dt, force,omega;
double exact;
// exact solution x=cos(omega*t)
FILE *fout;
fout=fopen("leapfrog.dat","w");
printf("Input dt: ");
scanf("%lf",&dt);
豆知識 Leapfrog
omega=sqrt(SPRING/MASS);
while (t<MAX_TIME) {
直訳すると「蛙跳び」，本来は子供の遊びとしての「馬跳び」．
t += dt;
x += dt*p/MASS;
座標の更新と運動量（速度）の更新を交互に行う様子から，こう
force = -SPRING * x;
p += dt*force;
名付けられたようです．
exact=cos(omega*t);
fprintf(fout,"%8.5f %8.5f %8.5f\n",t,x,exact);
☞ もっと系統的な数値積分法のお
話は，吉井先生の講義を。
}
分子シミュレーションスクール/TCCIウィンターカレッジ201４
} 岡崎)
fclose(fout);
return 0;
12
シミュレーションの全体像・概論：松本充弘
３．モンテカルロ法 vs. 分子動力学法
3.4 その他、さまざまな分子シミュレーション
原子・分子を「粒子」と考えて，その統計的性質を数値計算で調べるＭＣ法・ＭＤ法は，
1950年代以降，コンピュータの発達とともに急速に広まった．現在のフロンティアは・・・
より精密に
元のＭＣ法・ＭＤ法はポテンシャル関数が与えられていることが前提である．粒子配置が更新さ
れるたびに量子力学方程式を解いてエネルギーや力を精密に求めなおすというアイディアは，第
一原理ＭＤ法として結実している．
より複雑な系を
タンパク質や界面活性剤系に代表される，膨大な内部自由度を有する分子の集合体をなるべく
余分な近似なしに直接計算しようという機運が高まっている．
☞ 高橋先生，長岡先生の講義
より大規模な系を
物理現象の階層性を利用した粗視化シミュレーション，ブラウニアンダイナミクスシミュレーション，
散逸粒子動力学計算，マルチスケールシミュレーション，ハイブリッド計算などの手法の開拓が
さかんである．
☞ 篠田先生の講義
連続体力学からのアプローチ
メッシュフリー数値解析の１つとして分子シミュレーションの手法を
発展させた粒子シミュレーション法（MPS法，SPH法など）の応用が
広がっている．
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
13
シミュレーションの全体像・概論：松本充弘
１．本講の目的
２．分子シミュレーションとは
３．モンテカルロ法 vs. 分子動力学法
3.1 統計力学のエッセンス
3.2 モンテカルロ法
3.3 分子動力学法
3.4 その他、さまざまな分子シミュレーション
４．例題：単原子分子系の分子シミュレーション
4.1 分子間の相互作用
4.2 境界条件
4.3 分子動力学シミュレーション
4.4 モンテカルロシミュレーション
4.5 データ解析入門
５．計算効率を上げるには
5.1 物理学的な視点から
5.2 数値計算法の視点から
5.3 計算機の視点から
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
14
シミュレーションの全体像・概論：松本充弘
４．例題：単原子系の分子シミュレーション
一般論はこれぐらいにして，簡単な系を例にとって，実際のプログラミングやデータ
解析の感じをつかむことにしよう．他の講義を聴く際に，イメージがわくように・・・
自分たちでソフトウェアを作る方々へ：
ところで，プログラミング言語は何がいいか？
・コンピュータが正しく動くなら原理的には何でもＯＫ
浮動小数点演算ができる，構造化プログラミング，オブジェクト指向，…
・趣味的なプログラミングなら
Basic, Visual Basic, …
・過去の蓄積を重視するなら
Fortran77, Fortran90, …
・将来性や計算性能を考えると
C, C++, …
最近こんな本が出版
・その他 Javaなどでウェブアプリ化を目指す人も・・・
されました．おもしろ
そうですね．
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
15
シミュレーションの全体像・概論：松本充弘
４．例題：単原子系の分子シミュレーション
4.1 分子間の相互作用
粒子間に適切な相互作用ポテンシャルを仮定することが，最も重要である．
⇒ここの検討が不十分だと，想定している「もの」や「計算条件」に現実性がなくなる！
原理的な考え方
N原子系のポテンシャルエネルギーを，“Taylor展開” して考える．
☞ 三上先生の講義
高次の項：多体相互作用
２次の項：二体相互作用，pair-wise potential
１次の項：電場や重力場などの外場をあらわす．普通はゼロ
定数項：エネルギーの原点の選び方の話だから普通はゼロとして
よい
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
16
シミュレーションの全体像・概論：松本充弘
４．例題：単原子系の分子シミュレーション
4.1 分子間の相互作用
ここで扱う例題では，二体力相互作用のみとし，さらに，原子間距離だけの関数と
仮定する．その代表例が，Lennard-Jones相互作用である．
(i,j)と(j,i)は同じペアを表すので，２重に数えないように×1/2
電子雲の重なりに
由来する反発
誘起双極子による引力
（分散力）
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
17
シミュレーションの全体像・概論：松本充弘
４．例題：単原子系の分子シミュレーション
4.1 分子間の相互作用
ポテンシャルエネルギーが与えられれば，そのgradientから，各粒子にはたらく力
が計算できる！
粒子間の力を計算するコード例
for (i=0; i<NUM_ATOM; i++) {
for (j=i+1; j<NUM_ATOM; j++) { //作用反作用の法則を考慮
dx=posx[i]-posx[j];
dy=posy[i]-posy[j];
dz=posz[i]-posz[j];
r2=dx*dx+dy*dy+dz*dz;
ri06=1.0/(r2*r2*r2);
ri12=ri06*ri06;
eng=4*(ri12-ri06);
frc=4*(12*ri12-6*ri06)/r2;
eng_pot
frcx[i]
frcy[i]
frcz[i]
frcx[j]
frcy[j]
frcz[j]
} }
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
+=
+=
+=
+=
-=
-=
-=
eng;
frc*dx;
frc*dy;
frc*dz;
frc*dx;
frc*dy;
frc*dz;
作用反作用の法則を使うと
計算量が半分に
18
シミュレーションの全体像・概論：松本充弘
４．例題：単原子系の分子シミュレーション
4.2 境界条件
現実の系に比べてごく少数の原子・分子を使って，対象系の物性を調べるためには，
粒子を箱に閉じ込めることはできない．壁の影響を強く受けてしまうから．
周期境界条件 Periodic Boundary Conditions
・有限個の粒子が有限体積の計算セルに閉じ込められて
いる．
・計算セルは隙間なく３次元的に並べられていると仮定する．
・粒子は計算セルの端から何ら影響を受けずに運動するが，
計算セルから出た粒子は反対側の境界から同じ運動量
で戻ってくると考える．あるいは，隣のセルのゴースト粒
子が実体化すると考えても同じ．
補足：「壁の影響」を見たい場合も，もちろんある．
その時は，例えば２次元的な周期境界条件を考
えればよい．
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
19
シミュレーションの全体像・概論：松本充弘
４．例題：単原子系の分子シミュレーション
4.2 境界条件
周期境界条件を使う場合の，相互作用の計算に注意
粒子は計算セルに影響されずに自由に動く．ということは，隣の計算セル内の粒子と
相互作用することになる．したがって，実
粒子間の力を計算するコード例
質的に同じペアの相互作用を繰り返し足
for (i=0; i<NUM_ATOM; i++) {
for (j=i+1; j<NUM_ATOM; j++) { //作用反作用の法則
し合わせることがないよう，相互作用の計
dx=posx[i]-posx[j];
if (dx<-cxh) dx+=CELL_X;
算はセルサイズの半分以下で打ち切らな
else if (dx>cxh) dx-=CELL_X;
dy=posy[i]-posy[j];
ければならない (minimal image convention)．
if (dy<-cyh) dy+=CELL_Y;
else if (dy>cyh) dy-=CELL_Y;
dz=posz[i]-posz[j];
if (dz<-czh) dz+=CELL_Z;
else if (dz>czh) dz-=CELL_Z;
r2=dx*dx+dy*dy+dz*dz;
ri06=1.0/(r2*r2*r2);
ri12=ri06*ri06;
eng=4*(ri12-ri06);
frc=4*(12*ri12-6*ri06)/r2;
eng_pot
frcx[i]
frcy[i]
frcz[i]
frcx[j]
frcy[j]
frcz[j]
}
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
}
+=
+=
+=
+=
-=
-=
-=
eng;
frc*dx;
frc*dy;
frc*dz;
frc*dx;
frc*dy;
frc*dz;
20
シミュレーションの全体像・概論：松本充弘
４．例題：単原子系の分子シミュレーション
4.3 分子動力学シミュレーションのコード例
典型的なＭＤシミュレーションプログラムは次のような構造をしている：
int main( )
{
int step;
printf("Input Density :");
scanf("%lf",&density);
printf("Input Target Temperature :");
scanf("%lf",&t_target);
CELL_X=CELL_Y=CELL_Z=pow(NUM_ATOM/density, 1/3.0);
printf("NUM_ATOM=%d\n",NUM_ATOM);
printf("CELL_X=CELL_Y=CELL_Z=%f\n",CELL_X);
initial( );
initplot( );
for (step=0; step<=TOTAL_STEP; step++) {
force( );
move( );
if (step%SAVE_STEP==0) {
statistics(step);
}
}
finalplot( );
return 0;
}
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
21
シミュレーションの全体像・概論：松本充弘
ちょっとした実習
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
22
シミュレーションの全体像・概論：松本充弘
４．例題：単原子系の分子シミュレーション
4.4 モンテカルロシミュレーションのコード例
同じく，典型的なＭＣシミュレーションプログラムは次のような構造をしている：
int main( )
{
int step;
printf("Input Target Temperature :");
scanf("%lf",&t_target);
initial( );
initplot( );
for (step=0; step<=TOTAL_STEP; step++) {
update( );
if (step%SAVE_STEP==0) statistics(step);
}
finalplot( );
return 0;
}
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
23
シミュレーションの全体像・概論：松本充弘
４．例題：単原子系の分子シミュレーション
4.5 データ解析入門
手間のかかるプログラム開発と，コストのかかるシミュレーション実行が終わると，
いよいよ，結果の解析にとりかかる．元の研究目的に応じて，データ解析の手法
はさまざまであるが，至福の時のはず・・・
と，その前に．
あなたが目指すシミュレーションは次のどれにあてはまるだろうか？
平衡状態 equilibrium：集団平均，時間平均ともに可能
定常状態 steady state：時間平均が可能
非平衡状態 non-equilibrium
非定常状態 non-steady state：統計平均をとるのは厄介
9.5 nm
定常計算の例
平衡計算の例
温度勾配中の気泡移動
微小気泡の物性
分子シミュレーションスクール/TCCIウィンターカレッジ201４岡崎)
非定常計算の例
沸騰現象の初期過程
24
シミュレーションの全体像・概論：松本充弘
４．例題：単原子系の分子シミュレーション
4.5 データ解析入門
温度の話
☞ 奥村先生の講義
・MCシミュレーションなら普通，温度は我々が与えるべき外部パラメタ
・MDシミュレーションなら・・・
熱平衡が仮定できるなら，温度は運動エネルギーの平均値
しかし，これは我々が望んでいた温度とは異なるだろう．
方法 1) 速度スケーリングで、希望する温度に近づけていく
方法 2) 「熱浴」によって温度を制御する：能勢-Hoover熱浴，Langevin熱浴，…
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
25
シミュレーションの全体像・概論：松本充弘
４．例題：単原子系の分子シミュレーション
4.5 データ解析入門
圧力の話
☞ 奥村先生の講義
ビリアルによる圧力の表式
圧力一定の統計集団を作るためには，体積の揺らぎが必要である．
例：LJ流体の圧力の数密度依存性
補足
1) 系に異方性があれば，ビリアルは一般に
テンソル量となる．その各成分を局所的
に，かつ精度良く求めれば，界面張力な
どの情報が得られる．
2) 粒子間が引力的であれば，ビリアル項は
負になり，場合によっては大きな負圧
negative pressure を生じる．
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
26
シミュレーションの全体像・概論：松本充弘
４．例題：単原子系の分子シミュレーション
4.5 データ解析入門
● 分子シミュレーションを行えば，原子・分子の座標や速度など詳細なデータを求
めることができるが，データ量が膨大になるため，そのすべてをハードディスクな
どに保存することは難しい．
● どんなデータ解析を行うか，計画を立ててからシミュレーションを行わないと，必
要なデータが足りなくて，シミュレーションをやり直すような事態に陥る！
（ポイント）データの種類，データストアの頻度，ディスク容量
構造解析の例：
動径分布関数，構造因子，分子の配向性，・・・
ダイナミクス解析の例：
自己速度相関関数，平均２乗変位，・・・
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
27
シミュレーションの全体像・概論：松本充弘
１．本講の目的
２．分子シミュレーションとは
３．モンテカルロ法 vs. 分子動力学法
3.1 統計力学のエッセンス
3.2 モンテカルロ法
3.3 分子動力学法
3.4 その他，さまざまな分子シミュレーション
４．例題：単原子分子系の分子シミュレーション
4.1 分子間の相互作用
4.2 境界条件
4.3 分子動力学シミュレーション
4.4 モンテカルロシミュレーション
4.5 データ解析入門
５．計算効率を上げるには
5.1 物理学的な視点から
5.2 数値計算法の視点から
5.3 計算機の視点から
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
28
シミュレーションの全体像・概論：松本充弘
５．計算効率を上げるには
5.0 計算速度を支配するもの
計算屋の常識
1) 無駄な計算はしない
例：作用反作用の法則を使って計算量を半分に
＊しかし，並列計算の場合，あえて (i, j) と (j, i) をそれぞれ計算することもある．
2) 計算量の多い部分をつきとめ，集中的に対処する
プロファイラなどのツールを使って，関数やサブルーチンごとの時間を計測し，重要な
部分から効率向上を図ろう．
分子シミュレーションでは，通常，相互作用の計算が最も
時間がかかる．なぜなら，粒子数の２乗に比例するから．
（２体相互作用の場合）
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
29
シミュレーションの全体像・概論：松本充弘
５．計算効率を上げるには
5.1 物理学的な視点から
（１）相互作用のカットオフ (cutoff, truncation)
通常，相互作用は粒子間距離が大きくなるとともに減衰する．そこで，予め，あ
る距離 rc を決めておき，それより遠くの相互作用は計算しないことにする．
留意点
1) カットオフ距離の選び方
2) カットオフ時の補正
3) クーロン力など減衰が遅い相互作用
はカットできない！
が無視できないから
☞ Ewald法（三上先生の講義）
4) でも，カットオフだけではあまり時間
の短縮にはならない
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
30
シミュレーションの全体像・概論：松本充弘
５．計算効率を上げるには
5.1 物理学的な視点から
（２）大規模系での相互作用計算のサボり方
・帳簿法 Book keeping method
それぞれの粒子について，相互作用する可能性の
ある相手粒子をあらかじめリストアップしておく．リス
トは粒子の移動状況に応じて「時々」更新する．
・セル分割法 Cell division method
計算系を rc 程度のサブセルに分割しておく．相互
作用は，１つのサブセル内，および隣り合うサブセル
の粒子間のみを考えればよい．
さらに大規模な系や長距離相互作用がある系に対
しては
・粒子-セル法 Particle-particle & particle-cell method
・多重極子展開法 Multi-pole expansion method
・・・
☞ Ewald法（三上先生の講義）
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
31
シミュレーションの全体像・概論：松本充弘
５．計算効率を上げるには
5.1 物理学的な視点から
こうした工夫により，計算時間が粒子数にほぼ比例するようになれば大成功！
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
32
シミュレーションの全体像・概論：松本充弘
５．計算効率を上げるには
5.2 数値計算法の視点から
ＭＤ法において，なるべく時間刻み Dt を大きくしよう
・シンプレクティック積分法
エネルギー保存則を満足しつつ，系統的に数値積分の精度を上げることができる．
☞ 吉井先生の講義
・マルチタイムステップ法
特徴的な時間が大きく異なる複数の運動がある場合，普通なら，最も速い運動に合
わせて Dt を与えなければならない．速い運動には短いDt を，遅い運動には長いDt
を使うと，それほど精度を落とさずに計算できるはず…
⇒とくにマルチスケール, 粗視化シミュレーションに有効
☞ 篠田先生の講義
（例）圧力波による気泡崩壊のシミュレーション
気泡から離れたバルク液体領域は，20Dtで流体
計算をおこなっている
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
33
シミュレーションの全体像・概論：松本充弘
５．計算効率を上げるには
5.3 計算機の視点から
（１）高速なハードウェアを使おう
2014年6月の TOP500 リストから
1) 京コンピュータは現在，第4位．
2) 東京工業大のTSUBAMEが13位, 核融研のheliosが30位, …
3) 中国が1位，14位，…
上位はすべて数万コア以上の並列分散型，GPUやCELLが大活躍
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
34
シミュレーションの全体像・概論：松本充弘
理化学研究所
計算科学研究機構のウェブページから
☞安藤先生がご専門・・・
【計算ノード群の構成】
ノード数：88,128
ピーク演算性能：10.62PFLOPS
メモリ総容量： 1.26PB（ノード当り16GB）
計算ノード間ネットワーク：6次元メッシュ/トーラス
帯域： 3次元の正負各方向にそれぞれ5GB/s x 2（双方向）
【プロセッサ構成】
8コア構成、各コア256本の浮動小数点レジスタを備えたスーパー
スカラ方式
・SIMD拡張（積和演算器2個 x 2セット）
コア共有の2次キャッシュ（6MB，12way）
・ハードウェアバリア機構
・プリフェッチ機構
・セクタキャッシュ機能
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
35
シミュレーションの全体像・概論：松本充弘
５．計算効率を上げるには
5.3 計算機の視点から
（２）並列計算機をうまく使おう
並列計算機の基礎知
☞安藤先生の講義
○分散メモリ型
ハードウェアのコストパフォーマンスが良い
並列計算システムの主流
（ＰＣクラスターからスパコンまで）
データ管理は原則としてユーザーが面倒を見る
GPUやCELLもこの仲間
記述言語：PVM，MPI，CUDA，OpenCL，・・・
○共有メモリ型
メモリを共有する特別な仕掛けが必要
大規模な並列計算には不向き
マルチコアプロセッサでは使える技術
（Quad Core など）
記述言語：OpenMP が便利
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
36
シミュレーションの全体像・概論：松本充弘
OpenMP を使った簡単な実習
#include <omp.h>
#include <stdio.h>
int main()
{
#pragma omp parallel
printf("Hello from thread %d, nthreads %d\n", omp_get_thread_num(), omp_get_num_threads());
return 0;
}
簡単な実習
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
37
シミュレーションの全体像・概論：松本充弘
OpenMP を使ったMDプログラムの並列化：素朴な例
for (i=0; i<NUM_ATOM; i++) {
fix=fiy=fiz=0.0;
#pragma omp parallel for private (dx,dy,dz,r2,ri06,ri12,eng,frc) reduction(+:eng_pot,fix,fiy,fiz)
for (j=i+1; j<NUM_ATOM; j++) {
dx=posx[i]-posx[j];
if (dx<-cxh) dx+=CELL_X; else if (dx>cxh) dx-=CELL_X;
並列化前のforce計算部分
dy=posy[i]-posy[j];
for (i=0; i<NUM_ATOM; i++) {
if (dy<-cyh) dy+=CELL_Y; else if (dy>cyh) dy-=CELL_Y;
for (j=i+1; j<NUM_ATOM; j++) {
dz=posz[i]-posz[j];
dx=posx[i]-posx[j];
if (dz<-czh) dz+=CELL_Z; else if (dz>czh) dz-=CELL_Z;
if (dx<-cxh) dx+=CELL_X; else if (dx>cxh) dx-=CELL_X;
r2=dx*dx+dy*dy+dz*dz;
dy=posy[i]-posy[j];
if (dy<-cyh) dy+=CELL_Y; else if (dy>cyh) dy-=CELL_Y;
dz=posz[i]-posz[j];
if (dz<-czh) dz+=CELL_Z; else if (dz>czh) dz-=CELL_Z;
r2=dx*dx+dy*dy+dz*dz;
ri06=1.0/(r2*r2*r2);
ri12=ri06*ri06;
if (r2<RCUT*RCUT) {
ri06=1.0/(r2*r2*r2);
ri12=ri06*ri06;
eng=4*(ri12-ri06);
frc=4*(12*ri12-6*ri06)/r2;
}
}
}
}
eng_pot += eng;
fix += frc*dx;
fiy += frc*dy;
fiz += frc*dz;
frcx[j] -= frc*dx;
frcy[j] -= frc*dy;
frcz[j] -= frc*dz;
eng=4*(ri12-ri06);
frc=4*(12*ri12-6*ri06)/r2;
eng_pot += eng;
frcx[i] += frc*dx;
frcy[i] += frc*dy;
frcz[i] += frc*dz;
frcx[j] -= frc*dx;
frcy[j] -= frc*dy;
frcz[j] -= frc*dz;
frcx[i]+=fix;
frcy[i]+=fiy;
frcz[i]+=fiz;
分子シミュレーションスクール/TCCIウィンターカレッジ201４
}}
岡崎)
38
シミュレーションの全体像・概論：松本充弘
最近の話題 Xeon Phi コプロセッサー
http://ascii.jp/elem/000/000/818/818030/ より
既存のデータ解析とビッグデータの大きな差である処理速度や解像度。このギャップを埋めるに
は、既存のサーバーアーキテクチャでは難しい。ハイパースケールな拡張性を実現すべくインテ
ルが開発したのが、「Xeon Phiコプロセッサー製品ファミリー」である。・・・ Xeon Phiは、超並列
処理を実現するHPC（High Performance Computing）向けのコプロセッサーユニットだ。最大61の
プロセッサーコアとローカルメモリをPCIボードに搭載し、メインCPUからの処理をオフロードするコ
プロセッサーとして動作する。
NECの試算では、一般的な 2ソケット
サーバー120台を用いていた処理を、
1/10の台数で行なえるという。また、
設置面積は約1/3、消費電力は約1/6
に削減できる。さらにXeon Phi の採
用によって、独自の開発環境やスキ
ルが不要となるため、アプリケーション
の移行コストも大幅に削減することが
可能だ。
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
39
シミュレーションの全体像・概論：松本充弘
ここまでで，分子シミュレーションに関する基
礎事項の説明／復習はすべて終わりました．
このスクールで行われる他の講義への興味をか
きたてることができたなら幸いです…
１．本講の目的
２．分子シミュレーションとは
３．モンテカルロ法 vs. 分子動力学法
3.1 統計力学のエッセンス
3.2 モンテカルロ法
3.3 分子動力学法
3.4 その他，さまざまな分子シミュレーション
４．例題：単原子分子系の分子シミュレーション
4.1 分子間の相互作用
4.2 境界条件
4.3 分子動力学シミュレーション
4.4 モンテカルロシミュレーション
4.5 データ解析入門
５．計算効率を上げるには
5.1 物理学的な視点から
5.2 数値計算法の視点から
5.3 計算機の視点から
４章，５章の主なサンプルプログラムは
http://www.mitsuhiromatsumoto.mech.kyoto-u.ac.jp/
からダウンロードできます．
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
40
シミュレーションの全体像・概論：松本充弘
分子シミュレーションの教科書が数多く出版されています．ぜひもっと深く学んでいってください．
本スクールのロードマップ
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
41
シミュレーションの全体像・概論：松本充弘
１．本講の目的
２．分子シミュレーションとは
３．モンテカルロ法 vs. 分子動力学法
3.1 統計力学のエッセンス
3.2 モンテカルロ法
3.3 分子動力学法
3.4 その他，さまざまな分子シミュレーション
４．例題：単原子分子系の分子シミュレーション
4.1 分子間の相互作用
4.2 境界条件
4.3 分子動力学シミュレーション
4.4 モンテカルロシミュレーション
4.5 データ解析入門
５．計算効率を上げるには
5.1 物理学的な視点から
5.2 数値計算法の視点から
5.3 計算機の視点から
４章，５章の主なサンプルプログラムは
http://www.mitsuhiromatsumoto.mech.kyoto-u.ac.jp/
からダウンロードできます．
分子シミュレーションスクール/TCCIウィンターカレッジ201４
岡崎)
42