初期速度の設定

正規分布を持つ速度

分布分布になるようにめますはどのようなえてもよい のですがでは分布になるのでするために分布 えておきます分布るために Box-Muller いま u0   u1   0 から 1 のとき以下いるとS  分布となりま
    ∘ ---------
S =   - 2log u0cos(2πu1 )
(1)

これで分布られますがシミュレションではしないように ゼロになるようします1N  のときM  α(α = x,y,z )
          N
M   =  1-∑  mv
  α    N i=1   i
(2)

となるので
N∑
   m (vi - M α ) = 0
i=1
(3)

となりvi  vi - M α  いたゼロになります

分布を調べる

先ず、乱数の分布を求めるプログラムを作ってみます。C言語ではrand() で0〜RAND_MAXの間の整数値の乱数を発生させることができます。 整数を実数に変換するキャスト演算子(double)を利用すると、 0〜1の間の乱数を発生できます。

double d;
d = (double)rand()/RAND_MAX;

更にこれをN倍すれば、0〜Nの間の乱数も発生できます。

次に分布を表示するプログラムを作ってみます。 分布を調べるためには適当な幅を決めて、 乱数がその幅に入った数を数えるようにします。

乱数の分布

プログラムにすると以下のようになります。

#include <stdio.h>
#include <stdlib.h>
main()
{
  int i;
  double d;
  int bunpu[5];
  srand(1);
  for(i = 0; i < 5; i++){
    bunpu[i] = 0;
  }
  for(i = 0; i < 10; i++){
    d = (double)rand()/RAND_MAX*5;
    printf("d = %f\n",d);
    if(d < 1.0){
      bunpu[0]++;
    } else if(d < 2.0){
      bunpu[1]++;
    } else if(d < 3.0){
      bunpu[2]++;
    } else if(d < 4.0){
      bunpu[3]++;
    } else {
      bunpu[4]++;
    }
  }
  for(i = 0; i < 5; i++){
    printf("bunpu[%d]=%d\n",i,bunpu[i]);
  }
}

演習

上のプログラムでは分布を調べるためにif文を複数利用していて効率が良くない。 if文を使用せず分布を調べられるようにせよ。

分布の結果を数字で見てもばらつき具合がよくわからないので、 簡単な図を作ることにします。表示部分を以下のようにすると、 グラフのようなものが作れます。

  for(i = 0; i < 5; i++){
    printf("%d: ",i);
    for(j = 0; j < bunpu[i]; j++){
      printf("*");
    }
    printf("\n");
  }  

発生させる乱数を増やしていくと、 徐々に均等な分布になっていくことがわかります。 以下は乱数を10万個発生させたときの分布の例です。

bunpu[0]=20098
bunpu[1]=20072
bunpu[2]=20102
bunpu[3]=19939
bunpu[4]=19789
0: ********************
1: ********************
2: ********************
3: *******************
4: *******************
rand.c rand2.c rand3.c rand4.c rand5.c rand6.c

演習

Box-Muller法を使用して正規分布の乱数を発生させ、 -5.0〜5.0の範囲でのその分布を調べ、図示せよ。
ただし、分布を入れる配列変数の数は20とせよ (int bunpu[5]; -> int bunpu[20];)

プログラム作成のヒント
・ファイルの上方に#include <math.h>を入れるとM_PIにπ(=3.141592...)が入ります。
・ √xはsqrt(x)、log x は log(x)で計算できます。

全運動量をゼロにする

で説明したように、 速度から全運動量の和を粒子数で割った値を引くと、 全運動量をゼロにできます。

