j x j j j + 1 l j l j = x j+1 x j, n x n x 1 = n 1 l j j=1 H j j + 1 l j l j E

高分子計算科学B

自由連結鎖

1 2 3 4 5 6 7 8 9 高分子のモデルとして自由連結鎖を考える．これは，ボンドの結合長が一定値であるが， 隣り合う二つ のボンドベクトルの成す角が全く任意であるような鎖である． j 番目の原子団の位置を xjとし，j 番目の原子団 j + 1 番目の原子団の位置を結ぶベクトルを, ljとする と， lj = xj+1− xj, (1) 今，高分子が n 個の原子団からなるとき，末端の原子団同士を結ぶベクトルは次のようになる． xn− x1 = n_∑−1 j=1 lj (2) 考えている系が1本の分子鎖のみから成るとき，系のエネルギー H は, この分子鎖の今フォーメー ション のみによって決まる． 原子同士の重なりを無視し， j 番目の原子団と j + 1 番目の原子団を繋ぐベクトル lj が， 電気双極子を 持っており，その向きが ljと平行なとき，外部から電場 E が加え られた場合には，この系のエネルギー

は， H =− n_∑−1 j=1 aE· lj, (3) となる，ただし，a は定数である．

1次元の自由連結鎖

簡単のために，自由連結鎖を1次元空間で考える． 式(1)は， lj = xj+1− xj, (4) 式(2)は， xn− x1 = n_∑−1 j=1 lj, (5) と書き直すことができる．ただし， lj =±l (l > 0) である．また，式(3)は， H =−ϵ n_∑−1 j=1 lj, (6) と書き直すことができる．ここで， ϵ は，電場と双極子の大きさによって決まる定数である．

模擬試験1

ex _{を x = 0 のまわりでテーラー展開すると，} ex = ∞ ∑ j=0 xj j! (7) = 1 + x +1 2x 2₊1 6x 3_{+· · · +} 1 j!x j_{+· · · , (8)} である．ただし，0! = 1 とする．x = 1 とすると， e = 1 + 1 +1 2 + 1 6 +· · · + 1 j! +· · · (9) = ∞ ∑ 1 , (10)

となる．これを，次の式で近似する. e = n ∑ j=0 1 j!. (11) n を, 1 から 15 まで変化させ，それぞれの n について， 式(11)の値を出力するプログラムを作れ．

メトロポリスの方法

統計力学の結論の一つに次のものがある． 体積一定の系が，温度 T の環境にある場合を考える． 系の 微視状態に番号を振り， j 番目の微視状態のエネルギーが Hj で表されるとき， j 番目の微視状態が現れ る確立 Ψj は， Ψj = 1 Z exp ( − Hj kBT ) , (12) で表される．ただし，kBはボルツマン定数である． また, Z =∑ j exp ( − Hj kBT ) , (13) である．式(13)において，∑ j は全ての微視状 態についての和を表す． 数値計算によって，式(12)の確率分布に従うように微視状態を発生 することができれば，その微視状態 についてエネルギー，圧力など各種物理量 の平均を取ることで，温度 T における熱力学的な物理量の平 均を得ることが できる． ある母集団について，その母集団に関連した平均や和を計算するために，計算 機の中で乱数を発生さ せ，その母集団を近似的に生成するさいに使われる方法に モンテカルロ法がある．モンテカルロ法の中で も，式(12)に従う統 計力学的な母集団を生成する方法をメトロポリス法と呼ぶ．

詳細釣り合いの原理

式(12)の確立で現れる微視状態を生成する際に，基本となるのは詳 細釣り合いの原理(principle of detailed balance)である．微視状態 a から 微視状態 b へ遷移する確率をPab とすると，平衡状態にて， 式 (12)の確立に従い微視状態が現れる場合， exp ( −Ha kBT ) Pab  = exp ( − Hb kBT ) Pba, (14) が成り立つ．これが詳細釣り合いの原理と呼ばれる関係である．

