配布されたコードは一応、コンパイルと実行が可能となるようにしてあるが、 目的の計算ができるようにはなっていない。まず、コンパイルと実行が可能 であることを確認する。
実行すると計算結果は、出力されたはずであるが、原子(粒子)は運動せず、 エネルギーもゼロとなっている。 これはまだ、粒子に働く力の部分にゼロが指定されているためであり、 モニタリングでエネルギーを計算すべき部分でそれを行っていないからである。
そこでプログラムコードを完成させるために、ヒントを書いたのでそれを 参照して、プログラムコードを完成させる。
do i=1,natom
f(i,1)=0.0d0
f(i,2)=0.0d0
f(i,3)=0.0d0
enddo
の部分である。0.0d0の代わりに計算された値が入るようにすればよい。
各粒子は、周りの粒子から力を受けるので、doループを2重ループに して計算すればよいことが、思いつくであろう。例えば、
do i=1,natom
f(i,1)=0.0d0
f(i,2)=0.0d0
f(i,3)=0.0d0
enddo
do i=1,natom
do j=i+1,natom
f1=(x-component of force) <----この部分は各自考えよ
f2=(y-component of force) <----この部分は各自考えよ
f3=(z-component of force) <----この部分は各自考えよ
f(i,1)=f(i,1)+f1
f(i,2)=f(i,2)+f2
f(i,3)=f(i,3)+f3
f(j,1)=f(j,1)-f1
f(j,2)=f(j,2)-f2
f(j,3)=f(j,3)-f3
enddo
enddo
のようなスタイルのコードを書くことになる。
ただし、(x-component of force)
の部分には、手計算で求めておいた力の式を書くことになる。
ただし、ここでは作用・反作用の法則を顕に用いている。
内側のdo文のカウンターjは、i+1から
始まっていることに注意する。
c kinetic energy
ek=0.0d0
や
c potential energy
ev=0.0d0
といった部分があるが、ekとevに正しい、運動エネルギーと
ポテンシャルエネルギーを代入する。
運動エネルギーは、do文の１重ループを用いて計算できるであろう。
例えば、
ek=0.0d0
do i=1,natom
ek=ek+b1(i,1)**2+b1(i,2)**2+b1(i,3)**2
enddo
ek=ek*0.5d0*am
といったスタイルになるであろう。ポテンシャルエネルギーについては、
do文の２重ループを用いることになる。(各自自分で書いてみよう。)
ただし、同じ粒子間の相互作用を計算しないようにする。
力の計算では、同じ粒子間の組が現れないように内側のループのカウンター
を全ての粒子に対して行わないようにする。if文を使う
手もあるが、計算量が２倍になる。