MoDeNa  1.0
Software framework facilitating sequential multi-scale modelling
rdf_hs.f90
1 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
2 ! This module ...........
3 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
4 Module rdf_variables
5 
6  implicit none
7  save
8 
9  real, dimension(0:20) :: fac(0:20)
10 
11 End Module rdf_variables
12 
13 
14 
15 
16 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
17 ! SUBROUTINE rdf_int
18 !
19 ! this subroutine calculates the radial pair correlation function of
20 ! hard spheres (mseg=1) or tangent-sphere hard chains (with "mseg"
21 ! spherical segments). The pair correlation functions are for the
22 ! Percus-Yevick closure and the analytic solution was proposed by
23 ! Y. Tang, B.C.-Y. Lu, J. Chem. Phys. 105.1996.8262.
24 ! The density "dense" is the packing fraction, the radius "radius" is
25 ! dimensionless in terms of the hard-sphere diameter. The output is
26 ! the segment-segment radial distribution function g(r).
27 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
28 
29 SUBROUTINE rdf_int ( dense, mseg, radius, rdf )
30 
31  USE rdf_variables
32  IMPLICIT NONE
33 
34  !-----------------------------------------------------------------------------
35  REAL, INTENT(IN) :: dense
36  REAL, INTENT(IN) :: mseg
37  REAL, INTENT(IN) :: radius
38  REAL, INTENT(OUT) :: rdf
39 
40  !-----------------------------------------------------------------------------
41  INTEGER :: n
42  REAL :: rg, omega, radind, cfun, facul
43  !-----------------------------------------------------------------------------
44  ! OPEN (10,file='rdf.xlo')
45 
46  ! determine the faculty of integer numbers from 0 to 20. If this
47  ! subroutine is called repeatedly, it makes sense to do this only once.
48 
49  DO n = 0, 20
50  fac(n) = facul(n)
51  END DO
52 
53  omega = ( mseg - 1.0 ) / mseg
54  rg = 0.0
55 
56  DO n = 0,4
57 
58  radind = radius - 1.0 - REAL(n)
59  rg = rg + (-12.0*dense)**n *( (1.0+dense/2.0-omega*(1.0-dense)) &
60  *cfun(dense,omega,2,n,(n+1),radind) &
61  +(1.0+dense*2.0-(1.0+dense/2.0)*omega +(1.0-dense)/2.0*omega**2) &
62  *cfun(dense,omega,1,n,(n+1),radind) &
63  +((1.0+2.0*dense)*omega-(1.0-dense)*omega**2 ) &
64  *cfun(dense,omega,0,n,(n+1),radind) &
65  -((1.0+dense/2.0)*omega-(1.0-dense)/2.0*omega**2 ) &
66  *cfun(dense,omega,1,n,(n+1),(radind-1.0)) &
67  -((1.0+2.0*dense)*omega-(1.0-dense)*omega**2 ) &
68  *cfun(dense,omega,0,n,(n+1),(radind-1.0)) )
69 
70  END DO
71 
72  rdf = rg / radius
73  ! write (*,*) 'r, g(r)',radius,rg/radius
74 
75 END SUBROUTINE rdf_int
76 
77 
78 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
79 !
80 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
81 
82 REAL FUNCTION lnfac(xx)
83 
84  IMPLICIT NONE
85 
86  !-----------------------------------------------------------------------------
87  REAL, INTENT(IN) :: xx
88  !-----------------------------------------------------------------------------
89  INTEGER :: j
90  REAL, SAVE :: cof(6)
91  REAL, SAVE :: stp
92  REAL :: ser,tmp,x,y
93  !-----------------------------------------------------------------------------
94  DATA cof,stp/76.18009172947146,-86.50532032941677, &
95  24.01409824083091,-1.231739572450155,.1208650973866179e-2, &
96  -.5395239384953e-5,2.5066282746310005/
97 
98  x = xx
99  y = x
100 
101  tmp = x + 5.5
102  tmp = ( x + 0.5 ) * log( tmp ) - tmp
103  ser = 1.000000000190015
104  DO j = 1, 6
105  y = y + 1.0
106  ser = ser + cof(j) / y
107  END DO
108  lnfac = tmp + log( stp * ser / x )
109 
110 END FUNCTION lnfac
111 
112 
113 
114 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
115 !
116 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
117 
118 REAL FUNCTION facul(n)
119 
120  IMPLICIT NONE
121 
122  !-----------------------------------------------------------------------------
123  INTEGER, INTENT(IN) :: n
124 
125  !-----------------------------------------------------------------------------
126  INTEGER :: j
127  INTEGER, SAVE :: ntop
128  REAL, SAVE :: table(165)
129  REAL :: lnfac,mult
130  !-----------------------------------------------------------------------------
131  DATA ntop,table(1) / 0, 1.0 /
132 
133  IF ( n < 0 ) THEN
134  write (*,*) ' Negative factorial in FUNCTION fac'
135  stop
136  ELSE IF ( n <= ntop ) THEN
137  facul = table(n+1)
138  ELSE IF ( n <= 165 ) THEN
139  DO j = ntop+1, n
140  mult = REAL(j)
141  table(j+1)=mult*table(j)
142  END DO
143  ntop = n
144  facul = table(n+1)
145  ELSE
146  facul = log( lnfac( REAL(n)+1.0 ) )
147  END IF
148 
149 END FUNCTION facul
150 
151 
152 
153 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
154 !
155 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
156 
157 COMPLEX FUNCTION anroot(dens,omeg,m)
158 
159  IMPLICIT NONE
160 
161  !-----------------------------------------------------------------------------
162  REAL, INTENT(IN) :: dens
163  REAL, INTENT(IN OUT) :: omeg
164  INTEGER, INTENT(IN) :: m
165 
166  !-----------------------------------------------------------------------------
167  REAL :: coeff(3),qq,rr,aa,bb,ohm, pi
168  REAL, SAVE :: roh2
169  COMPLEX, SAVE :: anr0,anr1,anr2
170  !-----------------------------------------------------------------------------
171  DATA roh2,anr0,anr1,anr2 / 0.0, (0.0, 0.0), (0.0, 0.0), (0.0, 0.0) /
172 
173 
174  pi = 3.14159265359
175 
176  IF (dens /= roh2) THEN
177  roh2 = dens
178  coeff(1)=(6.0*dens-(1.0-dens)*omeg) /(1.0-dens)
179  coeff(2)=(18.0*dens**2 - 6.0*dens*(1.0-dens)*omeg) /((1.0-dens)**2)
180  coeff(3)=(-12.0*dens*(1.0+2.0*dens) +12.0*dens*(1.0-dens)*omeg) /((1.0-dens)**2)
181  qq = ( coeff(1)**2 - 3.0*coeff(2) ) / 9.0
182  rr = ( 2.0*coeff(1)**3 - 9.0*coeff(1)*coeff(2) + 27.0*coeff(3) ) / 54.0
183 
184  IF ( rr**2 < qq**3 ) THEN
185 
186  ohm = acos( rr/qq**(3.0/2.0) )
187  anr0=-2.0* sqrt(qq)* cos( ohm/3.0 )-coeff(1)/3.0
188  anr1=-2.0* sqrt(qq)* cos((ohm+2.0*pi)/3.0)-coeff(1)/3.0
189  anr2=-2.0* sqrt(qq)* cos((ohm-2.0*pi)/3.0)-coeff(1)/3.0
190 
191  ELSE
192 
193  aa = - ( abs(rr) + sqrt( rr**2 -qq**3 ) )**(1.0/3.0)
194  aa = rr/ abs(rr) * aa
195  IF ( aa /= 0.0 ) THEN
196  bb = qq/aa
197  ELSE
198  bb = 0.0
199  END IF
200 
201  anr0 = cmplx( aa+bb -coeff(1)/3.0 , 0.0 )
202  anr1 = cmplx((-(aa+bb)/2.0 -coeff(1)/3.0), (sqrt(3.0)/2.0*(aa-bb)))
203  anr2 = cmplx((-(aa+bb)/2.0 -coeff(1)/3.0), -(sqrt(3.0)/2.0*(aa-bb)))
204 
205  END IF
206 
207  END IF
208 
209  IF (m == 0) THEN
210  anroot = anr0
211  ELSE IF (m == 1) THEN
212  anroot = anr1
213  ELSE IF (m == 2) THEN
214  anroot = anr2
215  ELSE
216  WRITE (*,*) ' Only 3 roots do exist. Asked for: ',m+1
217  stop
218  END IF
219 
220 
221 
222  ! write (*,*) '1',anroot,(anroot**3
223  ! & + coeff(1)*anroot**2 + coeff(2)*anroot + coeff(3))
224  ! DO 2 i = 1,3
225  ! polyn =anroot**3 +coeff(1)*anroot**2 +coeff(2)*anroot+coeff(3)
226  ! deriv =3.0*anroot**2 +2.0*coeff(1)*anroot+coeff(2)
227  ! realrt = anroot
228  ! imagrt = AIMAG(anroot)
229  ! realpt = polyn/deriv
230  ! imagpt = AIMAG(polyn/deriv)
231  ! realrt = realrt - realpt
232  ! imagrt = imagrt - imagpt
233  !c anroot = anroot - polyn/deriv
234  ! anroot = CMPLX(realrt,imagrt)
235  ! dummy = anroot
236  ! write (*,*) '2',polyn/deriv,anroot
237  ! 2 CONTINUE
238  ! IF (m.eq.2) stop
239 
240 END FUNCTION anroot
241 
242 
243 
244 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
245 !
246 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
247 
248 COMPLEX FUNCTION afun (dense,om,an1,an2,ak,aalpha)
249 
250  USE rdf_variables
251  IMPLICIT NONE
252 
253  !-----------------------------------------------------------------------------
254  INTEGER :: i,j,an1,an2,ak
255  REAL :: dense,om
256  COMPLEX :: aalpha, af
257  !-----------------------------------------------------------------------------
258 
259  af = 0.0
260 
261  DO i = 0, an2
262  DO j = 0, an2
263  IF ((i+j) <= an2 .AND. (i+j+j) >= (ak-an1)) THEN
264 
265  af = af + fac(an2)*fac(i+an1+j+j) &
266  /(fac(i)*fac(j)*fac(an2-i-j)*fac(i+j+j+an1-ak)) &
267  * (1.0+dense/2.0-(1.0-dense)/2.0*om)**REAL(i) &
268  * (1.0+2.0*dense-(1.0-dense)*om)**REAL(an2-i-j) &
269  * ((1.0-dense)**2 *om/(12.0*dense))**REAL(j) &
270  * aalpha**REAL(an1+i+j+j-ak)
271 
272  END IF
273  END DO
274  END DO
275 
276  afun = af
277 
278 END FUNCTION afun
279 
280 
281 
282 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
283 !
284 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
285 
286 COMPLEX FUNCTION bfun (dense,om,bn1,bn2,bn3,bi,bal)
287 
288  USE rdf_variables
289  IMPLICIT NONE
290 
291  !-----------------------------------------------------------------------------
292  INTEGER :: k,l,bn1,bn2,bn3,bi,bal,beta,gamma
293  REAL :: dense, om, imsum,factor
294  COMPLEX :: talpha, anroot, afun, bf,comfra
295  !-----------------------------------------------------------------------------
296 
297  IF (bal == 0) THEN
298  beta = 1
299  gamma= 2
300  ELSE IF (bal == 1) THEN
301  beta = 0
302  gamma= 2
303  ELSE IF (bal == 2) THEN
304  beta = 0
305  gamma= 1
306  ELSE
307  WRITE(*,*) ' no proper value for bal in BFUN !!'
308  stop
309  END IF
310 
311  ! write(*,*) ' BFUN ',dense,om,bn1,bn2,bn3,bi,bal
312  bf = (0.0, 0.0)
313  imsum = 0.0
314 
315  DO k = 0, (bn3-bi)
316  DO l = 0, (bn3-bi-k)
317  talpha = anroot(dense,om,bal)
318  ! bf = bf
319  ! & +(-1.0)**REAL(bn3-bi-k)*fac(bn3-1+l)*fac(bn3+bn3-1-bi-k-l)
320  ! & /( fac(k)*fac(l)*fac(bn3-bi-k-l)*fac(bn3-1)*fac(bn3-1) )
321  ! & *AFUN(dense,om,bn1,bn2,k,talpha)
322  ! & /( (talpha
323  ! & - anroot(dense,om,beta))**REAL(bn3+l)
324  ! & *(talpha
325  ! & - anroot(dense,om,gamma))**REAL(bn3+bn3-bi-k-l) )
326 
327  factor=(-1.0)**REAL(bn3-bi-k)*fac(bn3-1+l)*fac(bn3+bn3-1-bi-k-l) &
328  /( fac(k)*fac(l)*fac(bn3-bi-k-l)*fac(bn3-1)*fac(bn3-1) )
329  comfra= afun(dense,om,bn1,bn2,k,talpha) /( (talpha &
330  - anroot(dense,om,beta))**REAL(bn3+l) *(talpha &
331  - anroot(dense,om,gamma))**REAL(bn3+bn3-bi-k-l) )
332  imsum = imsum + factor*aimag(comfra)
333  bf = bf + factor * comfra
334 
335  END DO
336  END DO
337 
338  bfun = 1.0 / (1.0-dense)**REAL(bn3+bn3) * bf
339 
340 END FUNCTION bfun
341 
342 
343 
344 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
345 !
346 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
347 
348 REAL FUNCTION cfun (dense,om,cn1,cn2,cn3,rad)
349 
350  USE rdf_variables
351  IMPLICIT NONE
352 
353  !-----------------------------------------------------------------------------
354  INTEGER :: cn1,cn2,cn3, i, alph
355  REAL :: dense,om,rad,fact
356  COMPLEX :: anroot,bfun,cf,root,swing,compl
357  !-----------------------------------------------------------------------------
358 
359  cf = cmplx( 0.0, 0.0 )
360  IF (abs(rad) < 1.e-10) rad = 0.0
361 
362  DO alph = 0,2
363  DO i = 1,cn3
364  IF (rad >= 0.) THEN
365  root = anroot(dense,om,alph)
366  swing = root*rad
367  compl =bfun(dense,om,cn1,cn2,cn3,i,alph)*cexp(swing)
368  fact = rad**REAL(i-1) / fac(i-1)
369  cf = cf + compl * fact
370  END IF
371  END DO
372  END DO
373  cfun = cf
374 
375  IF ( (cn1+cn2) < (cn3+cn3+cn3) ) THEN
376  cfun = cf
377  ELSE IF ( (cn1+cn2) == (cn3+cn3+cn3) ) THEN
378  IF (rad == 0.0) THEN
379  cfun = cf + (1.0+dense/2.0)**REAL(cn2) /((1.0-dense)**REAL(cn3+cn3))
380  WRITE (*,*) 'in Schleife',cfun
381  ELSE
382  cfun = cf
383  END IF
384  ELSE
385  WRITE(*,*) ' no proper case in CFUN !!'
386  stop
387  END IF
388 
389 END FUNCTION cfun