N 個の粒子系の場合

N 個の粒子から成る系の場合，ハミルトニアン(エネルギー) H は， 3N 個の粒子の座標 rN = (r1, r2,· · · rN) と，それに共 役な 3N 個の運動量 pN = (p1, p2,· · · pN) の関数とし て， H = N ∑ j=1 p2_j 2m+ U (r1, r2,· · · , rN), (15) と 表すことができる．ただし，U はポテンシャルエネルギーである． なお， rj = (xj, yj, zj), pj = (pxj, pyj, pzj) であることに注意する．ただし， xj, yj, zj はそれぞれ，j 番目の粒子の位置の x座標，y座 標, z座標．また，それ らに共役な運動量を pxj, pyj, pzj にて表している． この場合，pN_{, r}N _{のを決めると一つの微視状態が決まる． 微視状態が， r}N _{= (r1, r2,· · · , r} N), pN = (p₁, p₂,· · · , p_N) と r + drN = (r1 + dr1, r2+ dr2,· · · , rN + drN), p + dpN = (p1+ dp1, p2 + dp₂,· · · , p_N + dp_N) の間に見出される確率は， exp ( − H kBT ) drN_dpN ₍₁₆₎ にて与えられる．密度分布 exp ( − H kBT ) をカノニカル分布という． 系の体積と温度が一定で，ある力学量 Q が，粒子の座標 rN _{だけの関数であれば， Q の熱力学平均} < Q > は， < Q >= ∫ ∫ Q exp ( − H kBT ) drN_dpN ∫ ∫ exp ( − H kBT ) drN_dpN = ∫ Q exp ( − U kBT ) drN ∫ exp ( − U kBT ) drN , (17) にて与えられる．

N 粒子系の場合のアルゴリズム

適当な初期条件を始めとして，次のような過程の繰り返しで次々と rN _{を更新していく．} (1) k回目の更新で，系が状態 t をとるとする． これを，k番目の状態が t であると言うことにする． (2) この状態 t から転移する状態 m を次のように生成する． N 個の粒 子から 1 つをランダムに選び出す．この分子を aとし，位置座標の成分, xa, yb, zb を次のよう

に少しだけ変化させる． xa → xa+ δx(1− 2ξ1), ya → ya+ δy(1− 2ξ2), za → za+ δz(1− 2ξ3). ここで，δx, δy, δz は，x−, y−, z− 方向について 決めた微小な定数である．また， ξj は，[0, 1] の範囲の一 様乱数である． すなわち，図1で示されるように，でたらめに分子を選んで，その分 子を現在の位置を中 心とした，8δx× δy× δzの直 方体のなかにでたらめに移動させるということを行い，状態 t から状態 m を 作りだす．(3) 微視状態 t から m へ，詳細釣り合いの関係式(14) を満たす形で遷移させる．その一つの 方法は，まず，状態 t のエネルギー Utと，状態 m のエネルギー Um を比べる． Ut≥ Um ならば，k + 1 番目の状態として m を取る． Ut< Umならば，サイコロを振って，m に移るかどうかを決めるようなこ とをやる．すなわち， [0, 1] の 一様乱数 η を発生させ， η < exp ( −Um− Ut kBT ) (18) ならば，m を k + 1番目の状態として採用し， η > exp ( −Um− Ut kBT ) (19) ならば，t のまま，これを k + 1 番目の状態とする． (4) (1)へ戻る． この(1)から(4)までの操作を粒子の数 N 回繰り返すことをサンプリングの単 位として，1モンテカルロ (MC)時間と呼ぶことがある．

1次元の自由連結鎖のモンテカルロシミュレーション

熱統計力学の分野でモンテカルロシミュレーションといえば，メトロポリス法の ことを指すと考えてよ い．ここでは，先に説明した1次元の自由連結鎖について， メトロポリス法で微視状態を生成する方法を 記述する． 式(4)‒(6)において，lj =±1と仮定する． 適当な初期条件を設定しそれを 0 番目の微視状態とする．例

