初期配置の設定

面心立方格子

シミュレーションの初期配置としては結晶構造がよく用いられます。 基本的にシミュレーションの結果は初期配置に依存しません。 従って、どのような初期配置からはじめてもよいのですが、 原子の間の距離が極端に近い場合は、 想定以上の強い斥力が働くのでシミュレーションが正しく動かなくなります。 初期配置として結晶構造を使用する場合はそのような心配はなくなります。

アルゴンの結晶構造は面心立方格子(Face Centered Cubic lattice;FCC) ですので、初期配置としてはこの構造を使用します。

面心立方格子

周期境界条件

分子シミュレーションではよく境界条件として周期境界条件を使用します。 シミュレーションを行う空間が立方体や直方体のセルを使用しますが、 場合周期境界条件ではそのセルの上下、前後、左右にセルのコピーを配置します。 こうすると、 例えばセルの下側に粒子が出て行った場合はその粒子はセルの上側から入ってきます。 実際にシミュレーションを行うセルを基本セル、 コピーしたセルをイメージセルと呼びます。

周期境界条件

粒子の配置

周期境界条件の下で粒子が面心立方格子を作るように配置する方法を考えてみます。 面心立方格子の単位格子に入っている実際の粒子の数は4つですので、 下の左の図のように粒子を配置すれば、 周期境界条件下では下の右の図のようになります。

実際の粒子配置 周期境界を考えた粒子配置

には4 るのでρ  とするとρ  きさl
    4
ρ = -3
    l
(1)

けるのでl
    ∘ --
l =  34-
      ρ
(2)

となりますって4 つのlh = h∕2  いて

(0,0, 0)
(l,l ,0)
 h  h
(0,lh, lh)

(lh,0, lh)
けます

x, y,z  をそれぞれNl  とするとNp
Np =  4N 3l
(3)

シミュレションセルL
L =  lNl
(4)

ります

プログラミング

面心立方格子の単位格子の内の粒子の数をFCC_PARTILCE_NUM、 面心立方格子の単位格子の数 Nl をlattice_size、全粒子数 Np  をpnに設定します。

enum {
  FCC_PARTICLE_NUM = 4
};
int lattice_size = 2;
int pn;
pn = FCC_PARTICLE_NUM*lattice_size*lattice_size*lattice_size;

enum は列挙型の変数を定義するためのもので、 整数の定数を設定するのによく使用されます。 上記の例では、整数型でlattice_sizeとpnを定義し、lattice_sizeに2を代入し、 粒子数pnを式(3)に従って計算しています。

単位格子の一辺の長さl をside_unit、シミュレーションセルの一辺の長さL をsideとすると、 それぞれの計算式は以下のようになります。

side_unit = pow(FCC_PARTICLE_NUM / density, 1./3.);
side = side_unit*lattice_size;

densityには数密度が入っています。

演習
前回の演習の結果を利用し て、上の式で計算される、side_unitとsideを表示するプログラムを作成せよ。

粒子の座標を入れる変数には実数型(double)の配列変数cdを使用します。 配列の大きさは全粒子数pnの3倍にして、x,y,z座標を原子ごとに保存します。

例として、1次元状に並んだ粒子を作ってみます。

#include <stdio.h>
main()
{
  enum {
    N = 4
  };
  int i;
  double cd[N*3];
  double l = 1.0;
  for(i = 0; i < N; i++){
    cd[i*3]   = l*i;
    cd[i*3+1] = 0;
    cd[i*3+2] = 0;
    printf("%d %f %f %f\n",i,cd[i*3],cd[i*3+1],cd[i*3+2]);
  }
}

上記のプログラムをコンパイルし実行すると以下のような結果になります。

0 0.000000 0.000000 0.000000
1 1.000000 0.000000 0.000000
2 2.000000 0.000000 0.000000
3 3.000000 0.000000 0.000000

2次元状に粒子を並べるのプログラムは以下のようになります。

#include <stdio.h>
main()
{
  enum {
    N = 4
  };
  int i,j,k;
  double cd[N*N*3];
  double l = 1.0;
  for(i = 0; i < N; i++){
    for(j = 0; j < N; j++){
      k = (i*N+j);
      cd[k*3]   = l*i;
      cd[k*3+1] = l*j;
      cd[k*3+2] = 0;
      printf("%d %f %f %f\n",k,cd[k*3],cd[k*3+1],cd[k*3+2]);
    }
  }
}

上記のプログラムを実行すると以下のような結果となります。

