p208

スピノーダル分解およびブロックコポリマー系の相分離動力学のプログラム

ここまでセル力学系のモデルを作ってきた．100 x 100位の二次元格子(周期境界条件をおく)の上で十分何が起こるかは見られるから，ラップトップで簡単に出来上がったモデルを鑑賞することが出来る．

 Cプログラムは sp.c である．

 以下は同じものにより詳細な説明を付け足したもの.

ただし，このプログラムではf(x)として1.3 tanh x ではなくAxに単純化している．つまり，ユニバーサリティをたてにとって，区分的線形写像にして計算を単純化している．より詳しくいうと，x が (-1/A, 1/A)ではf(x) = Ax, この区間の外では x < 0 ならば f(x) = -1, x > 0ならばf(x) = +1. Aの値は各runごとに聞いてくるからキーボードで入れる．1.2くらいでいい．あまり大きくするとfreezeしてしまう．（以下のプログラムの行のそろえ方はもとのプログラムのように整然とはweb上では見られないようである．本物のプログラムはもう少しきれい)

/* spinodal decomposition and block co-polymer system */

#include <stdio.h>

#define L 50 /* lattice size */格子の大きさ

#define D 0.35 /* diffusion constant */拡散係数

#define R 0.1 /* initial amplitude */初期のオーダパラメタの０のまわりの小さなゆらぎの振幅

#define T 3000 /* time limit */何回iterationするか

#define DT 50 /* time increment */図を画面に出すのは何ステップに一回か．

#define frandom() (float)(((float)random()/(float)( 1 << 30 ))/ 2.0 )

#define max 1. /* display range */

/* A : quench depth 　　　クエンチの設定,1から大きく外れるほどdeep quench

B : block copolymer parameter */ これはブロックコポリマーモデル(4.7.1)のCと同じもの．B > 0 が小さいほどポリマーは長くなる．

char *ctab[] = { "e" , " " , "." , "+" , "*" , "#" , "$" , "E" };　これはディスプレーを文字で速やかに表示する工夫，黒っぽい文字が+1に，空白に近い方が-1になるように設定してある．もちろんもっといい画像出力を使っていいがこれで十分．

main()

{

float x[L+ 2 ][L+ 2 ], xt[L+ 2 ][L+ 2 ], fx[L+ 2 ][L+ 2 ];

float fo1, fo2, Map, Diff, Force, A, B, sum;

double  atof();

int i, j, k;

char buff[ 128 ];

int seed, atoi();

unsigned int c;

 printf( "seed :" );

 gets( buff );

 seed = atoi(buff);

 printf( "A :" ); 原点の周りでのｆの傾き，適当に1より少し大きく，1.1とか1.2を試してみよ．プログラムの前の説明参照

 gets( buff);

 A = atof ( buff );

 printf( "B :" ); ポリマーの長さにあたる．スピノーダル分解では B=0. ブロックコポリマーの典型ならB = 0.01.

 gets( buff);

 B = atof (buff );

 srandom( seed );

/* initial condition */

for (i= 1 ;i<=L;i++)

 {

for (j= 1 ;j<=L;j++)

{

 x[i][j] = R*(frandom()- 0.5 );

}

 }

/* periodic boundary condition imposed */

for (i= 1 ;i<=L;i++)

 {

 x[i][L+ 1 ] = x[i][ 1 ];

 x[i][ 0 ] = x[i][L];

 }

for (i= 0 ;i<=L+ 1 ;i++) 

 {

 x[L+ 1 ][i] = x[ 1 ][i];

 x[ 0 ][i] = x[L][i];

 }

/* updating */

for (k= 0 ;k<=T;k++)

 { 

/*computation of force*/　これはIの計算

for (i= 1 ;i<=L;i++)

{

for (j= 1 ;j<=L;j++)

 {

Diff = x[i+ 1 ][j] + x[i- 1 ][j] + x[i][j+ 1 ] + x[i][j- 1 ];

Diff = 2.0 *Diff + x[i+ 1 ][j+ 1 ] + 

x[i+ 1 ][j- 1 ] + x[i- 1 ][j+ 1 ] + x[i- 1 ][j- 1 ];

Diff = D*(Diff/ 12.0 - x[i][j]); 

fo1 = A*x[i][j];

fo2 = (fo1< max)? fo1:max;

Map = (fo2>-max)? fo2:-max; 

fx[i][j] = Map + Diff - x[i][j];　　　これが I.

 }

}

/* periodic boundary condition for the force */

for (i= 1 ;i<=L;i++)

{

fx[i][L+ 1 ] = fx[i][ 1 ];

fx[i][ 0 ] = fx[i][L];

}

for (j= 0 ;j<=L+ 1 ;j++)

{

fx[L+ 1 ][j] = fx[ 1 ][j];

fx[ 0 ][j] = fx[L][j];

}

/* computation of increment */

for (i= 1 ;i<=L;i++)

{

for (j= 1 ;j<=L;j++)

 {

/* Force is negative actual increment */ (4.5.6)を計算している．

Force = fx[i+ 1 ][j] + fx[i- 1 ][j] + fx[i][j+ 1 ] + fx[i][j- 1 ];

Force = 2.0 *Force + fx[i+ 1 ][j+ 1 ] + 

fx[i+ 1 ][j- 1 ] +fx[i- 1 ][j+ 1 ] + fx[i- 1 ][j- 1 ];

Force = Force/ 12.0 - fx[i][j]; 

xt[i][j] = (1.- B)*x[i][j] - Force;

 }

}

/* periodic boundary condition imposed */

for (i= 1 ;i<=L;i++)

{

xt[i][L+ 1 ] = xt[i][ 1 ];

xt[i][ 0 ] = xt[i][L];

}

for (i= 0 ;i<=L+ 1 ;i++)

{

xt[L+ 1 ][i] = xt[ 1 ][i];

xt[ 0 ][i] = xt[L][i];

}

sum = 0. ;

for (i= 1 ;i<=L;i++)

{

for (j= 1 ;j<=L;j++)

 {

sum = sum + xt[i][j];

 }

}

/* up dating - data swap */

memcpy( x, xt, (L+ 2 )*(L+ 2 )* sizeof ( float ) );

/* display */

if (k%DT== 0 )

{

 printf( "\n\n\ntime= %6d\n" ,k);

 printf( "sum = %6.2f \n" , sum);

for (i= 1 ;i<=L;i++)

 {

if (i% 2 == 0 )

 {

for (j= 1 ;j<=L;j++) 

 {

c = (x[i][j] + max)* 3. /max+ 1. ;

printf( "%s" , ctab[c]);

/* if(x[i][j] >0.)

 {

 printf("$");

 }

 else

 {

 printf(" ");

 } */

 }

 printf( "\n" );

}

 }

for (i= 20 ;i<= 26 ;i++) これは実際のオーダパラメタをサンプルしている．うまくいっている限りいらない．コメントアウトしていい．

 {

if (i% 2 == 0 )

 {

for (j= 48 ;j<=L;j++) 

 {

printf( "%f " , x[i][j]);

 }

 printf( "\n" );

}

 } 

}

 }

}