t

m

Figure 1: N 粒子系のメトロポリス法での座標の更新 えば，全 ての lj を 1 とする． (1) k 番目の更新で，系が微視状態 t をとるとする． (2) この状態 t から転移する状態 m を次のように生成する. 1 から n− 1 までの数をでたらめに選ぶ． 選んだ数を a とする．laは，+1 か −1 である，la が +1 な ら_{−1 に，−1 なら+1 になるように，la}に _{−1 をかける．こうして状 態 t から作った状態を m とする．す} なわち，ボンドをでたらめに選んで， そのボンドの向きを逆にした状態を m とする． (3) まず，状態 t のエネルギー Htと，状態 m のエネルギー Hmを比べる． Ht≥ Hm ならば，k + 1 番目の状態として m を取る． Ht< Hm ならば，[0, 1] の一様乱数 η を発生させ， η < exp ( −Hm− Ht kBT ) (20) ならば，m を k + 1番目の状態として採用し， η > exp ( −Hm− Ht kBT ) (21) ならば，t のまま，それを k + 1 番目の状態とする． (4) (1)へ戻る．

prog6.c

prog6.c として次のプログラムを作り，コンパイルおよび実行せよ． 1 #include <stdio.h> 2 #include <stdlib.h> 3 4 #define N (20) 5 6 double uran(){ 7 return((double)rand()/RAND_MAX); 8 } 9 10 int main() 11 { 12 float rn[N]; 13 int l[N]; 14 int j; 15 16 for(j=0;j<N;j++){ 17 rn[j] = uran(); 18 if(rn[j] < 0.5){ 19 l[j] = -1; 20 }else{ 21 l[j] = 1; 22 } 23 24 printf("%f,%2d\n",rn[j],l[j]); 25 } 26 return(0); 27 } 28

mc.c

次のプログラムは1次元自由連結鎖に，エネルギーとして 式(6) の H を用いた場合のモンテカルロシミュ レーションの一例である． ただし， ϵ = 6.0× 10−22_{[J], k} B= 1.380662×−23[J/K] とした． 1 #include <stdio.h> 2 #include <stdlib.h> 3 #include <math.h> 4 5 #define N (100) 6 #define EP (6.0e-22) /* [J] */

7 #define KB (1.380662e-23) /* The Botzmann's constant [J/K] */ 8 #define MCTIME (1000) 9 10 double uran(){ 11 return((double)rand()/RAND_MAX); 12 } 13 14 int main() 15 { 16 int l[N],p[N+1]; 17 int k,j,mct,lbefore,end; 18 float dene,t; 19 FILE *fp; 20 21 printf("Temperature ="); 22 scanf("%f",&t); 23 printf("\n"); 24 25 for(j=0;j<N;j++){ 26 l[j] = 1; 27 } 28 29 for(j=0;j<N+1; j++){ 30 p[j] = 0;

31 } 32 33 for(mct=0;mct<MCTIME;mct++){ 34 for(k=0;k<N;k++){ 35 j = (N+1)*uran(); 36 lbefore = l[j]; 37 l[j] *= (-1);

38 dene = - EP*l[j] - (-EP*lbefore); 39 if(dene > 0){ 40 if(exp(-dene/(KB*t)) < uran()){ 41 l[j] = lbefore; 42 } 43 } 44 } 45 46 for(j=0,end=0;j<N;j++){ 47 end += l[j]; 48 } 49 p[(end+N)/2]++; 50 51 printf("%d\n",mct); 52 } 53 54 fp = fopen("result.txt","w"); 55 56 for(j=0;j<(N+1);j++){ 57 fprintf(fp,"%d\t%e\n",2*j-N,((float)p[j])/MCTIME); 58 } 59 60 fclose(fp); 61 62 return(0); 63 } 64

