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 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196
| !****************************************************************************
!
! =================MD for the Lennard-jones particules=======================
!
! -------------------------Input Parametre as follows------------------------
!
!
! n Number of particles (must bet a multiple of 4)
! x,y,z The positions
! vx,vy,vz Velocities
! T Reduced Temperature
! rho Densite
! cut Cutoff of the potential in sigma units
! dt Basic time step
! npasos Number of integration steps
! aL Side length of the cubical Boxe in sigma units
! dum Is the seed
! ngr Is the size of the histogram
! hgr Is the width of bin
!
!
!****************************************************************************
implicit real*8(a-h,o-z)
parameter (n=256)
dimension x(n),y(n),z(n)
dimension vx(n),vy(n),vz(n)
dimension fx(n),fy(n),fz(n)
dimension hist(1000),gdr(1000)
real*4 r
common // hgr,cutr2,ngr
!****************************************************************************
!
!==================Definition of the simulation parametrs====================
!
!****************************************************************************
rho = 0.83134
T = 0.722
npasos = 2000
aL = 6.76284
cut = 2.5
pi=4d0*datan(1d0)
!
call init(n,x,y,z,aL,aL2)
!****************************************************************************
aL2 = aL/2d0
aL4 = aL/4d0
! cut2 = (2.5)**2
! cutr2= (aL/2)**2
ngr = 1000
hgr = aL2/ngr
! dt = 0.01d0
!
!
do i=1,1000
hist(i)=0
end do
!
! call grr (ngr,x,y,z,aL,aL,aL,aL2,aL2,aL2)
!****************************************************************************
! ............................Histogram......................................
!****************************************************************************
data hr/0.05d0/
do i=1,n-1
do j=i+1,n
xx=x(i)-x(j)
yy=y(i)-y(j)
zz=z(i)-z(j)
call mic (xx,yy,zz,xlx,yly,zlz,xll,yll,zll)
r=xx*xx+yy*yy+zz*zz
rr=sqrt(r)
if(rr.lt.xll) then
bin=int(rr/hr) + 1
if(bin.lt.ngr) then
hist(bin)=hist(bin)+2
endif
endif
enddo
enddo
!****************************************************************************
!normalise radial distribution function
!****************************************************************************
const=4.0*pi*rho/3.0
do bin=1,ngr
rg=(bin-1)*hgr+hgr/2d0
dvol=const*((rg+hgr)**3-(rg)**3)
gdr(bin)=real(hist(bin))/real(npasos)/real(n)/dvol
write(*,*)'distance=', rg
write(*,*)'radial distribution function=',gdr(bin)
! write(*,*) 'histogram=',hist(bin)
end do
!****************************************************************************
! open files
!****************************************************************************
! OPEN (1,file = 'r.dat',form='unformatted')
! OPEN (2,file = 'r,gdr.dat',form='unformatted')
! OPEN (3,file = 'temp.dat',form='unformatted')
! OPEN (4,file = 'u.dat',form='unformatted')
! OPEN (8,FILE = 'ap.DAT',form='unformatted')
! write(2,*) r,gdr
! OPEN (10,FILE = 'donnees.DAT')
end
!****************************************************************************
!
!///////////////////////////////THE SUBROUTINES//////////////////////////////
!
!****************************************************************************
!
subroutine init(lmn,x,y,z,aL,aL4)
!****************************************************************************
!
! ========================Generation of fcc lattice =========================
!
!****************************************************************************
implicit real*8 (a-h,o-z)
parameter (num=256)
dimension x(num),y(num),z(num)
m=0
do lg=0,1
do i=0,3
do j=0,3
do k=0,3
m=m+1
x(m) = i*aL4+lg*aL4*0.5d0
y(m) = j*aL4+lg*aL4*0.5d0
z(m) = k*aL4
end do
end do
end do
end do
do lg=1,2
do i=0,3
do j=0,3
do k=0,3
m=m+1
x(m) = i*aL4+(2-lg)*aL4*0.5d0
y(m) = j*aL4+(lg-1)*aL4*0.5d0
z(m) = k*aL4+aL4*0.5d0
end do
end do
enddo
end do
end
subroutine mic (xx,yy,zz,alx,aly,alz,alx2,aly2,alz2)
!****************************************************************************
!
! =========================Minimum image convention =========================
!
!****************************************************************************
implicit real*8 (a-h,o-z)
if(abs(dble(xx)).ge.alx2)then
if(dble(xx).lt.alx2)then
xx=xx+alx
else
xx=xx-alx
endif
endif
if(abs(dble(yy)).ge.aly2)then
if(dble(yy).lt.aly2)then
yy=yy+aly
else
yy=yy-aly
endif
endif
if(dabs(dble(zz)).ge.alz2)then
if(dble(zz).lt.alz2)then
zz=zz+alz
else
zz=zz-alz
endif
endif
end |
Partager