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 |