65 このプログラムでは，結果が result.txt というファイルに出力されるようになっている． 温度はプログラ ムの実行中に入力するようになっている． 温度を，10, 100, 1000 として行った結果を $ mv result.txt 10.txt などによって，名前を result.txt から，10.txt, 100.txt, 1000.txt に変更して保存し，各温度の結果 をグラフ上にプロットして比較せよ．

mc.cについて

1‒3 行目: このプログラムで使用する関数，定数の情報を含むファイルを指定する． 10 ‒ 12 行目: [0, 1]の一様乱数を発生させる関数 uran() の定義．10行目を見ると， この関数が double 型の値を返し すことが分かる. 14 ‒ 63行目: メイン関数．C言語のプログラムは最初に，メイン関数から実行される． 16行目: 1次元自由連結鎖の j 番目のボンドベクトルの向きを表すための int型（整数型）の配列変数 l[N]の宣 言，これにより， 変数，l[0] . . . l[N− 1] の整数型の変数を使う準備ができる．また， 末端の位置がある場 所に来た回数を計算するための整数型の配列変数， p[N+1] の宣言をしている．これにより, p[0] . . . p[N], の整数型の変数を使う準備ができた． 17 行目: 整数型の各種カウンターの宣言，および，でたらめにボンドを選んでひっくり返 す際の，ひっくり返す前 の値を格納しておく変数 lbefor と，末端の位 置を計算して格納する変数 end の宣言がしてある． 18 行目: float型（実数型）の変数として，ボンドをひっくり返す際に必要なエネ ルギーを格納する dene, および 温度を格納する変数 t を宣言す る． 19 行目: 結果を書き込むファイルを使う際に必要な, FILE 型のポインタ変数 fp の宣言． 21 ‒ 23 行目: 温度を標準入力（キーボード）から入力し，変数 t に格納する手続き． 25 ‒ 27行目: ボンドの向きを初期化する．for文は，for の後の(初期値設定; ループ実行判断式;ループ付加実行式)に 続く，中括弧 { }の中を繰り返 し実行する．ループ実行判断式が成り立つ場合には中括弧内の作業を実行

し，実 行するごとに，ループ付加実行式が実行される．ここでは， j を最初に 0にし，l[j] = 1 を実行し， その後，j をj + 1にする即ち, j の数を一つ増やして，j < N ならば，再度ループを実行するとい作業を繰 り 返す．最終的に全ての l[j] の値を1にする. 29 ‒ 31行目: jを 0 から N まで変化させながら，p[j] を 0 に初期化 する．p[j] は，末端の位置がどこに来るか，そ の確率を調べるために使 う配列変数である． 33 ‒ 52行目: モンテカルロシミュレーションを実際に行っている部分である．たった20行であ る．最初に mct を 0 に し，33 行目の { から，52 行目の } までの作業を繰り返す．一回作業をするごとに， mct の値を 1 ずつ増加 させ，mct が MCTIME になるまで繰り返す． 34 ‒ 44行目: 次の操作を N 回繰り返す．この N 回が 1 モンテカルロ時間に自由連結鎖 に行う作業である． 35 行目にて，j に 0 から N までの整数をでたらめに代入する．すなわち， でたらめにボンドを選んで， そのボンドの向きが l[j]となる． l[j] は 1 か −1 の値しか取らない．36行目にて，lbefore に，現在の l[j] の値を記録する．37行目にて，l[j] の示す向 きを逆転する．すなわち，l[j] が 1 であればこれを −1 とし， −1 であれば +1 になるように，l[j] に −1 を掛けたものを l[j] に代入する． l[j] *= (-1); は， l[j] = l[j]*(-1); と同じことである． 38行目にて，ボンドを逆転するのに必要なエネルギーを計算する． 逆転後の値を lj とし，逆転前の値 を l′ j とする． 式(6)より，逆転により変化しているのは lj だけであることをこうりょすると，逆転前のエネ ルギー H′_{と， 逆転後のエネルギー H の差は，} H− H′ = −ϵ(l0+· · · + lj+· · · + lN−1) −(−ϵ)(l0+· · · + l′_j+· · · + lN−1) = −ϵlj− (−ϵl′j), (22) である．38行目で dene にこの値を代入している．

