MoDeNa  1.0
Software framework facilitating sequential multi-scale modelling
DFT-utilities.f90
1 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
2 ! Module DFT_MODULE
3 !
4 ! This module contains parameters and variables for DFT calculations.
5 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
6 
7 Module dft_module
8 
9  use parameters, only: nc
10  implicit none
11  save
12 
13  INTEGER, PARAMETER :: ndft = 4000
14  INTEGER, PARAMETER :: r_grid = 200
15  INTEGER :: kmax, den_step
16  LOGICAL :: shift, wca, mft
17  REAL :: rc, rg, dzr, dzp, tau_cut
18  REAL :: d_hs, dhs_st, z3t_st
19  REAL :: z_ges
20  REAL, DIMENSION(r_grid) :: x1a
21  REAL, DIMENSION(NDFT) :: x2a
22  REAL, DIMENSION(r_grid,NDFT) :: ya, y1a, y2a, y12a
23  REAL, DIMENSION(r_grid,NDFT,4,4) :: c_bicub
24  REAL :: fres_temp
25 
26  REAL, DIMENSION(r_grid) :: x1a_11
27  REAL, DIMENSION(NDFT) :: x2a_11
28  REAL, DIMENSION(r_grid,NDFT) :: ya_11, y1a_11, y2a_11, y12a_11
29  REAL, DIMENSION(r_grid,NDFT,4,4) :: c_bicub_11
30  REAL, DIMENSION(r_grid) :: x1a_12
31  REAL, DIMENSION(NDFT) :: x2a_12
32  REAL, DIMENSION(r_grid,NDFT) :: ya_12, y1a_12, y2a_12, y12a_12
33  REAL, DIMENSION(r_grid,NDFT,4,4) :: c_bicub_12
34  REAL, DIMENSION(r_grid) :: x1a_22
35  REAL, DIMENSION(NDFT) :: x2a_22
36  REAL, DIMENSION(r_grid,NDFT) :: ya_22, y1a_22, y2a_22, y12a_22
37  REAL, DIMENSION(r_grid,NDFT,4,4) :: c_bicub_22
38 
39 End Module dft_module
40 
41 
42 
43 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
44 ! SUBROUTINE DFT_RAD_INT
45 !
46 ! This subroutine integrates the kernel of the perturbation theory in
47 ! radial direction (more precisely in distance coordinate r^ hat). A
48 ! non-mean field approach is taken, using a radial distribution function
49 ! (currently at a Percus-Yevick level).
50 !
51 ! The first and second order contributions to the perturbation theory
52 ! are calculated - although the second order contribution is rarely used.
53 !
54 ! ToDo: comment the subroutine and remove the GOTO constructs
55 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
56 
57 SUBROUTINE dft_rad_int (i,j,ih,zz1,rhop,f_int_r,my_int_r, &
58  f_int2_r,my_int2_r,my_int3_r)
59 
60  USE dft_module
61  IMPLICIT NONE
62 
63  !-----------------------------------------------------------------------------
64  INTEGER, INTENT(IN) :: i
65  INTEGER, INTENT(IN) :: j
66  INTEGER, INTENT(IN OUT) :: ih
67  REAL, INTENT(IN) :: zz1
68  REAL, INTENT(IN) :: rhop(-ndft:ndft)
69  REAL, INTENT(OUT) :: f_int_r
70  REAL, INTENT(OUT) :: my_int_r
71  REAL, INTENT(OUT) :: f_int2_r
72  REAL, INTENT(OUT) :: my_int2_r
73  REAL, INTENT(OUT) :: my_int3_r
74 
75  !-----------------------------------------------------------------------------
76  INTEGER :: k
77 
78  REAL :: fint0, fint0_2, fint1, fint1_2
79  REAL :: myint0, myint0_2, myint0_3, myint1
80  REAL :: myint1_2, myint1_3
81  REAL :: dg_drho, dg_dr
82  REAL :: rad, xg, rdf, rho_bar, ua, rs
83  REAL :: analytic1, analytic2, tau_rs
84  LOGICAL :: shortcut
85  !-----------------------------------------------------------------------------
86 
87  shortcut = .true.
88  fint0 = rc * tau_cut ! first order term
89  fint0_2 = rc * tau_cut*tau_cut ! 2nd order
90  myint0 = rc * tau_cut ! first order term
91  myint0_2 = rc * tau_cut*tau_cut ! 2nd order
92  myint0_3 = 0.0 ! 2nd order
93 
94  f_int_r = 0.0
95  f_int2_r = 0.0
96  my_int_r = 0.0
97  my_int2_r= 0.0
98  my_int3_r= 0.0
99  !-------------this block only speeds up the integration
100  IF ( shortcut ) THEN
101  rs = max( rg, zz1 ) ! +dzr
102  IF ( rs > rc ) WRITE (*,*) 'error !!!!'
103  analytic1 = 0.4*rs**-10 - 0.4*rc**-10 - rs**-4 + rc**-4
104  analytic2 = 16.0/22.0 * (rs**-22 - rc**-22 ) &
105  -2.0*rs**-16 +2.0*rc**-16 +1.6*rs**-10 - 1.6*rc**-10
106  f_int_r = f_int_r + analytic1
107  f_int2_r = f_int2_r + analytic2
108  my_int_r = my_int_r + analytic1
109  my_int2_r= my_int2_r+ analytic2
110  IF ( rs == zz1 ) GO TO 10
111  tau_rs = 4.0 * ( rs**-12 - rs**-6 )
112  fint0 = rs * tau_rs
113  fint0_2 = rs * tau_rs*tau_rs
114  myint0 = rs * tau_rs
115  myint0_2= rs * tau_rs*tau_rs
116  rad = rs ! the simple integration scheme: set to rc
117  k = 0 + nint( (rc-rs)/dzr ) ! in simple scheme: set to 0
118  ELSE
119  rad = rc
120  k = 0
121  END IF
122  !-------------
123  !rad = rc
124  !k = 0
125 
126 4 CONTINUE
127  rad = rad - dzr
128  k = k + 1
129  ua = 4.0 * ( rad**-12 -rad**-6 )
130  xg = rad / d_hs
131  rho_bar = ( rhop(i) + rhop(j) )/2.0
132  rdf = 1.0
133  dg_drho = 0.0
134  IF ( rad <= rg ) THEN
135  CALL bi_cub_spline (rho_bar,xg,ya,x1a,x2a,y1a,y2a,y12a, &
136  c_bicub,rdf,dg_drho,dg_dr,den_step,ih,k)
137  ! write (*,*) 'vor ',rdf,dg_drho
138  ! CALL BI_CUB_SPLINE2 (rho_bar,xg,ya,x1a,x2a,y1a,
139  ! & rdf,dg_drho,den_step,ih,k)
140  ! write (*,*) 'nach',rdf,dg_drho
141  ! read (*,*)
142  END IF
143 
144  fint1 = rdf * rad * ua
145  fint1_2 = rdf * rad * ua * ua
146  myint1 = rad * (rdf + 0.5*rhop(i)*dg_drho) * ua
147  myint1_2 = rdf * rad * ua * ua
148  myint1_3 = dg_drho * rad * ua * ua
149  ! intf(k) = fint1
150  f_int_r = f_int_r + dzr * (fint1 + fint0)/2.0
151  f_int2_r = f_int2_r + dzr * (fint1_2 + fint0_2)/2.0
152  my_int_r = my_int_r + dzr * (myint1 + myint0)/2.0
153  my_int2_r= my_int2_r+ dzr * (myint1_2 + myint0_2)/2.0
154  my_int3_r= my_int3_r+ dzr * (myint1_3 + myint0_3)/2.0
155 
156  fint0 = fint1
157  fint0_2 = fint1_2
158  myint0 = myint1
159  myint0_2= myint1_2
160  myint0_3= myint1_3
161  IF ( zz1 >= 1.0 .AND. rad-dzr+1.e-8 >= zz1 ) GO TO 4 ! integration down to ABS(zz1)
162  IF ( zz1 < 1.0 .AND. rad-dzr+1.e-8 >= 1.0 ) GO TO 4 ! integration down to ABS(zz1) but for r^hat<1, g(r)=0 (so stop at r^hat=1)
163 
164  ! IF (k.GT.30) THEN
165  ! stepno = k
166  ! CALL SPLINE_PARA (dzr,intf,utri,stepno)
167  ! CALL SPLINE_INT (f_int_r,dzr,intf,utri,stepno)
168  ! ENDIF
169 
170 10 CONTINUE
171 
172  analytic1 = 4.0/10.0*rc**-10 - rc**-4
173  analytic2 = 16.0/22.0*rc**-22 - 2.0*rc**-16 + 1.6*rc**-10
174  f_int_r = f_int_r + analytic1
175  f_int2_r = f_int2_r + analytic2
176  my_int_r = my_int_r + analytic1
177  my_int2_r= my_int2_r+ analytic2
178 
179 END SUBROUTINE dft_rad_int
180 
181 
182 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
183 !
184 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
185 
186 SUBROUTINE f_dft ( I1_dft, I2_dft )
187 
188  USE eos_variables, ONLY: nc, pi, rho, parame
189  USE dft_module
190  IMPLICIT NONE
191 
192  !-----------------------------------------------------------------------------
193  REAL, INTENT(OUT) :: i1_dft
194  REAL, INTENT(OUT) :: i2_dft
195 
196  !-----------------------------------------------------------------------------
197  INTEGER :: k,ih
198  ! REAL :: z3
199  REAL :: ua, ua_c, ua_2, ua_c_2, rm
200  REAL :: int10, int11, int20, int21
201  REAL :: dg_drho
202  REAL :: rad, xg, rdf, rho_st, msegm
203  REAL :: sig_ij
204  REAL :: dg_dr, dzr_org !,rdf_d
205  ! REAL :: intgrid(0:NDFT),intgri2(0:NDFT)
206  !-----------------------------------------------------------------------------
207 
208  !-----------------------------------------------------------------------------
209  ! constants
210  !-----------------------------------------------------------------------------
211  msegm = parame(1,1)
212  rho_st = rho * parame(1,2)**3
213 
214  ua_c = 4.0 * ( rc**-12 - rc**-6 )
215  ua_c_2 = ua_c * ua_c
216  rm = 2.0**(1.0/6.0)
217 
218  int10 = rc*rc* ua_c
219  int20 = rc*rc* ua_c_2
220  ! intgrid(0)= int10
221  ! intgri2(0)= int20
222 
223 
224  sig_ij = parame(1,2)
225 
226 
227  i1_dft = 0.0
228  i2_dft = 0.0
229  rad = rc
230  !dzr = dzp / 2.0 ! this line is obsolete. dzr is defined in DFT-nMF2 (dimensionless)
231  dzr_org= dzr
232  k = 0
233  ih = 85
234 
235  DO WHILE ( rad-dzr+1.e-9 >= 1.0 )
236 
237  rad = rad - dzr
238  ! IF (rad <= 8.0) dzr = dzp
239  ! IF (rad <= rg) dzr = dzp/2.0
240  k = k + 1
241  xg = rad / dhs_st
242  ua = 4.0 * ( rad**-12 - rad**-6 )
243  ua_2 = ua * ua
244  rdf = 1.0
245  dg_drho = 0.0
246  IF ( rad <= rg ) THEN
247  CALL bi_cub_spline (rho_st,xg,ya,x1a,x2a,y1a,y2a,y12a, &
248  c_bicub,rdf,dg_drho,dg_dr,den_step,ih,k)
249  END IF
250 
251  int11 = rdf*rad*rad* ua
252  int21 = rdf*rad*rad* ua_2
253  i1_dft= i1_dft + dzr*(int11+int10)/2.0
254  i2_dft= i2_dft + dzr*(int21+int20)/2.0
255  int10 = int11
256  int20 = int21
257 
258  END DO
259 
260  dzr = dzr_org
261 
262  ! stepno = k
263  ! CALL SPLINE_PARA (dzr,intgrid,utri,stepno)
264  ! CALL SPLINE_INT (I1,dzr,intgrid,utri,stepno)
265 
266  ! caution: 1st order integral is in F_EOS.f defined with negative sign
267  i1_dft= - i1_dft - ( 4.0/9.0 * rc**-9 - 4.0/3.0 * rc**-3 )
268 
269  ! CALL SPLINE_PARA (dzr,intgri2,utri,stepno)
270  ! CALL SPLINE_INT (I2,dzr,intgri2,utri,stepno)
271 
272  i2_dft = i2_dft + 16.0/21.0 * rc**-21 - 32.0/15.0 * rc**-15 + 16.0/9.0 * rc**-9
273 
274 
275 END SUBROUTINE f_dft
276 
277 
278 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
279 !
280 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
281 
282 SUBROUTINE rdf_matrix (msegm)
283 
284  USE parameters, ONLY: pi
285  USE dft_module
286  IMPLICIT NONE
287 
288  !-----------------------------------------------------------------------------
289  REAL, INTENT(IN) :: msegm
290 
291  !-----------------------------------------------------------------------------
292  INTEGER :: i, k = 0
293  REAL :: rdf, rad, xg, rho_rdf, z3
294  !-----------------------------------------------------------------------------
295 
296 
297  DO i = 1, den_step
298 
299  ! rho_rdf= rhob(1) + (rhob(2)-rhob(1))*REAL(i)/den_step
300  rho_rdf = 1.e-5 + (1.0) * REAL(i-1) / REAL(den_step-1) ! segment density, rho_s*sigma**3
301  rho_rdf = rho_rdf / msegm ! molar density, rho_m*sigma**3
302  rad = rc
303  k = 0
304 
305  do while ( rad-dzr+1.e-9 >= 0.95 )
306 
307  rad = rad - dzr
308  k = k + 1
309  xg = rad / d_hs
310  z3 = rho_rdf * msegm * pi/6.0 * d_hs**3
311  ya(i,k) = 1.0
312  IF ( xg <= rg .AND. z3 > 0.0 ) CALL rdf_int ( z3, msegm, xg, rdf )
313  IF ( xg <= rg ) ya(i,k) = rdf
314  ! ya(i,k) = y(x1a(i), x2a(k)) with x1a: density-vector, x2a: r-vector
315  x1a(i) = rho_rdf
316  x2a(k) = xg
317 
318  end do
319 
320  END DO
321 
322  if ( xg > 1.0 ) stop 'rdf_matrix: 0.95*sigma is too high for lower bound'
323 
324  WRITE (*,*) ' done with calculating g(r)',d_hs
325 
326  kmax = k
327  CALL bicub_derivative ( ya, x1a, x2a, y1a, y2a, y12a, den_step, kmax )
328  CALL bicub_c ( ya, x1a, x2a, y1a, y2a, y12a, c_bicub, den_step, kmax )
329 
330 END SUBROUTINE rdf_matrix
331 
332 
333 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
334 !
335 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
336 
337 SUBROUTINE bi_cub_spline2 ( rho_rdf, xg, ya, x1a, x2a, y1a, &
338  rdf, dg_drho, i_max, ih, k )
339 
340  USE dft_module, ONLY: ndft, r_grid
341  IMPLICIT NONE
342 
343  !-----------------------------------------------------------------------------
344  REAL, INTENT(IN) :: rho_rdf
345  REAL, INTENT(IN OUT) :: xg
346  REAL, INTENT(IN) :: ya(r_grid,ndft)
347  REAL, INTENT(IN) :: x1a(r_grid)
348  REAL, INTENT(IN OUT) :: x2a(ndft)
349  REAL, INTENT(IN) :: y1a(r_grid,ndft)
350  REAL, INTENT(OUT) :: rdf
351  REAL, INTENT(OUT) :: dg_drho
352  INTEGER, INTENT(IN OUT) :: i_max
353  !INTEGER, INTENT(IN OUT) :: k_max
354  INTEGER, INTENT(OUT) :: ih
355  INTEGER, INTENT(IN OUT) :: k
356 
357  !-----------------------------------------------------------------------------
358  REAL :: as,bs,cs,ds,rho_t,del_rho
359  !REAL :: y2a(r_grid,NDFT),y12a(r_grid,NDFT), c_bicub(r_grid,NDFT,4,4),c(4,4)
360  !-----------------------------------------------------------------------------
361 
362  IF ( rho_rdf < x1a(1) ) THEN
363  dg_drho = 0.0
364  rdf = 1.0
365  RETURN
366  END IF
367  IF ( x1a(ih) <= rho_rdf.AND.rho_rdf < x1a(ih+1) ) GO TO 10
368  IF ( ih > 2 ) THEN
369  IF ( x1a(ih-1) <= rho_rdf.AND.rho_rdf < x1a(ih) ) THEN
370  ih = ih - 1
371  GO TO 10
372  END IF
373  END IF
374  ! write (*,*) 'in ',ih
375  CALL hunt(x1a,i_max,rho_rdf,ih)
376  ! write (*,*) 'out',ih
377 10 CONTINUE
378  IF ( (x2a(k) > (xg + 1.e-10)) .OR. (x2a(k) < (xg - 1.e-10)) ) &
379  WRITE (*,*) 'error in bi-cubic-spline',x2a(k),xg
380 
381 
382  del_rho = x1a(ih+1) - x1a(ih)
383 
384  as = ya(ih,k)
385  bs = y1a(ih,k)
386  cs = 3.0*(ya(ih+1,k)-as)/del_rho**2 - (2.0*bs+y1a(ih+1,k))/del_rho
387  ds = -2.0*(ya(ih+1,k)-as)/del_rho**3 + (bs+y1a(ih+1,k))/del_rho**2
388 
389  rho_t = rho_rdf - x1a(ih)
390  rdf = as + bs*rho_t + cs*rho_t*rho_t + ds*rho_t**3
391  dg_drho = bs + 2.0*cs*rho_t + 3.0*ds*rho_t*rho_t
392 
393 END SUBROUTINE bi_cub_spline2
394 
395 
396 
397 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
398 !
399 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
400 
401 SUBROUTINE bi_cub_spline ( rho_rdf, xg, ya, x1a, x2a, y1a, y2a, y12a, &
402  c_bicub, rdf, dg_drho, dg_dr, i_max, ih, k )
403 
404  USE dft_module, ONLY: ndft, r_grid
405  IMPLICIT NONE
406 
407  !-----------------------------------------------------------------------------
408  REAL, INTENT(IN OUT) :: rho_rdf
409  REAL, INTENT(IN OUT) :: xg
410  REAL, INTENT(IN) :: ya(r_grid,ndft)
411  REAL, INTENT(IN) :: x1a(r_grid)
412  REAL, INTENT(IN) :: x2a(ndft)
413  REAL, INTENT(IN) :: y1a(r_grid,ndft)
414  REAL, INTENT(IN) :: y2a(r_grid,ndft)
415  REAL, INTENT(IN) :: y12a(r_grid,ndft)
416  REAL, INTENT(IN) :: c_bicub(r_grid,ndft,4,4)
417  REAL, INTENT(OUT) :: rdf
418  REAL, INTENT(OUT) :: dg_drho
419  REAL, INTENT(OUT) :: dg_dr
420  INTEGER, INTENT(IN OUT) :: i_max
421  !INTEGER, INTENT(IN OUT) :: k_max
422  INTEGER, INTENT(OUT) :: ih
423  INTEGER, INTENT(IN) :: k
424 
425  !-----------------------------------------------------------------------------
426  INTEGER :: m, n
427 
428  REAL :: y(4),y1(4),y2(4),y12(4),x1l,x1u,x2l,x2u
429  REAL :: c(4,4)
430  !-----------------------------------------------------------------------------
431 
432  IF ( rho_rdf < x1a(1) ) THEN
433  dg_drho = 0.0
434  dg_dr = 0.0
435  rdf = 1.0
436  RETURN
437  END IF
438  IF ( x1a(ih) <= rho_rdf .AND. rho_rdf < x1a(ih+1) ) GO TO 10
439  IF ( ih > 2 ) THEN
440  IF ( x1a(ih-1) <= rho_rdf .AND. rho_rdf < x1a(ih) ) THEN
441  ih = ih - 1
442  GO TO 10
443  END IF
444  END IF
445  ! write (*,*) 'in ',ih
446  CALL hunt ( x1a, i_max, rho_rdf, ih )
447  ! write (*,*) 'out',ih
448 10 CONTINUE
449  IF ( x2a(k) /= xg ) THEN
450  ! write (*,*) 'error bi-cubic-spline',k,x2a(k),xg
451  ! DO k=1,NDFT
452  ! write (*,*) k,x2a(k)
453  ! ENDDO
454  ! stop
455  END IF
456 
457 
458 
459  y(1) = ya(ih,k)
460  y(2) = ya(ih+1,k)
461  y(3) = ya(ih+1,k+1)
462  y(4) = ya(ih,k+1)
463 
464  y1(1) = y1a(ih,k)
465  y1(2) = y1a(ih+1,k)
466  y1(3) = y1a(ih+1,k+1)
467  y1(4) = y1a(ih,k+1)
468 
469  y2(1) = y2a(ih,k)
470  y2(2) = y2a(ih+1,k)
471  y2(3) = y2a(ih+1,k+1)
472  y2(4) = y2a(ih,k+1)
473 
474  y12(1) = y12a(ih,k)
475  y12(2) = y12a(ih+1,k)
476  y12(3) = y12a(ih+1,k+1)
477  y12(4) = y12a(ih,k+1)
478 
479  x1l = x1a(ih)
480  x1u = x1a(ih+1)
481  x2l = x2a(k)
482  x2u = x2a(k+1)
483 
484  DO m = 1, 4
485  DO n = 1, 4
486  c(m,n) = c_bicub( ih, k, m, n )
487  END DO
488  END DO
489  CALL bcuint ( x1l, x1u, x2l, x2u, rho_rdf, xg, c, rdf, dg_drho, dg_dr )
490 
491 END SUBROUTINE bi_cub_spline
492 
493 
494 
495 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
496 !
497 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
498 
499 SUBROUTINE bicub_derivative ( ya, x1a, x2a, y1a, y2a, y12a, i_max, k_max )
500 
501  USE dft_module, ONLY: ndft, r_grid
502  IMPLICIT NONE
503 
504  !-----------------------------------------------------------------------------
505  REAL, INTENT(IN) :: ya(r_grid,ndft)
506  REAL, INTENT(IN) :: x1a(r_grid)
507  REAL, INTENT(IN) :: x2a(ndft)
508  REAL, INTENT(OUT) :: y1a(r_grid,ndft)
509  REAL, INTENT(OUT) :: y2a(r_grid,ndft)
510  REAL, INTENT(OUT) :: y12a(r_grid,ndft)
511  INTEGER, INTENT(IN) :: i_max
512  INTEGER, INTENT(IN) :: k_max
513 
514  !-----------------------------------------------------------------------------
515  INTEGER :: i, k
516  !-----------------------------------------------------------------------------
517 
518 
519  DO i = 2, i_max-1
520  DO k = 2, k_max-1
521  y1a(i,k) = (ya(i+1,k)-ya(i-1,k)) / (x1a(i+1)-x1a(i-1))
522  y2a(i,k) = (ya(i,k+1)-ya(i,k-1)) / (x2a(k+1)-x2a(k-1))
523  y12a(i,k)= (ya(i+1,k+1)-ya(i+1,k-1)-ya(i-1,k+1)+ya(i-1,k-1)) &
524  /((x1a(i+1)-x1a(i-1))*(x2a(k+1)-x2a(k-1)))
525  END DO
526  END DO
527 
528  i = 1
529  DO k = 1, k_max
530  y1a(i,k) = 0.0
531  y2a(i,k) = 0.0
532  y12a(i,k)= 0.0
533  END DO
534 
535  k = 1
536  DO i = 1, i_max
537  y1a(i,k) = 0.0
538  y2a(i,k) = 0.0
539  y12a(i,k)= 0.0
540  END DO
541 
542 
543  i = i_max
544  DO k = 2, k_max-1
545  y1a(i,k) = (ya(i,k)-ya(i-1,k)) / (x1a(i)-x1a(i-1))
546  y2a(i,k) = (ya(i,k+1)-ya(i,k-1)) / (x2a(k+1)-x2a(k-1))
547  y12a(i,k)= (ya(i,k+1)-ya(i,k-1)-ya(i-1,k+1)+ya(i-1,k-1)) &
548  /((x1a(i)-x1a(i-1))*(x2a(k+1)-x2a(k-1)))
549  END DO
550 
551 
552  k = k_max
553  DO i = 2, i_max-1
554  y1a(i,k) = (ya(i+1,k)-ya(i-1,k)) / (x1a(i+1)-x1a(i-1))
555  y2a(i,k) = (ya(i,k)-ya(i,k-1)) / (x2a(k)-x2a(k-1))
556  y12a(i,k)= (ya(i+1,k)-ya(i+1,k-1)-ya(i-1,k)+ya(i-1,k-1)) &
557  /((x1a(i+1)-x1a(i-1))*(x2a(k)-x2a(k-1)))
558  END DO
559 
560  k = k_max
561  i = i_max
562  y1a(i,k) = (ya(i,k)-ya(i-1,k)) / (x1a(i)-x1a(i-1))
563  y2a(i,k) = (ya(i,k)-ya(i,k-1)) / (x2a(k)-x2a(k-1))
564  y12a(i,k)= (ya(i,k)-ya(i,k-1)-ya(i-1,k)+ya(i-1,k-1)) &
565  /((x1a(i)-x1a(i-1))*(x2a(k)-x2a(k-1)))
566 
567 END SUBROUTINE bicub_derivative
568 
569 
570 
571 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
572 !
573 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
574 
575 SUBROUTINE bicub_c ( ya, x1a, x2a, y1a, y2a, y12a, c_bicub, i_max, k_max )
576 
577  USE dft_module, ONLY: ndft, r_grid
578  IMPLICIT NONE
579 
580  !-----------------------------------------------------------------------------
581  REAL, INTENT(IN) :: ya(r_grid,ndft)
582  REAL, INTENT(IN) :: x1a(r_grid)
583  REAL, INTENT(IN) :: x2a(ndft)
584  REAL, INTENT(IN) :: y1a(r_grid,ndft)
585  REAL, INTENT(IN) :: y2a(r_grid,ndft)
586  REAL, INTENT(IN) :: y12a(r_grid,ndft)
587  REAL, INTENT(OUT) :: c_bicub(r_grid,ndft,4,4)
588  INTEGER, INTENT(IN) :: i_max
589  INTEGER, INTENT(IN) :: k_max
590 
591  !-----------------------------------------------------------------------------
592  INTEGER :: i, k, m, n
593  REAL :: y(4),y1(4),y2(4),y12(4),x1l,x1u,x2l,x2u
594  REAL :: c(4,4)
595  !-----------------------------------------------------------------------------
596 
597  DO i = 1, i_max-1
598  DO k = 1, k_max-1
599  y(1)=ya(i,k)
600  y(2)=ya(i+1,k)
601  y(3)=ya(i+1,k+1)
602  y(4)=ya(i,k+1)
603 
604  y1(1)=y1a(i,k)
605  y1(2)=y1a(i+1,k)
606  y1(3)=y1a(i+1,k+1)
607  y1(4)=y1a(i,k+1)
608 
609  y2(1)=y2a(i,k)
610  y2(2)=y2a(i+1,k)
611  y2(3)=y2a(i+1,k+1)
612  y2(4)=y2a(i,k+1)
613 
614  y12(1)=y12a(i,k)
615  y12(2)=y12a(i+1,k)
616  y12(3)=y12a(i+1,k+1)
617  y12(4)=y12a(i,k+1)
618 
619  x1l=x1a(i)
620  x1u=x1a(i+1)
621  x2l=x2a(k)
622  x2u=x2a(k+1)
623 
624  CALL bcucof(y,y1,y2,y12,x1u-x1l,x2u-x2l,c)
625  DO m=1,4
626  DO n=1,4
627  c_bicub(i,k,m,n)=c(m,n)
628  END DO
629  END DO
630 
631  END DO
632  END DO
633 
634 END SUBROUTINE bicub_c
635 
636 
637 
638 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
639 !
640 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
641 
642 !SUBROUTINE bcuint ( y, y1, y2, y12, x1l, x1u, x2l, x2u, x1, x2, c, ansy, ansy1, ansy2 )
643 SUBROUTINE bcuint ( x1l, x1u, x2l, x2u, x1, x2, c, ansy, ansy1, ansy2 )
644 
645  USE utilities
646  IMPLICIT NONE
647 
648  !-----------------------------------------------------------------------------
649  !REAL, INTENT(IN OUT) :: y(4)
650  !REAL, INTENT(IN OUT) :: y1(4)
651  !REAL, INTENT(IN OUT) :: y2(4)
652  !REAL, INTENT(IN OUT) :: y12(4)
653  REAL, INTENT(IN OUT) :: x1l
654  REAL, INTENT(IN OUT) :: x1u
655  REAL, INTENT(IN OUT) :: x2l
656  REAL, INTENT(IN OUT) :: x2u
657  REAL, INTENT(IN OUT) :: x1
658  REAL, INTENT(IN OUT) :: x2
659  REAL, INTENT(IN) :: c(4,4)
660  REAL, INTENT(OUT) :: ansy
661  REAL, INTENT(OUT) :: ansy1
662  REAL, INTENT(OUT) :: ansy2
663 
664  !-----------------------------------------------------------------------------
665  !U USES bcucof
666  INTEGER :: i
667  REAL :: t, u
668  !-----------------------------------------------------------------------------
669 
670  ! call bcucof ( y, y1, y2, y12, x1u-x1l, x2u-x2l, c )
671 
672  IF ( x1u == x1l .OR. x2u == x2l ) call paus ('bad input in bcuint')
673  t = (x1-x1l) / (x1u-x1l)
674  u = (x2-x2l) / (x2u-x2l)
675  ansy = 0.0
676  ansy2 = 0.0
677  ansy1 = 0.0
678  DO i = 4, 1, -1
679  ansy = t *ansy + ( (c(i,4)*u + c(i,3))*u+c(i,2) )*u + c(i,1)
680  ansy2 = t *ansy2 + ( 3.*c(i,4)*u+2.*c(i,3) )*u + c(i,2)
681  ansy1 = u *ansy1 + ( 3.*c(4,i)*t+2.*c(3,i) )*t + c(2,i)
682  END DO
683  ansy1 = ansy1 / (x1u-x1l)
684  ansy2 = ansy2 / (x2u-x2l)
685 
686 END SUBROUTINE bcuint
687 
688 
689 
690 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
691 !
692 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
693 
694 SUBROUTINE bcucof(y,y1,y2,y12,d1,d2,c)
695 
696  IMPLICIT NONE
697 
698  !-----------------------------------------------------------------------------
699  REAL, INTENT(IN) :: y(4)
700  REAL, INTENT(IN) :: y1(4)
701  REAL, INTENT(IN) :: y2(4)
702  REAL, INTENT(IN) :: y12(4)
703  REAL, INTENT(IN) :: d1
704  REAL, INTENT(IN) :: d2
705  REAL, INTENT(OUT) :: c(4,4)
706 
707  !-----------------------------------------------------------------------------
708  INTEGER :: i,j,k,l
709  REAL :: d1d2,xx,cl(16),wt(16,16),x(16)
710  SAVE wt
711  DATA wt/1,0,-3,2,4*0,-3,0,9,-6,2,0,-6,4,8*0,3,0,-9,6,-2,0,6,-4,10* &
712  0,9,-6,2*0,-6,4,2*0,3,-2,6*0,-9,6,2*0,6,-4,4*0,1,0,-3,2,-2,0,6,-4, &
713  1,0,-3,2,8*0,-1,0,3,-2,1,0,-3,2,10*0,-3,2,2*0,3,-2,6*0,3,-2,2*0, &
714  -6,4,2*0,3,-2,0,1,-2,1,5*0,-3,6,-3,0,2,-4,2,9*0,3,-6,3,0,-2,4,-2, &
715  10*0,-3,3,2*0,2,-2,2*0,-1,1,6*0,3,-3,2*0,-2,2,5*0,1,-2,1,0,-2,4, &
716  -2,0,1,-2,1,9*0,-1,2,-1,0,1,-2,1,10*0,1,-1,2*0,-1,1,6*0,-1,1,2*0, &
717  2,-2,2*0,-1,1/
718  !-----------------------------------------------------------------------------
719 
720  d1d2 = d1 * d2
721  DO i = 1, 4
722  x(i) = y(i)
723  x(i+4) = y1(i)*d1
724  x(i+8) = y2(i)*d2
725  x(i+12) = y12(i)*d1d2
726  END DO
727  DO i = 1, 16
728  xx = 0.0
729  DO k = 1, 16
730  xx = xx + wt(i,k) * x(k)
731  END DO
732  cl(i) = xx
733  END DO
734  l = 0
735  DO i = 1, 4
736  DO j = 1, 4
737  l = l + 1
738  c(i,j) = cl(l)
739  END DO
740  END DO
741 
742 END SUBROUTINE bcucof
743 
744 
745 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
746 ! SUBROUTINE hunt
747 !
748 ! Given an array xx(1:n), and given a value x, returns a value jlo
749 ! such that x is between xx(jlo) and xx(jlo+1). xx(1:n) must be
750 ! monotonic, either increasing or decreasing. jlo=0 or jlo=n is
751 ! returned to indicate that x is out of range. jlo on input is taken
752 ! as the initial guess for jlo on output.
753 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
754 
755 SUBROUTINE hunt ( xx, n, x, jlo )
756 
757  IMPLICIT NONE
758 
759  !-----------------------------------------------------------------------------
760  INTEGER, INTENT(IN) :: n
761  INTEGER, INTENT(OUT) :: jlo
762  REAL, INTENT(IN) :: xx(n)
763  REAL :: x
764 
765  !-----------------------------------------------------------------------------
766  INTEGER :: inc,jhi,jm
767  LOGICAL :: ascnd
768  !-----------------------------------------------------------------------------
769 
770  ascnd = xx(n) >= xx(1)
771  IF( jlo <= 0 .OR. jlo > n ) THEN
772  jlo = 0
773  jhi = n + 1
774  GO TO 3
775  END IF
776  inc = 1
777  IF( x >= xx(jlo) .EQV. ascnd ) THEN
778 1 jhi = jlo + inc
779  IF ( jhi > n ) THEN
780  jhi = n + 1
781  ELSE IF ( x >= xx(jhi) .EQV. ascnd ) THEN
782  jlo = jhi
783  inc = inc + inc
784  GO TO 1
785  END IF
786  ELSE
787  jhi = jlo
788 2 jlo = jhi - inc
789  IF ( jlo < 1 ) THEN
790  jlo = 0
791  ELSE IF ( x < xx(jlo) .EQV. ascnd ) THEN
792  jhi = jlo
793  inc = inc + inc
794  GO TO 2
795  END IF
796  END IF
797 3 IF (jhi-jlo == 1 ) THEN
798  IF ( x == xx(n)) jlo = n - 1
799  IF ( x == xx(1) ) jlo = 1
800  RETURN
801  END IF
802  jm = ( jhi + jlo ) / 2
803  IF ( x >= xx(jm) .EQV. ascnd ) THEN
804  jlo = jm
805  ELSE
806  jhi = jm
807  END IF
808  GO TO 3
809 END SUBROUTINE hunt
810 
811 
812 
813 
814 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
815 !
816 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
817 
818 SUBROUTINE dens_inv_coeff (rhob)
819 
820  USE basic_variables
821  USE dft_module
822  USE utilities
823  IMPLICIT NONE
824 
825  !-----------------------------------------------------------------------------
826  REAL, INTENT(IN) :: rhob(2)
827 
828  !-----------------------------------------------------------------------------
829  INTEGER :: i, no_step_l, no_step_h
830  REAL :: my_eta_div,den_st,den_min,den_max, den_div
831  REAL :: lnphi(np,nc)
832  !-----------------------------------------------------------------------------
833  REAL :: rho_array1(200),rho_array2(200), &
834  my_rho_array1(200),my_rho_array2(200)
835  COMMON /deninv/ rho_array1,rho_array2,my_rho_array1, &
836  my_rho_array2,my_eta_div,den_min,den_max, no_step_l,no_step_h
837  !-----------------------------------------------------------------------------
838 
839  den_min = rhob(2) * z3t_st * 0.8
840  den_max = rhob(1) * z3t_st * 1.1
841  den_div = 0.04
842 
843  no_step_l = 80
844  no_step_h = 200
845 
846  my_eta_div = 0.0
847  IF (den_min < den_div) THEN
848  DO i=1,no_step_l
849  dense(1) = den_min + (den_div-den_min) * REAL(i-1)/REAL(no_step_l-1)
850  densta(1)= dense(1)
851  rho_array1(i) = dense(1)
852  CALL fugacity (lnphi)
853  my_rho_array1(i) = exp( lnphi(1,1) ) *dense(1)/z3t_st / parame(1,2)**3 ! myrho = my_res(T,rho) + ln(rho)
854  END DO
855  my_eta_div = my_rho_array1(no_step_l)
856  END IF
857 
858  den_st = max( den_min, den_div )
859  DO i = 1, no_step_h
860  dense(1) = den_st + (den_max-den_st) * REAL(i-1)/REAL(no_step_h-1)
861  densta(1)= dense(1)
862  rho_array2(i) = dense(1)
863  CALL fugacity (lnphi)
864  my_rho_array2(i) = lnphi(1,1) + log( dense(1)/z3t_st /parame(1,2)**3 ) ! myrho = my_res(T,rho) + ln(rho)
865  IF ( i >= 2 ) THEN
866  IF ( (my_rho_array2(i)- my_rho_array2(i-1)) < 0.0 ) THEN
867  call paus ('DENS_INV_COEFF: derivative negative')
868  END IF
869  END IF
870  END DO
871 
872 END SUBROUTINE dens_inv_coeff
873 
874 
875 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
876 !
877 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
878 
879 REAL FUNCTION dens_invers ( rhob, my_rho )
880 
881  USE basic_variables
882  USE dft_module
883  IMPLICIT NONE
884 
885  !-----------------------------------------------------------------------------
886  REAL, INTENT(IN) :: rhob(2)
887  REAL, INTENT(IN OUT) :: my_rho
888 
889  !-----------------------------------------------------------------------------
890  INTEGER :: i, no_step_l, no_step_h
891  REAL :: my_eta_div, den_min, den_max, eta
892  !-----------------------------------------------------------------------------
893  REAL :: rho_array1(200),rho_array2(200), &
894  my_rho_array1(200),my_rho_array2(200)
895  COMMON /deninv/ rho_array1,rho_array2,my_rho_array1, &
896  my_rho_array2,my_eta_div,den_min,den_max, no_step_l,no_step_h
897  !-----------------------------------------------------------------------------
898 
899  ! den_min = rhob(2)*z3t_st * 0.8
900  ! den_max = rhob(1)*z3t_st * 1.1
901  ! den_div = 0.04
902 
903  ! no_step_l = 40
904  ! no_step_h = 100
905 
906  IF ( exp(my_rho) < my_eta_div ) THEN
907  DO i = 2, no_step_l
908  IF ( my_rho_array1(i) > exp(my_rho) ) THEN
909  eta = rho_array1(i-1) + (rho_array1(i)-rho_array1(i-1)) &
910  / (my_rho_array1(i)-my_rho_array1(i-1)) &
911  * (exp(my_rho)-my_rho_array1(i-1))
912  dens_invers = eta/z3t_st
913  GO TO 1
914  END IF
915  END DO
916  WRITE (*,*) 'error 1',exp(my_rho),my_eta_div
917  ! call paus (' ')
918  END IF
919 1 CONTINUE
920 
921  IF ( exp(my_rho) >= my_eta_div ) THEN
922  DO i = 2, no_step_h
923  IF ( my_rho_array2(i) > my_rho ) THEN
924  eta = rho_array2(i-1)+(rho_array2(i)-rho_array2(i-1)) &
925  / (my_rho_array2(i)-my_rho_array2(i-1)) *(my_rho-my_rho_array2(i-1))
926  dens_invers = eta/z3t_st
927  GO TO 2
928  END IF
929  END DO
930  WRITE (*,*) 'error 2',my_rho,my_eta_div
931  ! call paus (' ')
932  END IF
933 2 CONTINUE
934 
935  IF ( dens_invers < (rhob(2)*0.8) ) THEN
936  WRITE (*,*) 'lower bound',dens_invers,my_rho
937  dens_invers = rhob(2)
938  END IF
939  IF ( dens_invers > (rhob(1)*1.1) ) THEN
940  WRITE (*,*) 'upper bound',dens_invers,my_rho
941  dens_invers = rhob(1)
942  END IF
943 
944 END FUNCTION dens_invers
parameters
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains constant...
Definition: modules.f90:6
dft_module
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW Module DFT_MODULE This module...
Definition: modules.f90:272
eos_variables
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:120
basic_variables
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:29