読者です 読者をやめる 読者になる 読者になる

lochtext

http://us.battle.net/sc2/en/profile/2312833/1/lochtext/

FortranをCっぽく書く

授業で分子動力学(※)シミュレーションをFortran実装するという課題が出たのでやってみた。本当はFortran2003の機能を使ってC++っぽくオブジェクト指向で書きたかったが途中で挫折した。

とりあえずCっぽくはなったと思うので公開してみる。
他の人が読めないFortranコードを書く研究者が減るといいなぁ。

(※)FCC構造に置いたアルミニウム原子のシミュレーション。

参考:
Verlet法の計算手順
ドウジンテイスウ.log — Fortran2003の試用ーgfortran

module particles
implicit none

real(8) :: dt = 1.0d-15
real(8) :: m = 1.67d-27 * 27.0d0
real(8) :: a0 = 4.05d-10

integer :: step
integer :: endStep = 2000

type :: particle
  real(8) :: x(3),xnew(3),xold(3) !座標
  real(8) :: v(3) !速度
  real(8) :: f(3) !力
  real(8) :: r(108) !原子間距離(allocatable is better)
end type particle

end module particles

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
program MD
use particles
implicit none

type(particle) :: ps(108)
integer i

call psInit(ps)
call forceUpdate(ps)
call positionUpdateInit(ps)

do step=1,endStep

  call forceUpdate(ps)
  call positionUpdate(ps)
  call velocityUpdate(ps)
  call writeTrajectory(ps)

  call goToNextStep(ps)

end do

stop
end program

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine psInit(ps)
use particles
  type(particle):: ps(108)
  integer i,j,k

  do i =1,3
    do j = 1,3
      do k = 1,3
        ps(((i-1)*9+(j-1)*3+k)*4-3)%x(1)=(i-1)*a0
        ps(((i-1)*9+(j-1)*3+k)*4-3)%x(2)=(j-1)*a0
        ps(((i-1)*9+(j-1)*3+k)*4-3)%x(3)=(k-1)*a0
      end do
    end do
  end do

  do i = 1, 108, 4

  ps(i+1)%x(1)=ps(i)%x(1)+a0/2.0d0
  ps(i+1)%x(2)=ps(i)%x(2)+a0/2.0d0
  ps(i+1)%x(3)=ps(i)%x(3)+0.0d0
  ps(i+1)%v(1)=0.0d0
  ps(i+1)%v(2)=0.0d0
  ps(i+1)%v(3)=0.0d0

  ps(i+2)%x(1)=ps(i)%x(1)+0.0d0
  ps(i+2)%x(2)=ps(i)%x(2)+a0/2.0d0
  ps(i+2)%x(3)=ps(i)%x(3)+a0/2.0d0
  ps(i+2)%v(1)=0.0d0
  ps(i+2)%v(2)=0.0d0
  ps(i+2)%v(3)=0.0d0

  ps(i+3)%x(1)=ps(i)%x(1)+a0/2.0d0
  ps(i+3)%x(2)=ps(i)%x(2)+0.0d0
  ps(i+3)%x(3)=ps(i)%x(3)+a0/2.0d0
  ps(i+3)%v(1)=0.0d0
  ps(i+3)%v(2)=0.0d0
  ps(i+3)%v(3)=0.0d0

  end do

  return
end subroutine psInit

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine positionUpdateInit(ps)
use particles
implicit none
  type(particle) :: ps(108)
  integer :: i,j

  do i = 1,108
    do j = 1,3
      ps(i)%xold(j) = ps(i)%x(j)
      ps(i)%x(j) = ps(i)%x(j)+dt*ps(i)%v(j)+dt*dt/(2.0d0*m)*ps(i)%F(j)
    end do
  end do

  return
end subroutine positionUpdateInit

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine forceUpdate(ps)
use particles
implicit none
  type(particle) :: ps(108)

  real(8) :: eps = 0.27d0 * 1.602d-19
  real(8) :: alpha = 1.16d10
  real(8) :: r0 = 3.25d-10

  integer :: i,j,k

  do i=1,108
    do j = 1,3
      ps(i)%f(j) = 0.0d0
      do k=1,108
        if(i /= k) then
          ps(i)%r(k) = dsqrt((ps(i)%x(1) - ps(k)%x(1))*(ps(i)%x(1) - ps(k)%x(1))+&
                              (ps(i)%x(2) - ps(k)%x(2))*(ps(i)%x(2) - ps(k)%x(2))+&
                              (ps(i)%x(3) - ps(k)%x(3))*(ps(i)%x(3) - ps(k)%x(3)))
          ps(i)%f(j) = ps(i)%f(j) + (-1.0d0*eps)*(-2.0d0*alpha*dexp(-2.0d0*alpha*(ps(i)%r(k)-r0))&
            +2.0d0*alpha*dexp(-1.0d0*alpha*(ps(i)%r(k)-r0)))*(ps(i)%x(j)-ps(k)%x(j))/ps(i)%r(k)
        else
          ps(i)%r(k) = 0.0d0
        end if
      end do
    end do
  end do

  return
end subroutine forceUpdate

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine positionUpdate(ps)
use particles
implicit none
  type(particle) :: ps(108)

  integer i,j

  do i = 1,108
    do j = 1,3
      ps(i)%xnew(j) = 2.0d0*ps(i)%x(j)-ps(i)%xold(j)+dt*dt/m*ps(i)%f(j)
    end do
  end do

  return
end subroutine positionupdate

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine velocityUpdate(ps)
use particles
implicit none
  type(particle) :: ps(108)

!!!!未実装
  
  return
end subroutine velocityUpdate

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine writeTrajectory(ps)
use particles
implicit none
  type(particle) :: ps(108)

  integer :: i
if(mod(step,1000) == 0) then
  print '("step=",i10)',step

  do i = 1,108
    print '("x=",e10.2e2,",y=",e10.2e2,",z=",e10.2e2)',&
	ps(i)%xnew(1),ps(i)%xnew(2),ps(i)%xnew(3)
  end do

end if

  return
end subroutine writeTrajectory

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine goToNextStep(ps)
use particles
implicit none
  type(particle) :: ps(108)

  integer :: i,j

  do i=1,108
    do j=1,3
      ps(i)%xold(j)=ps(i)%x(j)
      ps(i)%x(j)=ps(i)%xnew(j)
    end do
  end do

  return
end subroutine goToNextStep