deneが正の場合，ボンドの向きの逆転によりエネルギーが増加するため， p.5 の方式に従い，ボルツマ ン因子 exp ( −H− H′ kBT ) , (23) の確率で転移をさせる方針をとる．39 行目で H− H′ _{の値の正負を判断し，正 の場合は，40 行目で [0, 1]} の一様乱数と (23) を比較し，一様乱数の 方が大きければ，41 行目で l[j] を元の値に戻すようにする． 46 ‒ 48行目: 末端の位置を求めている． 末端の位置を x とすると，全ての lj が −1 のとき，x = −N である． このとき，一つのボンドをひっく り返すと，x = N− 2 になる．一つのボンド を反転すると x は 2 変化するので，x = −N, −N + 2, −N + 4, . . . , N−4, N −2, N の値をとる．j = (x+N)/2 とすれば，x = −N, −N +2, −N +4, . . . , N −4, N −2, N に対し，j = 0, 1, 2, . . . , N − 2, N − 1, N となる． 49 行目で，x = 2j− N の時，p[j] に 1 を加えることになる． この操作は，1モンテカルロ時間に1 回行う．すなわち，p[j] には，1モ ンテカルロ時間ごとに，末端がどこにあるか調べ，末端が x ならば, p[(xN)/2]+ に 1 を加える．こうしておくと，後で，末端が， x = 2j− N をどれくらいの 頻度で訪れたか を知ることができる． 51 行目: モンテカルロ時間の表示． 54 ‒ 60 行目: ファイル result.txt に結果を出力する手続きである． 57 行目にて, 末端の位置 2j− N と，末端がそ の位置を取った回数をサンプリ ング数で割ったもの, p[j]/MCTIME が出力してある． 62行目: このメイン関数を終了します．

2DmcGL.c

次のプログラムは2次元自由連結鎖で，エネルギーとして H =−ϵ n_∑−1 j=1 lxj, (24)

を用いた場合のモンテカルロシミュレーションを行うものである． ただし, lxj は, j番目のボンドの x-座標 成分である. 可視化部分は，rothex.c を参考に書きかえた． 1 #include <stdio.h> 2 #include <stdlib.h> 3 #include <GL/glut.h> 4 #include <math.h> 5 #include <time.h> 6 7 #define WIDTH 600 8 #define HIGHT 600 9 #define RADIUS 100 10 #define LENGTH 3 11 #define X0 300 12 #define Y0 300 13 #define DEGSTEP 10 14 15 #define N (100) 16 #define EP (6.0e-22) /* [J] */

17 #define KB (1.380662e-23) /* The Botzmann's constant [J/K] */ 18 #define MCTIME (1000)

19 #define DL (0.5) 20

21 void display(void);

22 void reshape(int w, int h);

23 void keyboard(unsigned char key, int x, int y); 24 void ginit(int* pargc, char** argv);

25 void idle(void); 26

27 int Degshift=0; 28 int Zoom=1,Stop=1;

29 float Ep,T, Lx[N], Ly[N]; 30

