- #1

svartak

- 5

- 0

i'm a beginner in programming. i wrote a program to calculate velocity correlation function. i do not get any compilation errors but the lines after line number 42 are not being read. could someone tell me what is it that I'm doing wrong?

thanks a lot!

svartak

!program to calculate velocity of all atoms at all iterations in the given trajectory and calculation of velocity correlation function

!velocity at t=0 is written as v(i,2) i.e. velocity calculated at iteration 2

parameter natm=144, idim=1000

real*8 x(natm),y(natm),z(natm),x_old(natm),y_old(natm),z_old(natm)

real*8 vx(natm,idim),vy(natm,idim),vz(natm,idim),dt,proj

real*8 vdot(natm,idim),vmag0, Niteration

character*80 filename

character*2 ch

c call getarg(1, filename)

open(1, file='traj.xyz', status='OLD')

open(2, file='velocity-all.dat', status='unknown')

open(3, file='input-data', status='unknown')

read(3,*)dt !timestep

read(3,*)Nsample !highest iteration number for sampling region

read(3,*)Nmax !last iteration number in the trajectory

read(3,*)dnsample !interval for the selection of sampling points

read(3,*)Nstep !stepsize for interval between v(t0) and v(t)

j=0

do

read(1,*)Natom

read(1,*)

j=j+1 !loop for all iterations

if (j.gt.idim) then

write(*,*)'idim is too small.'

stop

end if

do i = 1, Natom !loop for all atoms in each iteration

read(1,*,end=100)ch, x(i), y(i),z(i) !i'th atom

if (j.ge.2) then

vx(i,j) = (x(i) - x_old(i))/dt !velocity will be in angstrom/fs

vy(i,j) = (y(i) - y_old(i))/dt

vz(i,j) = (z(i) - z_old(i))/dt

end if

x_old(i)=x(i)

y_old(i)=y(i)

z_old(i)=z(i)

write(2,*)vx(i,j),vy(i,j),vz(i,j)

end do

end do

proj=0

Niteration=0

do dN = Nsample+1, Nmax-Nsample, Nstep !dN=v(t)-v(t0)

do k = 1, Nsample, dnsample !loop over different initial velocity points

write(2,*)Niteration

do i = 1, Natom !loop over all atoms

vdot(i,k)=vx(i,k)*vx(i,k+dN)+vy(i,k)*vy(i,k+dN)

& +vz(i,k)*vz(i,k+dN)

vmag0=vx(i,k)*vx(i,k)+vy(i,k)*vy(i,k)+vz(i,k)*vz(i,k)

proj = (proj + (vdot(i,k)/vmag0))

write(2,*)proj

write(2,*)vdot(i,k)

write(2,*)vmag0

end do

end do

end do

vcf = proj/Natom/Niteration

write(2,*)vcf

c end do

100 end