3.3.1 mdcode.fの若干の説明

ここで新しく登場した重要な事項のみを列挙する。

• （１）2次元の配列宣言

        dimension a1(natomx,3),a2(natomx,3),a3(natomx,3)
  

  大きさnatomxおよび3の大きさの2次元配列a1,a2,a3の３つを使用することを宣言す ると同時にその配列の大きさを指定する。ここでは、その大きさは1000および3に なっており、あらかじめ大きさが確定している必要がある。natomxの代わりに、 1000と書いても同じであるが、このようにしておけば、変更する場合に便利である。

  FORTRANでは、配列a1のデータは、

        a1(1,1)
        a1(2,1)
        a1(3,1)
        a1(4,1)
         .
         .
        a1(natomx,1)
        a1(1,2)
        a1(2,1)
        a1(3,1)
         .
         .
        a1(natomx,natomx)
  

  といった順に並んでいる。通常、a1(natomx+1,1)を参照すると 、a1(1,2)が得られる。

• （２）if文、ifプロック if文は、通常、ifとthenとelseifとelseとendifの組で使用される。 条件によりプログラムの動作を分岐させたいときに使用する。 配布時のmdcode.fでは、変数natomがnatomxよりも大きいときだけ実行する ようにしている。

        if( natom .gt. natomx )then
          write(06,*)' error #### natom > natomx ####'
          stop
        endif
  

• （３）read文

        read(15,*)(a1(i,1),a1(i,2),a1(i,3),i=1,natom)
  

  粒子の初期座標を読み込んでいる。この例は 15番(装置番号15)のファイルから a1(i,1),a1(i,2),a1(i,3)に左から$x,y,z$座標を読み込んでいる。

• （１２）アルゴリズム
  mdcode.fでは、データをファイル(14)および(15)から読み込んで、 分子動力学シミュレーションを行い、結果を出力するプログラムと なっている。出力用の(16)には、エネルギーの時間変化が出力される。 時間積分にはベルレ法が使用される。 以下の部分がベルレ法での粒子位置の更新をする部分である。

  c
  c     coordinates
  c
          do i=1,natom
            a3(i,1)=2.0d0*a2(i,1)-a1(i,1)+f(i,1)*dt**2/am
            a3(i,2)=2.0d0*a2(i,2)-a1(i,2)+f(i,2)*dt**2/am
            a3(i,3)=2.0d0*a2(i,3)-a1(i,3)+f(i,3)*dt**2/am
          enddo
  c
  c     velocities
  c
          do i=1,natom
            b1(i,1)=(a3(i,1)-a1(i,1))/(2.0d0*dt)
            b1(i,2)=(a3(i,2)-a1(i,2))/(2.0d0*dt)
            b1(i,3)=(a3(i,3)-a1(i,3))/(2.0d0*dt)
          enddo
  c
  

  粒子の座標についても出力するのがよいが、可視化のために 新しいソフトウエアが必要なので、この授業では行わない。 XcrysDenと呼ばれるソフトウエアが広く使われている。