32 double uran(){ 33 return((double)rand()/RAND_MAX); 34 } 35 36 37 int mcrun() 38 { 39 int j; 40 float dene;

41 float lxbefore, lybefore; 42 float dx, dy, r; 43 44 j = (N-1)*uran()+1; 45 lxbefore = Lx[j]; 46 lybefore = Ly[j]; 47 48 do{ 49 dx = 2.0*DL*(uran()-0.5); 50 dy = 2.0*DL*(uran()-0.5); 51 }while( dx*dx + dy*dy > DL*DL); 52 53 /* 54 Lx[j] = Lx[j] + dx; 55 Ly[j] = Ly[j] + dy;

56 */ 57 58 Lx[j] = dx; 59 Ly[j] = dy; 60 61 r = sqrt(Lx[j]*Lx[j] + Ly[j]*Ly[j]); 62 Lx[j] = Lx[j]/r; 63 Ly[j] = Ly[j]/r; 64

66 dene = dene*EP; 67 68 if(dene > 0){ 69 if(exp(-dene/(KB*T)) < uran()){ 70 Lx[j] = lxbefore; 71 Ly[j] = lybefore; 72 } 73 } 74 75 return(0); 76 } 77 78 void idle(void) 79 { 80 if(Stop == 1)return; 81 mcrun(); 82 glutPostRedisplay(); 83 } 84 85 void display(void) 86 { 87 int x, y, j; 88 89 glClear(GL_COLOR_BUFFER_BIT); 90 glPointSize(1); 91 glBegin(GL_LINE_STRIP); 92 glColor3f(1.0, 1.0, 1.0); /* white */ 93 x = 0; 94 y = 0; 95 glVertex2f(X0+x, Y0+y); 96 for(j=1;j<N;j++){ 97 if(j%2 == 0){ 98 glColor3f(1.0, 1.0, 1.0); 99 }else{

100 glColor3f(1.0, 1.0, 0.0); 101 } 102 x = x+Zoom*LENGTH*Lx[j]; 103 y = y+Zoom*LENGTH*Ly[j]; 104 glVertex2f(X0+x, Y0+y); 105 } 106 107 glEnd(); 108 glutSwapBuffers(); 109 } 110

111 void reshape(int w, int h) 112 {

113 glViewport(0, 0, (GLsizei)w, (GLsizei)h ); 114 glMatrixMode(GL_PROJECTION);

115 glLoadIdentity();

116 gluOrtho2D(0.0, (GLdouble)w, 0.0, (GLdouble) h); 117 }

118

119 void keyboard(unsigned char key, int x, int y) 120 { 121 switch(key) { 122 case 's': 123 Stop *= -1; 124 break; 125 case 'i': 126 if(Zoom == 0){ 127 Zoom=1; 128 }else{ 129 Zoom *= 2; 130 } 131 break; 132 case 'o': 133 Zoom /= 2;

134 break; 135 case 'e': 136 printf("\nElectric field ="); 137 scanf("%f",&Ep); 138 Stop = 1; 139 break; 140 case 't': 141 printf("\nTemperature ="); 142 scanf("%f",&T); 143 Stop = 1; 144 break; 145 case 'q': 146 case 'Q': 147 exit(0); 148 default: 149 break; 150 } 151 } 152

153 void ginit(int* pargc, char** argv) 154 {

155 glutInit(pargc,argv);

156 glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB); 157 glutInitWindowSize(WIDTH,HIGHT);

158 glutInitWindowPosition(100,50);

159 glutCreateWindow("Numerical Recipes for Polymer Science"); 160 glClearColor(0.0,0.0,0.0,0.0); /* background color */ 161 glShadeModel(GL_FLAT); 162 glutDisplayFunc(display); 163 glutReshapeFunc(reshape); 164 glutKeyboardFunc(keyboard); 165 glutIdleFunc(idle); 166 } 167

168 int main(int argc, char* argv[]) 169 { 170 int j; 171 ginit(&argc,argv); 172 173 printf("Temperature ="); 174 scanf("%f",&T); 175 printf("\n"); 176 printf("Electric field ="); 177 scanf("%f",&Ep); 178 179 for(j=1;j<N;j++){ 180 Lx[j] = 0.0; 181 Ly[j] = 1.0; 182 } 183 184 glutMainLoop(); 185 return 0; 186 } 187