Bonjour
Je suis débutant en programmation Fortran.
J'ai écrit un programme en langage fortran, de Dynamique Moléculaire, j'ai essayé de déterminer la fonction radiale g(r), mais la partie de l'histogramme me donne des valeurs nulles, pour plus de détaille voici le code.
Merci

Code : Sélectionner tout - Visualiser dans une fenêtre à part
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