プログラムにすると以下のようになります。(pnは粒子数)

  srand(1);
  /* 速度の設定 */
  for(i = 0; i < pn*3; i++){
    u0 = (double)rand()/RAND_MAX;
    u1 = (double)rand()/RAND_MAX;
    vl[i] = sqrt(-2*log(u0))*cos(2*M_PI*u1);
  }

  /* 速度の和の計算 */
  for(i = 0; i < 3; i++){
    sum[i] = 0;
  }
  for(i = 0; i < pn*3; i += 3){
    for(j = 0; j < 3; j++){
      sum[j] += vl[i+j];
    }
  }
  /* 速度の和をゼロになるようにする */
  for(i = 0; i < 3; i++){
    sum[i] /= pn;
  }
  for(i = 0; i < pn*3; i += 3){
    for(j = 0; j < 3; j++){
      vl[i+j] -= sum[j];
    }
  }  

  /* 速度の和がゼロになっているかどうかの確認 */
  for(i = 0; i < 3; i++){
    sum[i] = 0;
  }
  for(i = 0; i < pn*3; i += 3){
    for(j = 0; j < 3; j++){
      sum[j] += vl[i+j];
    }
    printf("%3d % f % f % f\n",i/3,vl[i],vl[i+1],vl[i+2]);
  }
  printf("sum % f % f % f\n",sum[0],sum[1],sum[2]);

rand()はゼロからRAND_MAXの間の乱数を発生させる関数です。 RAND_MAXの値はstdlib.hの中で定義されています。 (double)rand()/RAND_MAXとするとゼロから1の間の乱数が得られます。 (double)はキャスト演算子と呼ばれるもので、 これにより変数をdouble型に変換することができます。 M_PIはmath.hの中で定義されており円周率の値を得ることができます。

演習

上記のプログラムを参考に、 実際にアルゴンのプログラムで正規分布を持つ初期速度を発生させよ。

argon_vel0.c

温度の計算

T  v
 i
3        1 ∑N  1
-kBT  =  ---   -mv2i
2        N i=1 2
(4)

となるので
            N
       -1--∑     2
kBT  = 3N     mv i
           i=1
(5)

となりはこのからめることがますプログラムとなっているの
         N
T′ = -1--∑  v′2
     3N  i=1  i
(6)

できますただしゼロにしているので1 っていますって 使する
               N∑
T′ = ----1----    v′2
     3(N  - 1) i=1  i
(7)

となります

この式を使った温度を計算するプログラムは以下のようになります。

  vl2 = 0;
  for(i = 0; i < pn; i++){
    vl2 += vl[i*3]*vl[i*3] + vl[i*3+1]*vl[i*3+1] + vl[i*3+2]*vl[i*3+2];
  }
  temp =  vl2/(3*(pn-1));

得られた温度とシミュレーションを行うときの設定温度を

  printf("%f %f\n",target_temp,temp);

のように比較すると、結果は

0.722000 1.097547

となり、設定温度と現在の温度は異なっていることがわかります。

演習

速度から温度を計算し上記の結果を得るプログラムを作成せよ。

argon_vel0.c

温度の補正

現在 ′
T  ′
T0   しくするためにはx  じますこのと きの
          N
----1----∑      ′2
3(N  - 1)    (xv i)
          i=1
(8)

となりこれは 2  ′
x T なので
T ′0 = x2T ′
(9)

からx
    ∘ ---
      T0′
x =   T ′
(10)

られます

演習

上記の式を参考に速度を規格化し、 実際に設定温度と等しくなっていることを確認せよ。

速度分布の確認

最後に得られた速度の分布を確認してみます。

void calc_velocity_distribution(double *vl, int pn)
{
  int i,j;

  enum {
    N = 20
  };
  int vel_dis[N];
  double dv;
  double dis_min,dis_max;

  dis_min = -2;
  dis_max =  2;
  dv = (dis_max - dis_min)/N;
  
  for(i = 0; i < N; i++){
    vel_dis[i] = 0;
  }
  for(i = 0; i < pn*3; i++){
    j = (vl[i]-dis_min)/dv;
    if(j >= 0 && j < N){
      vel_dis[j]++;
    }
  }
  for(i = 0; i < N; i++){
    printf("% f %f\n",
           dv*(i+0.5)+dis_min,
	   (double)vel_dis[i]/pn*3);
  }
}

この関数では-2から2の間の速度を20等分して、その分布を求めています。 速度分布をlattice_numが2と20の場合の二通りで求めた結果は

速度分布
md_program3.c

となります。面心立方格子の単位格子の数lattice_numと全粒子数との関係は

lattice_numpn
232
2032000

となります。前粒子数が32の場合の速度分布はばらつきが大きくなっていますが、 32000の場合は速度分布が正規分布になっていることがわかります。