1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80
| program md
call initialize(x,y,z,vx,vy,vz)
print*,x(1)
end program md
subroutine initialize(x,y,z,vx,vy,vz)
implicit none
integer n,nn,i,j,k,nb
real*8 s,s5,rho,fn,k0,cu,sigma,rb,temp,xm,ym,zm
parameter(nn=6,n=nn**3,rho=0.8,sigma=1.0,temp=1.35)
real*8, dimension(1,n)::x,y,z,vx,vy,vz
print*, "computing"
S=(real(N)/rho)**(1.0/3.0)
s5=s/2.0
fn=float(n)
k0=0.5*(3.0*fn-4.0)*temp
cu=2.5*sigma
rb=(S/4)**2
! initialize the momentum of the MD box
xm=0.0
ym=0.0
zm=0.0
nb=1
do i=0,NN-1
do j=0,NN-1
do k=0,NN-1
! use a lattice instead of random for positions
x(nb)=(0.5+real(i))*S/real(nn)
y(nb)=(0.5+real(j))*S/real(nn)
z(nb)=(0.5+real(k))*S/real(nn)
vx(nb)=(-1)**nb*(nb+k)
vy(nb)=(-1)**nb*(nb+i)
vz(nb)=(-1)**(nb+1)*(nb+j)
xm=xm+vx(nb)/fn
ym=ym+vy(nb)/fn
zm=zm+vz(nb)/fn
nb=nb+1
enddo
enddo
enddo
nb=nb-1
return
end subroutine initialize |
Partager