0 0.000000 0.000000 0.000000
1 0.000000 1.000000 0.000000
2 0.000000 2.000000 0.000000
3 0.000000 3.000000 0.000000
4 1.000000 0.000000 0.000000
5 1.000000 1.000000 0.000000
6 1.000000 2.000000 0.000000
7 1.000000 3.000000 0.000000
8 2.000000 0.000000 0.000000
9 2.000000 1.000000 0.000000
10 2.000000 2.000000 0.000000
11 2.000000 3.000000 0.000000
12 3.000000 0.000000 0.000000
13 3.000000 1.000000 0.000000
14 3.000000 2.000000 0.000000
15 3.000000 3.000000 0.000000

演習
粒子を3次元状に並べるプログラムを作成せよ。(単純立方格子)



次に、粒子の配置が面心立方格子となるようなプログラムを考えます。粒子の配置にあるように面心立方格子では、 単純立方格子の格子点に4つずつ粒子を配置すればよいので、 粒子の座標を表示するプログラムは3重のfor文を使用して以下のようになります。

i = 0;
for(iz = 0; iz < lattice_size; iz++){
  for(iy = 0; iy < lattice_size; iy++){
    for(ix = 0; ix < lattice_size; ix++){
      cd[i]   = side_unit*ix;
      cd[i+1] = side_unit*iy;
      cd[i+2] = side_unit*iz;
      i += 3;
      cd[i]   = side_unit*(ix+0.5);
      cd[i+1] = side_unit*(iy+0.5);
      cd[i+2] = side_unit*iz;
      i += 3;
      cd[i]   = side_unit*(ix+0.5);
      cd[i+1] = side_unit*iy;
      cd[i+2] = side_unit*(iz+0.5);
      i += 3;
      cd[i]   = side_unit*ix;
      cd[i+1] = side_unit*(iy+0.5);
      cd[i+2] = side_unit*(iz+0.5);
      i += 3;
    }
  }
}

演習
上記のサンプルと前の結果を利用して、 実際のアルゴンで面心立方格子状の結晶を作成せよ。

argon_fcc.c

粒子の配置の表示

粒子配置を表示するプログラムcdviewを利用して、実際に 粒子の配置を確認してみましょう。先ず、表示をcdviewが表示できる形式に変更 します。以下はプログラムの粒子を表示する部分です。

printf("' box_sx=0 box_sy=0 box_sz=0 box_ex=%f box_ey=%f box_ez=%f\n",
       side,side,side);
printf("' box_wt=0.007 r1=0.5\n");
for(i = 0; i < pn; i++){
  printf("%d 1 %f %f %f\n",i,cd[i*3],cd[i*3+1],cd[i*3+2]);
}
プログラム例

演習で作成したプログラムの粒子座標を表示する部分を上記のものに変更し、実 行すると以下のような結果が得られます。

' box_sx=0 box_sy=0 box_sz=0 box_ex=3.395154 box_ey=3.395154 box_ez=3.395154
' box_wt=0.007 r1=0.5
0 1 0.000000 0.000000 0.000000
1 1 0.848789 0.848789 0.000000
2 1 0.848789 0.000000 0.848789
3 1 0.000000 0.848789 0.848789
[中略]
29 1 2.546366 2.546366 1.697577
30 1 2.546366 1.697577 2.546366
31 1 1.697577 2.546366 2.546366

この結果は粒子座標表示プログラムcdview を使って表示できる形式になっています。

1行目はシミュレーションセルを表示するための情報で box_ex,box_ey,box_ezはシミュレーションセルサイズに設定されています。 2行目以降は粒子番号、粒子種、x座標、y座標、z座標となっています。 粒子種は原子を表示するときに色を設定するもので1の場合は緑になります。

Qt Creatorを使ってプログラムを作っている場合には、 「アプリケーション出力」 に表示された粒子の位置情報をコピー&ペーストでメモ帳などに貼り付けて、 例えばファイル名fcc.cdvでファイルに保存します。

コマンドラインでのプログラムを作っている場合、 画面に出力された結果をファイルにするには、'>'(リダイレクト) を使用します。例えば、実行プログラムがargon.exeの場合

 argon > fcc.cdv

とすると、 今まで画面に出力されていた結果がファイルfcc.cdvに出力されるようになります。

cdviewを 「ダウンロードとインストール」 からダウンロードしてして展開し、 cdview3.exeを現在の作業フォルダーにコピーして、 作った座標ファイルをcdview3にドラッグ&ドロップするか、 コマンドラインから以下のようなコマンド

 cdview3 fcc.cdv

を実行すると

cdviewの表示結果

のようにcdviewが粒子を表示してくれます。 マウスやキーボードででいろいろ操作が出来ますので、 詳しいcdviewの使い方はcdviewのページ を参照してください。('p'キーを押すと周期境界条件を反映した表示なります。)