.. index:: molecular dynamics

========================
多体系２
========================

.. toctree::
   :maxdepth: 1

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

数が多くなると
===========================

天体の計算の時にも，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

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

.. code-block:: fortran
    :linenos:
    :emphasize-lines: 25,26

        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というモジュールの中で，
以下のように宣言されている．

.. code-block:: fortran
    :linenos:
    :emphasize-lines: 18

    !+++++++++++++++++++++++++++++++++++++++++++++++
    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行目のように

.. code-block:: fortran

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

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

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

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

.. figure:: ../img/345MD_164sphere.png
   :scale: 70 %
   :align: center
   :alt: 164sphere

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

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

.. figure:: ../img/346MD_164ene.PNG
   :scale: 100 %
   :align: center
   :alt: Energy

動画をYoutubeに掲載した．
`>>Youtube4 <http://youtu.be/lbTezjSLx_Q>`_

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

------------------------------------------

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