多体系２

　12gの炭素に含まれている炭素原子の数は約6.02×10²³ 個である． 分子動力学法で本格的に何か現実の系を模擬して計算しようとする場合，これに準ずる数の，というわけにはいかないのだが， それなりに多くの粒子で計算を行うことになる．

数が多くなると

天体の計算の時にも，5つの粒子で計算を行うために配列を使うことにしたが， より数が増えていく場合には位置や速度の初期条件をプログラムコード中に書くと，それだけで プログラムの大部分を占めることになり非常にみづらい．そもそも一つ一つの位置，速度を数値で与えるのも大変である．

初期条件は別プログラムで生成する．初期条件を記述したファイルを， ＭＤプログラムで読み込むことにする． たとえば以下のようなファイルを用意したとする．

ndim          3
nbdy        164
mass  1.0000
mass  1.0000
mass  1.0000
.
.
.
.
cord -0.8889 -2.3094 -4.8990
cord  1.1111 -2.3094 -4.8990
cord -1.8889 -0.5774 -4.8990
cord  0.1111 -0.5774 -4.8990
.
.
.
end

これを読み込むために，以下のようなサブルーチンを用意した．

    subroutine read_ini
        use var_trajectory
        use param_particle

        integer,parameter::     unit_ini = 10
        integer                 nline
        integer                 nm,nc
        integer                 i,rstart,rend
        double precision        xx
        character(len=80)       inline

        open(unit=unit_ini,file='ini.dat',status='old')
        do
            read(unit_ini,'(a80)') inline
            if (inline(1:3) == 'end') then
                exit
            else if (inline(1:4) == 'ndim') then
                read(inline(7:16),'(i10)')ndim
            else if (inline(1:4) == 'nbdy')then
                read(inline(7:16),'(i10)')nbody
            end if
        end do
        close(unit_ini)

        allocate (mass(nbody))
        allocate (coord(ndim,nbody))

        open(unit=unit_ini,file='ini.dat',status='old')
        nm = 1
        nc = 1
        do
            read(unit_ini,'(a80)') inline
            if (inline(1:3) == 'end') then
                exit
            else if (inline(1:4) == 'mass') then
                read(inline(6:12),'(f7.4)') mass(nm)
                nm = nm + 1
            else if (inline(1:4) == 'cord') then
                do i = 1,ndim
                    rstart = 6 + 8*(i-1)
                    rend = rstart + 6
                    read(inline(rstart:rend),'(f7.4)') coord(i,nc)
                end do
                nc = nc + 1
            end if
        end do

        close(unit_ini)
    end subroutine read_ini

初期条件ファイルには，最初の4文字で後に続く内容の区別をしている． サブルーチンではこの4文字を手がかりにif文判定で情報読み込んでいく． 整数を読み込むのか小数点実数を読み込むのかによって，書式（i10とかf7.4とか）を使い分ける．

ここで一つ問題となるのが配列の確保である． プログラムの中で必要となる配列は，次元や粒子数によって大きさが決まるものが多いが， その次元（ndim），粒子数（nbody）も別ファイルから読み込むことにしてある． これらの読み込みが済むまで配列の大きさが決まらないのである．

割り付け配列を使う必要がある． たとえばcoordという配列はvar_trajevtoryというモジュールの中で， 以下のように宣言されている．

!+++++++++++++++++++++++++++++++++++++++++++++++
module cal_ene

 module globals

double precision,parameter::    pi = acos(-1.0d0)

integer                         ndim
integer                         nbody

end module globals
!+++++++++++++++++++++++++++++++++++++++++++++++
module var_trajectory

use globals
implicit none

double precision,allocatable::  coord(:,:)

end module var_trajectory

allocatable属性というのがついている． このように宣言しておいた上で，read_iniサブルーチン内の25,26行目のように

allocate (mass(nbody))
allocate (coord(ndim,nbody))

nbody,ndimの値を先に読み込んでからmass,coordをallocate（’割り当てる’の意味）して，配列の大きさを確定している．

以下は，上記の方法で行った計算の例である．

164個のアルゴン原子をＬＪポテンシャルのvan der Waals半径(\(R=2^{1/6}\sigma\) )に基づいて面心立方状に 配置し10⁷ ステップ計算した．

ちなみに，私の普段使いノートPCで約2時間半も時間がかかった．

やはり数が多くても全エネルギー（運動エネルギー＋ポテンシャルエネルギー）は一定である．

動画をYoutubeに掲載した． >>Youtube4

速度の大小によって色をつけている． ポテンシャルの谷底の間隔で結晶配置をつくったので，その構造を維持したまま振動をしているだけである．

---

さて，この先本格的に分子動力学法計算を進めていくためには， 計算時間を短縮する工夫や境界条件の導入， 物理量の計算，温度や圧力の制御などが必要である． 分子動力学法は，ここまでに述べた本体部分よりも， この先に必要となるそれ以外の部分がよほど難しく手間がかかると思う．