ファイルfcode-bungy.f
c2--567--1--------2--------3---------4---------5---------6---------772
c
program fcodebungy
c
c copyright Tatsuki Oda (2009.05.01)
c
implicit double precision (a-h,o-z)
implicit integer (i-n)
c
parameter (nx=1000)
c
integer number,i
integer method
dimension a1(nx),a2(nx),a3(nx)
dimension a(nx)
dimension s(nx)
dimension z(nx),v(nx),f(nx)
c
c
c2--567--1--------2--------3---------4---------5---------6---------772
c
c for file-open
c for input
open(unit=10,file='f-data.in')
c for output
open(unit=17,file='bungy.xy')
c
c
c2--567--1--------2--------3---------4---------5---------6---------772
c
c
c for data input
c
c input-data with (14)
c
read(10,*) nstep ! number of steps
read(10,*) dt ! Delta t
read(10,*) am ! mass of particle
read(10,*) ak0 ! force constant (k) of spring
read(10,*) alength ! length of spring
read(10,*) height ! height of stage
read(10,*) gamma ! friction of resistance
read(10,*) method ! algorithm of solving
c method =0 velret
c =1 Runge-Kutta (2th)
c =2 Runge-Kutta (4th)
c
g=9.80d0 ! acceleratio of gravity
c
c
c2--567--1--------2--------3---------4---------5---------6---------772
c
c
c for initial conditions
c
z(1)=height
v(1)=0.0d0 ! free diving
f(1)=-am*g
c
hml=height-alength
c
c for solving the equation of motion
c
if(method.eq.0)then
z(2)=z(1)+v(1)*dt+f(1)*dt**2*0.5d0/am
v(2)=v(1)+f(1)*dt/am
do i=2,nstep
f(i)=-am*g
eta=gamma*dt/(2.0d0*am)
etap1=eta+1.0d0
z(i+1)=2.0d0*z(i)/etap1-
& z(i-1)*(1.0d0-eta)/etap1+
& f(i)*dt**2/am/etap1
v(i)=(z(i+1)-z(i-1))/(2.0d0*dt)
enddo
else
write(*,*)' Runge-Kutta (2nd or 4th) not installed yet ! '
stop
endif
c
c
c for output of solutions
c
write(17,800)(i,dt*dble(i-1),z(i),v(i),f(i),i=1,nstep)
800 format(i5,4e22.12)
999 continue
stop
end
ファイルf-data.in
500 nstep ! number of steps 0.1 dt ! Delta t 100.0 am ! mass of particles (MKS unit) 80.0d0 ak0 ! force constant (k) of spring 40d0 alength ! length of spring 100d0 height ! height of stage 20d0 gamma ! friction of resistance 0 method ! algorithm (method=0:velret,1:R-K(2th),2:R-K(4th)) i
ファイルmdcode.f
c2--567--1---------2---------3---------4---------5---------6---------772
c
program mdcode
c
c copyright Tatsuki Oda (2006.10.09)
c
implicit double precision (a-h,o-z)
implicit integer (i-n)
c
parameter (natomx=1000)
c
integer number,nsum,i
dimension a1(natomx,3),a2(natomx,3),a3(natomx,3)
dimension b1(natomx,3),b2(natomx,3),b3(natomx,3)
dimension f(natomx,3)
c
c
open(unit=14,file='mddata.in')
open(unit=15,file='xyzdata.in')
open(unit=16,file='data-1.out')
open(unit=17,file='data.xyz')
c
c
c for input
c
c input-data with (14)
c
read(14,*) nstep ! number of steps
read(14,*) dt ! Delta t
read(14,*) am ! mass of particles
read(14,*) ak0 ! force constant of spring
read(14,*) r0 ! natural length of spring
c
c
c
c input-data with (15): initial coordinates of atoms
c
read(15,*) natom ! number of atoms
if(natom.gt.natomx)then
write(06,*)' error #### natom > natomx ####'
stop
endif
c
read(15,*)(a1(i,1),a1(i,2),a1(i,3),i=1,natom)
c
c
c set up the initial conditions
c
do i=1,natom
c
c forces deduced from the result of homeworks
c
f(i,1)=0.0d0
f(i,2)=0.0d0
f(i,3)=0.0d0
enddo
c
do i=1,natom
c
c we may assume no velocity at the initial state
c
b1(i,1)=0.0d0
b1(i,2)=0.0d0
b1(i,3)=0.0d0
enddo
c
c
do i=1,natom
a2(i,1)=a1(i,1)+b1(i,1)*dt+f(i,1)*dt**2/am*0.5d0
a2(i,2)=a1(i,2)+b1(i,2)*dt+f(i,2)*dt**2/am*0.5d0
a2(i,3)=a1(i,3)+b1(i,3)*dt+f(i,3)*dt**2/am*0.5d0
enddo
c
c
c starting a simulation
c
do istep=2,nstep
c
do i=1,natom
c
c forces deduced from the result of homeworks
c
f(i,1)=0.0d0
f(i,2)=0.0d0
f(i,3)=0.0d0
enddo
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
c
c monitering
c
c kinetic energy
c
ek=0.0d0
c
c potential energy
c
ev=0.0d0
c
c total energy
c
et=ek+ev
c
c
c distance of 2 particles
c
r12=sqrt((a2(i,1)-a2(i,1))**2+
& (a2(i,1)-a2(i,1))**2+
& (a2(i,1)-a2(i,1))**2)
c
c writing data of the energies:(16)
c
write(16,500)istep,dt*istep,ek,ev,et,r12
500 format(i7,5e15.6)
c
c
c
c
c treatment for the next step
c
do i=1,natom
a1(i,1)=a2(i,1)
a1(i,2)=a2(i,2)
a1(i,3)=a2(i,3)
a2(i,1)=a3(i,1)
a2(i,2)=a3(i,2)
a2(i,3)=a3(i,3)
enddo
c
c
enddo
c
c
stop
end
ファイルmddata.in
500 nstep ! number of steps 0.01 dt ! Delta t 100.0 am ! mass of particles 85.0d0 ak0 ! force constant of spring 40d0 r0 ! natural length of spring
ファイルxyzdata.in
2 0.0 0.0 0.0 50.0 0.0 0.0