MoDeNa  1.0
Software framework facilitating sequential multi-scale modelling
DFT-nMF-mixtures.f90
1 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
2 ! SUBROUTINE DFT_nMFT2
3 !
4 ! This subroutine performs a DFT calculation of the PCP-SAFT model using
5 ! the fundamental measure theory (FMT) for the hard-sphere fluid, the
6 ! TPT1-theory for chain formation and a first (or second) order perturba-
7 ! tion theory for the dispersive interactions.
8 !
9 ! Background to the code:
10 ! The functional derivative of d(F)/d(rho) is done in two parts, as
11 ! 'myrho(i) + mu_att(i)' , where 'mu_att' comprises the perturbation
12 ! theory (dispersive attraction) as well as the non-local part of the hard-
13 ! sphere term and chain term. For the hard-sphere term mu_att contains
14 ! contributions of the form mu_FMT = d(F^hs)/d(rho*) - (dF^hs/drho*)_LDA,
15 ! where the index 'hs' is for hard-sphere (analogous for the chain term)
16 ! and 'LDA' is for 'local density approximation', i.e. the values for the
17 ! local density.
18 !
19 ! The term 'myrho' contains the hard-sphere term and the chain term in
20 ! the LDA. Because the PC-SAFT dispersion term is not identical to the
21 ! Barker-Henderson (BH) perturtation theory (first or second order) for
22 ! LJ fluids, the term 'myrho' also contains the difference of the PC-SAFT
23 ! dispersion contribution and the BH perturbation theory used here the
24 ! LDA. This difference is obtained by using the flag {subtract1 = '1PT'}
25 ! the effect of the flag is that the BH-perturbation theory is subtracted
26 ! from the PC-SAFT model for calculating 'myrho'.
27 !
28 ! Caution: the indication for using the second order term in the pertur-
29 ! bation expansion is currently soley given by the flag {subtract1 = '2PT'}
30 ! The second order case needs to be re-tested before usage !
31 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
32 
33 SUBROUTINE dft_nmft_mix
34 
35  USE parameters, ONLY: pi, rgas, kbol
36  USE basic_variables
37  USE eos_variables, ONLY: fres, eta, dhs, mseg, uij, sig_ij, rho, x
38  USE starting_values
39  USE dft_module
40  use utilities
41  use eos_polar, only: f_polar, phi_polar
42  USE eos_numerical_derivatives
43 
44  IMPLICIT NONE
45 
46  !-----------------------------------------------------------------------------
47  INTEGER :: i, j, ik
48  INTEGER :: l, m
49  INTEGER :: kk, converg, read_file
50  INTEGER :: discret, fa(nc), outpt, irc(nc), irc_j, maxiter, ih
51  LOGICAL :: diagram
52  REAL, DIMENSION(-NDFT:NDFT) :: zp
53  REAL, DIMENSION(-NDFT:NDFT,2) :: rhop, rhop_o, r_corr
54  REAL, DIMENSION(-NDFT:NDFT) :: f_tot, f_corr
55  REAL :: f_l, f_r
56  REAL :: f_att
57  REAL, DIMENSION(-NDFT:NDFT,2) :: df_drho_tot
58  REAL, DIMENSION(2) :: df_drho_att
59  REAL, DIMENSION(2,0:nc) :: rhob
60  REAL, DIMENSION(nc) :: dhs_star
61  REAL :: sigij, dij
62  REAL :: f_disp_1pt, f_disp_pcsaft
63  REAL, DIMENSION(2) :: mu_disp_1pt, mu_disp_pcsaft
64  REAL :: zges(np), pbulk
65  REAL :: delta_st, tanhfac
66  REAL :: my0(nc)
67  REAL :: f_int_r, mu_int1_r, mu_int2_r
68  REAL :: zz1
69  REAL :: dz_local
70  REAL :: dev, maxdev
71  REAL :: delta_f, free
72  REAL :: surftens(0:240), st_macro(240)
73  REAL :: f01, f02, f03
74  ! REAL :: z3,z2,z1,z0,zms,z0_rk,z1_rk,z2_rk,mhc(2) ! temporary
75  ! REAL :: gij(2,2),gij_rk(2,2),dij_ab(2,2),tgij(2),lngij_rk(2) ! temporary
76  ! REAL :: mu_dsp(nc) ! temporary
77  INTEGER :: dum_i
78  REAL :: dum_r
79  CHARACTER*100 :: dum_text
80  REAL :: tc !, pc, rhoc, tsav, psav
81  REAL :: density(np), w(np,nc), lnphi(np,nc)
82  REAL :: damppic(2)
83  REAL :: box_l_no_unit
84  REAL, DIMENSION(-NDFT:NDFT, 2) :: lambda, rhobar
85  REAL, DIMENSION(-NDFT:NDFT) :: phi_dn0, phi_dn1, phi_dn2
86  REAL, DIMENSION(-NDFT:NDFT) :: phi_dn3, phi_dn4, phi_dn5
87  REAL, DIMENSION(0:5,-NDFT:NDFT) :: ni
88  REAL :: zs
89  REAL :: term1(nc), term2, i2_up(nc,nc), last_z = 0.0
90  REAL :: f_ch, df_drho_ch(nc)
91  REAL :: f_fmt, df_drho_fmt(nc)
92  REAL :: tc_l, tc_v, xi_save(np,nc), m_average
93  REAL :: f_assoc, mu_assoc(nc)
94  CHARACTER (LEN=4) :: char_len
95  ! REAL :: steps, end_x
96  REAL :: fres_polar, fdd, fqq, fdq
97  REAL :: mu_polar(nc), fdd_rk, fqq_rk, fdq_rk, z3_rk
98  !-----------------------------------------------------------------------------
99  ! note: variable ensemble_flag = 'tp' is set for calculating: mu_res=mu_res(T,p)
100  ! note: variable ensemble_flag = 'tv' is set for calculating: mu_res=mu_res(T,rho)
101  ! note: variable subtract1 = 'no' is set for regular calculations
102  ! note: variable subtract1 = '1PT' is set for calculating properties subtracting a
103  ! first order perturbation term (dispersion)
104 
105  num = 1
106  CALL set_default_eos_numerical
107 
108 
109  wca = .false.
110  shift = .false.
111  diagram = .true.
112 
113  CALL read_input
114 
115 
116  dhs(1:ncomp) = parame(1:ncomp,2) * ( 1.0 - 0.12*exp( -3.0*parame(1:ncomp,3)/t ) ) ! needed for rdf_matrix
117  dhs_star(1:ncomp) = dhs(1:ncomp)/parame(1:ncomp,2)
118  ! z3t_st = PI/6.0* parame(1,1) * d_hs**3
119 
120  IF ( ncomp > 2 ) THEN
121  write (*,*) 'SPECIFY ONLY ONE OR TWO COMPONENTS IN THE INPUT-FILE:'
122  write (*,*) ' ./input_file/INPUT.INP'
123  stop
124  END IF
125  OPEN (68,file = './output_file/DFT_profiles.xlo')
126  OPEN (69,file = './output_file/DFT_sigma.xlo')
127  OPEN (71,file = './output_file/DFT_iteration.xlo')
128  OPEN (72,file = './output_file/DFT_h.xlo')
129 
130 
131  !-----------------------------------------------------------------------------
132  ! The the cut-off distance rc is the max. distance for integrating inter-
133  ! actions. Beyond rc, tail-corrections are added.
134  !-----------------------------------------------------------------------------
135  rc = 9.0 ! dimensionless cut-off distance rc = r_c/sigma
136 
137 
138  !-----------------------------------------------------------------------------
139  ! basic definitions for calculating density profile,
140  !-----------------------------------------------------------------------------
141  ! grid-size dzp = zp(1)-zp(0)
142 
143  box_l_no_unit = 250.0 ! lenth of simulation box (Angstrom)
144  discret= 1000 ! number of spacial discretizations for profile
145  dzp = box_l_no_unit / REAL(discret) ! grid-distance (Angstrom)
146  fa(1:ncomp) = nint( parame(1:ncomp,2) / dzp + 1 ) ! number of steps per sigma (currently of component 1)
147  outpt = 100 ! number output-points = discret/outpt
148 
149 
150  !-----------------------------------------------------------------------------
151  ! definitions for the numerical algorithm
152  !-----------------------------------------------------------------------------
153 
154  maxiter = 200 ! maximum number of iterations per temp.step
155  maxdev = 1.e-6 ! maximum deviation
156  damppic(1) = 0.00001 ! damping of Picard-iteration
157  damppic(2) = 0.00001 ! damping of Picard-iteration
158 
159 
160  !-----------------------------------------------------------------------------
161  ! get a matrix of values for the pair correlation function g(r,rho)
162  !-----------------------------------------------------------------------------
163 
164  rg = 4.0
165  den_step = 40
166  !dzr = dzp / 2.0 ! dimensionless grid-distance for integrating
167  ! ! the Barker-Henderson attraction term
168  dzr = 0.02 ! dimensionless grid-distance for integrating
169  ! the Barker-Henderson attraction term
170  CALL rdf_matrix_mix
171 
172 
173  !-----------------------------------------------------------------------------
174  ! prepare for phase equilibrium calculation for given T
175  !-----------------------------------------------------------------------------
176 
177  ! Diagram_for_various_T_Loop: DO
178 
179  nphas = 2
180  outp = 0 ! output to terminal
181 
182 
183  !-----------------------------------------------------------------------------
184  ! calculate phase equilibrium for given T
185  !-----------------------------------------------------------------------------
186 
187  CALL start_var (converg) ! gets starting values, sets "val_init"
188 
189  IF ( converg /= 1 ) THEN
190  WRITE (*,*) 'no VLE found'
191  RETURN
192  END IF
193  ! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
194 !!$ n_unkw = ncomp ! number of quantities to be iterated
195 !!$ it(1) = 'p' ! iteration of temperature
196 !!$ it(2) = 'x21' ! iteration of mol fraction of comp.1 phase 2
197 !!$ sum_rel(1) = 'x12' ! summation relation: x12 = 1 - sum(x1j)
198 !!$ sum_rel(2) = 'x22' ! summation relation: x22 = 1 - sum(x2j)
199 !!$ end_x = 0.4
200 !!$ running = 'x11' ! xi(1,1) is running var. in PHASE_EQUILIB
201 !!$ steps = 5.0
202 !!$ ensemble_flag = 'tp' ! this flag is for: 'regular calculation'
203 !!$ CALL PHASE_EQUILIB( end_x, steps, converg )
204  ! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
205 
206  ! rhob(phase,0): molecular density
207  rhob(1,0) = dense(1) / ( pi/6.0* sum( xi(1,1:ncomp) * parame(1:ncomp,1) * dhs(1:ncomp)**3 ) )
208  rhob(2,0) = dense(2) / ( pi/6.0* sum( xi(2,1:ncomp) * parame(1:ncomp,1) * dhs(1:ncomp)**3 ) )
209  ! rhob(phase,i): molecular component density (with i=(1,...ncomp) ) in units (1/A^3)
210  rhob(1,1:ncomp) = rhob(1,0)*xi(1,1:ncomp)
211  rhob(2,1:ncomp) = rhob(2,0)*xi(2,1:ncomp)
212  WRITE (*,*) ' '
213  WRITE (*,*) 'temperature ',t, 'K, and p=', p/1.e5,' bar'
214  WRITE (*,*) 'x1_liquid ',xi(1,1),' x1_vapor', xi(2,1)
215  WRITE (*,*) 'densities ',rhob(1,0), rhob(2,0)
216  WRITE (*,*) 'dense ',dense(1), dense(2)
217 
218 
219  !-----------------------------------------------------------------------------
220  ! get density in SI-units (kg/m**3)
221  !-----------------------------------------------------------------------------
222  CALL si_dens ( density, w )
223 
224 
225  !-----------------------------------------------------------------------------
226  ! (re-)calculate residual chemical potential of both phases
227  !-----------------------------------------------------------------------------
228  ensemble_flag = 'tv' ! this flag is for: mu_res=mu_res(T,rho)
229  densta(1) = dense(1) ! Index 1 is for liquid density (here: packing fraction eta)
230  densta(2) = dense(2) ! Index 2 is for vapour density (here: packing fraction eta)
231  CALL fugacity (lnphi)
232  my0(1:ncomp) = lnphi(1,1:ncomp)! + LOG( rhob(1,1:ncomp) ) ! my0 = mu_res(T,rho_bulk_L) + ln(rho_bulk_l)
233  !zges(1) = p / ( RGAS*t*density(1) ) * mm_average/1000.0
234  !zges(2) = p / ( RGAS*t*density(2) ) * mm_average/1000.0
235  zges(1) = ( p * 1.e-30 ) / ( kbol*t* rhob(1,0) )
236  zges(2) = ( p * 1.e-30 ) / ( kbol*t* rhob(2,0) )
237 
238  pbulk = zges(1) * rhob(1,0) ! pressure p/kT (= Z*rho)
239  WRITE (*,*) ' '
240  WRITE (*,*) 'chem. potentials comp. 1', lnphi(1,1) + log( rhob(1,1) ), &
241  lnphi(2,1) + log( rhob(2,1) )
242  WRITE (*,*) 'chem. potentials comp. 2', lnphi(1,2) + log( rhob(1,2) ), &
243  lnphi(2,2) + log( rhob(2,2) )
244  WRITE (*,*) ' '
245 
246 
247  !-----------------------------------------------------------------------------
248  ! determine the critical temp. to xi of liquid and to xi of vapor
249  !-----------------------------------------------------------------------------
250 
251  num = 0
252  xi_save(:,:) = xi(:,:)
253  dense(1) = 0.15
254  WRITE (*,*) 'provide an estimate of the crit. Temp. of the mixture'
255  READ (*,*) t
256 
257  xif(1:ncomp) = xi_save(1,1:ncomp)
258  CALL heidemann_khalil
259  tc_l = t
260  WRITE (*,*) 'critical temperature to xi_liquid',tc_l
261 
262  xif(1:ncomp) = xi_save(2,1:ncomp)
263  ! dense(1) = 0.15
264  ! WRITE (*,*) 'provide an estimate of the crit. Temp. of the mixture'
265  ! READ (*,*) t
266  CALL heidemann_khalil
267  tc_v = t
268  WRITE (*,*) 'critical temperature to xi_vapor ',tc_v
269  num = 1
270 
271  ! tc_L = 600.0
272  ! WRITE (*,*) 'I have tentitativly set tc=700 '
273  ! pause
274 
275 
276  !tc = ( tc_L + tc_V ) / 2.0
277  tc = tc_l
278  xi(:,:) = xi_save(:,:)
279  densta(1:nphas) = val_conv(1:nphas)
280  dense(1:nphas) = val_conv(1:nphas)
281  t = val_conv(3)
282  WRITE (*,*) 'estimate of critical temperature:',tc,'K'
283  WRITE (*,*) ' '
284 
285  !tc = ( xi(1,1)*630.39 + xi(1,2)*551.92 &
286  ! + xi(2,1)*630.39 + xi(2,2)*551.92 ) / 2.0 - 10.0
287  !tc = ( xi(1,1)*191.406 + xi(1,2)*133.68 &
288  ! + xi(2,1)*191.406 + xi(2,2)*133.68 ) / 2.0
289 
290 
291  !-----------------------------------------------------------------------------
292  ! update z3t, the T-dependent quantity that relates eta and rho, as eta = z3t*rho
293  !-----------------------------------------------------------------------------
294  CALL perturbation_parameter
295 
296 
297  !-----------------------------------------------------------------------------
298  ! define initial density profile rhop(i)
299  ! and dimensionless space coordinates zp(i)
300  !
301  ! discret : number of grid-points within "the box"
302  ! irc : number of grid-points extending the the box to left and right
303  ! the box is extended in order to allow for numerical integration
304  ! 'irc' is determined by the cut-off distance 'rc'
305  !-----------------------------------------------------------------------------
306 
307  irc(1:ncomp) = nint(rc*parame(1:ncomp,2)/dzp) + 1
308  tanhfac = -2.3625*t/tc + 2.4728
309  ! tanhfac = 2.7*(1.0 - t/tc) ! this parameterization was published (Gross, JCP, 2009)
310 
311  irc_j = maxval( irc(1:ncomp) )
312  DO j = 1, ncomp
313  DO i = -irc_j, (discret+irc_j)
314  zp(i) = REAL(i) * dzp
315  END DO
316  DO i = -irc_j, (discret+irc_j)
317  !rhop(i,j) = TANH(-(zp(i)-zp(INT(discret/2))) / parame(j,2) *tanhfac) * (rhob(1,j)-rhob(2,j))/2.0 &
318  ! + (rhob(1,j)+rhob(2,j))/2.0
319  rhop(i,j) = ( tanh(-(zp(i)-zp(int(discret/2))) / parame(j,2) *tanhfac) + 1.0 ) * rhob(1,j)/2.0 &
320  - ( tanh(-(zp(i)-zp(int(discret/2))) / parame(j,2) *tanhfac) - 1.0 ) * rhob(2,j)/2.0
321  ! rhop(i,j) = rhob(1,j)
322  END DO
323  END DO
324 
325 
326  !-----------------------------------------------------------------------------
327  ! Optional: read starting profile of densities from file: DFT_profile.xlo
328  !-----------------------------------------------------------------------------
329  write (*,*) ' read starting profile of rho(z) from file: DFT_profile.xlo? (0: no, 1: yes)'
330  read (*,*) read_file
331  IF (read_file == 1 ) THEN
332  READ (68,*) dum_text
333  READ (68,*) dum_text
334  READ (68,*) dum_text
335  READ (68,*) dum_text
336  DO i = 0, discret-1
337  READ(68,*) dum_i, dum_r, rhop(i,1), rhop(i,2)
338  END DO
339  rewind(68)
340  END IF
341 
342 
343  !-----------------------------------------------------------------------------
344  ! Initialize the DENS_INV subroutine
345  !-----------------------------------------------------------------------------
346 
347  nphas = 1
348  ensemble_flag = 'tv' ! this flag is for: mu_res=mu_res(T,rho)
349 
350 
351 
352  !=============================================================================
353  ! Start iterating the density profile
354  !=============================================================================
355 
356  kk = 1
357  ih = 85
358 
359  dft_convergence_loop: DO
360 
361  !--------------------------------------------------------------------------
362  ! Getting auxilliary quantities along the profile
363  !--------------------------------------------------------------------------
364 
365  CALL aux_chain_mix ( discret, fa, dzp, zp, rhop, rhobar, lambda )
366 
367  CALL fmt_dens_mix ( discret, fa, dzp, zp, rhop, ni, &
368  phi_dn0, phi_dn1, phi_dn2, phi_dn3, phi_dn4, phi_dn5 )
369 
370 
371  !--------------------------------------------------------------------------
372  ! Start loop for density profile
373  !--------------------------------------------------------------------------
374 
375  DO i = 0, discret
376 
377  !-----------------------------------------------------------------------
378  ! hard sphere contribution
379  ! f_fmt is the Helmholtz energy density F_FMT/(VkT) = F_FMT/(NkT)*rho
380  ! dF_drho_fmt is the functional derivative of F_FMT/(kT) to rho
381  !-----------------------------------------------------------------------
382  zs = 1.0 - ni(3,i)
383  f_fmt = - ni(0,i)*log(zs) + ni(1,i)*ni(2,i)/zs - ni(4,i)*ni(5,i)/zs &
384  + (ni(2,i)**3 -3.0*ni(2,i)*ni(5,i)*ni(5,i)) *(ni(3,i)+zs*zs*log(zs)) &
385  /36.0/pi/zs/zs/ni(3,i)**2
386 
387  CALL df_fmt_drho_mix ( i, fa, dzp, zp, phi_dn0, &
388  phi_dn1, phi_dn2, phi_dn3, phi_dn4, phi_dn5, df_drho_fmt )
389 
390 
391 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
392 !!$ z0 = SUM( rhop(i,1:ncomp) * PI/6.0 * parame(1:ncomp,1) )
393 !!$ z1 = SUM( rhop(i,1:ncomp) * PI/6.0 * parame(1:ncomp,1) * dhs(1:ncomp) )
394 !!$ z2 = SUM( rhop(i,1:ncomp) * PI/6.0 * parame(1:ncomp,1) * dhs(1:ncomp)**2 )
395 !!$ z3 = SUM( rhop(i,1:ncomp) * PI/6.0 * parame(1:ncomp,1) * dhs(1:ncomp)**3 )
396 !!$ z0_rk = PI/6.0 * mseg(1)
397 !!$ z1_rk = PI/6.0 * mseg(1) * dhs(1)
398 !!$ z2_rk = PI/6.0 * mseg(1) * dhs(1)*dhs(1)
399 !!$ z3_rk = PI/6.0 * mseg(1) * dhs(1)**3
400 !!$ zms = 1.0 - z3
401 !!$ write(*,*) f_fmt,SUM( rhop(i,1:ncomp)*parame(1:ncomp,1) )*(3.0*z1*z2/zms &
402 !!$ + z2**3 /z3/zms/zms + (z2**3 /z3/z3-z0)*LOG(zms) )/z0
403 !!$ write(*,*) 'error A',dF_drho_fmt(1)- 6.0/PI* ( 3.0*(z1_rk*z2+z1*z2_rk)/zms + 3.0*z1*z2*z3_rk/zms/zms &
404 !!$ + 3.0*z2*z2*z2_rk/z3/zms/zms + z2**3 *z3_rk*(3.0*z3-1.0)/z3/z3/zms**3 &
405 !!$ + ((3.0*z2*z2*z2_rk*z3-2.0*z2**3 *z3_rk)/z3**3 -z0_rk) *LOG(zms) &
406 !!$ + (z0-z2**3 /z3/z3)*z3_rk/zms )
407 !!$ z0_rk = PI/6.0 * mseg(2)
408 !!$ z1_rk = PI/6.0 * mseg(2) * dhs(2)
409 !!$ z2_rk = PI/6.0 * mseg(2) * dhs(2)*dhs(2)
410 !!$ z3_rk = PI/6.0 * mseg(2) * dhs(2)**3
411 !!$ write(*,*) 'error A',dF_drho_fmt(2)- 6.0/PI* ( 3.0*(z1_rk*z2+z1*z2_rk)/zms + 3.0*z1*z2*z3_rk/zms/zms &
412 !!$ + 3.0*z2*z2*z2_rk/z3/zms/zms + z2**3 *z3_rk*(3.0*z3-1.0)/z3/z3/zms**3 &
413 !!$ + ((3.0*z2*z2*z2_rk*z3-2.0*z2**3 *z3_rk)/z3**3 -z0_rk) *LOG(zms) &
414 !!$ + (z0-z2**3 /z3/z3)*z3_rk/zms )
415 !!$ pause
416 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
417 
418  !-----------------------------------------------------------------------
419  ! chain term
420  !-----------------------------------------------------------------------
421  CALL df_chain_drho_mix ( i, fa, dzp, zp, rhop, lambda, rhobar, f_ch, df_drho_ch )
422 
423 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
424 !!$ z0 = SUM( rhop(i,1:ncomp) * PI/6.0 * parame(1:ncomp,1) )
425 !!$ z1 = SUM( rhop(i,1:ncomp) * PI/6.0 * parame(1:ncomp,1) * dhs(1:ncomp) )
426 !!$ z2 = SUM( rhop(i,1:ncomp) * PI/6.0 * parame(1:ncomp,1) * dhs(1:ncomp)**2 )
427 !!$ z3 = SUM( rhop(i,1:ncomp) * PI/6.0 * parame(1:ncomp,1) * dhs(1:ncomp)**3 )
428 !!$ zms = 1.0 - z3
429 !!$ z1_rk = PI/6.0 * mseg(1) * dhs(1)
430 !!$ z2_rk = PI/6.0 * mseg(1) * dhs(1)*dhs(1)
431 !!$ z3_rk = PI/6.0 * mseg(1) * dhs(1)**3
432 !!$ DO l = 1, ncomp
433 !!$ DO j = 1, ncomp
434 !!$ dij_ab(l,j)=dhs(l)*dhs(j)/(dhs(l)+dhs(j))
435 !!$ gij(l,j) = 1.0/zms + 3.0*dij_ab(l,j)*z2/zms/zms + 2.0*(dij_ab(l,j)*z2)**2 /zms**3
436 !!$ gij_rk(l,j) = z3_rk/zms/zms &
437 !!$ + 3.0*dij_ab(l,j)*(z2_rk+2.0*z2*z3_rk/zms)/zms/zms &
438 !!$ + dij_ab(l,j)**2 *z2/zms**3 *(4.0*z2_rk+6.0*z2*z3_rk/zms)
439 !!$ if (l == j ) tgij(l) = gij(l,j)
440 !!$ if (l == j ) lngij_rk(l) = gij_rk(l,j) / gij(l,j)
441 !!$ END DO
442 !!$ END DO
443 !!$ write(*,*) f_ch,SUM( rhop(i,1:ncomp) *(1.0- mseg(1:ncomp)) *LOG(tgij(1:ncomp)) )
444 !!$
445 !!$ mhc(1) = ( 1.0-mseg(1)) * LOG( gij(1,1) ) + SUM( rhop(i,1:ncomp) * (1.0-mseg(1:ncomp)) * lngij_rk(1:ncomp) )
446 !!$ write(*,*) 'error 1', dF_drho_ch(1) - mhc(1)
447 !!$ z1_rk = PI/6.0 * mseg(2) * dhs(2)
448 !!$ z2_rk = PI/6.0 * mseg(2) * dhs(2)*dhs(2)
449 !!$ z3_rk = PI/6.0 * mseg(2) * dhs(2)**3
450 !!$ DO l = 1, ncomp
451 !!$ DO j = 1, ncomp
452 !!$ gij_rk(l,j) = z3_rk/zms/zms &
453 !!$ + 3.0*dij_ab(l,j)*(z2_rk+2.0*z2*z3_rk/zms)/zms/zms &
454 !!$ + dij_ab(l,j)**2 *z2/zms**3 *(4.0*z2_rk+6.0*z2*z3_rk/zms)
455 !!$ if (l == j ) lngij_rk(l) = gij_rk(l,j) / gij(l,j)
456 !!$ END DO
457 !!$ END DO
458 !!$ mhc(2) = ( 1.0-mseg(2)) * LOG( gij(2,2) ) + SUM( rhop(i,1:ncomp) * (1.0-mseg(1:ncomp)) * lngij_rk(1:ncomp) )
459 !!$ write(*,*) 'error 2',dF_drho_ch(2) - mhc(2)
460 !!$ pause
461 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
462 
463 
464  !-----------------------------------------------------------------------
465  ! dispersive attraction
466  !-----------------------------------------------------------------------
467  f_att = 0.0
468  df_drho_att(:) = 0.0
469  term2 = 0.0
470  i2_up(:,:) = 0.0
471  DO l = 1, ncomp
472  term1(l) = 0.0
473  DO m = 1, ncomp
474  f01 = 0.0
475  f02 = 0.0
476  f03 = 0.0
477  term2 = 0.0
478  sigij = sig_ij(l,m)
479  dij = ( dhs(l) + dhs(m) ) / 2.0
480  irc_j = ( irc(l) + irc(m) + 1 ) / 2
481  DO j = (i-irc_j), (i+irc_j) ! integration in z-coordinate
482  zz1 = abs( zp(j) - zp(i) ) ! distance z12 between 1 and 2
483  dz_local = dzp
484  IF ( zz1 < rc*sigij ) THEN
485  IF ( -( zp(j) - zp(i) ) >= rc*sigij - dzp ) dz_local = rc*sigij - zz1
486  CALL dft_rad_int_mix ( i, j, l, m, ih, zz1, rhop, sigij, dij, f_int_r, mu_int1_r, mu_int2_r )
487 
488  f_att = f_att + pi*mseg(l)*mseg(m) *uij(l,m)/t * rhop(i,l) * &
489  dz_local * ( rhop(j,m)* f_int_r + f01 ) / 2.0
490  term1(l)= term1(l) + 2.0*pi*mseg(l)*mseg(m) *uij(l,m)/t * &
491  dz_local * ( rhop(j,m)*mu_int1_r + f02 ) / 2.0
492  ! term2 = term2 + PI*mseg(l)*mseg(m) *uij(l,m)/t * rhop(i,l) * &
493  ! dz_local * ( rhop(j,m)* mu_int2_r + f03 ) / 2.0
494  i2_up(l,m) = i2_up(l,m) + dz_local * ( rhop(j,m)* mu_int2_r + f03 ) / 2.0
495  f01 = rhop(j,m)* f_int_r
496  f02 = rhop(j,m)* mu_int1_r
497  f03 = rhop(j,m)* mu_int2_r
498  last_z = zz1
499  ELSE
500  f01 = rc * sig_ij(l,m) * 4.0 * ( rc**-12 - rc**-6 ) * rhop(j,m)
501  f02 = rc * sig_ij(l,m) * 4.0 * ( rc**-12 - rc**-6 ) * rhop(j,m)
502  f03 = 0.0
503  END IF
504 
505  END DO
506  dz_local = rc*sigij - last_z
507  f_att = f_att + pi*mseg(l)*mseg(m) *uij(l,m)/t * rhop(i,l) * dz_local * f01
508  term1(l) = term1(l) + 2.0*pi*mseg(l)*mseg(m) *uij(l,m)/t * dz_local * f02
509 
510  END DO
511  END DO
512 
513  term2 = 0.0
514  DO l = 1, ncomp
515  DO m = 1, ncomp
516  term2 = term2 + pi * rhop(i,l) * mseg(l)*mseg(m) * uij(l,m)/t *i2_up(l,m)
517  END DO
518  END DO
519 
520 
521  DO l = 1, ncomp
522  df_drho_att(l) = term1(l) + term2 * pi / 6.0 * mseg(l) * dhs(l)**3
523  END DO
524 
525 
526 
527  !-----------------------------------------------------------------------
528  ! cut-off corrections
529  !-----------------------------------------------------------------------
530 
531  ! CALL cutoff_corrections_mix ( i, irc, rhop, rhob, f_att, dF_drho_att )
532 
533 
534 
535  !-----------------------------------------------------------------------
536  ! The integration of the DFT-integral can be done with cubic splines
537  !-----------------------------------------------------------------------
538  ! stepno = discret + irc*2
539  ! CALL SPLINE_PARA (dzp, intgrid, utri, stepno)
540  ! CALL SPLINE_INT (f_int_z, dzp, intgrid, utri, stepno)
541  ! f_int_z = f_int_z + rhob(1)*( 4.0/90.0*rc**-9 -1.0/3.0*rc**-3 )
542  ! f_int_z = f_int_z + rhob(2)*( 4.0/90.0*rc**-9 -1.0/3.0*rc**-3 )
543 
544 
545  !-----------------------------------------------------------------------
546  ! Calculation of 'dF_drho_att', which denotes non-local contributions
547  ! of d(F)/d(rho*). Apart from the dispersive attraction, it contains
548  ! non-local contributions of the hard-sphere term, and the chain term.
549  ! 'dF_drho_att' is defined as: dF_drho_att = d(F)/d(rho*) - (dF/drho*)_LDA
550  ! where LDA is for 'local density approximation', i.e. the values for
551  ! the local density.
552  !
553  ! F is the intrinsic Helmholtz energy
554  ! rho* = rhop(i) denotes the dimensionless density.
555  ! parame(1,3) is epsilon/k (LJ-energy parameter)
556  ! parame(1,2) is sigma (LJ-size parameter)
557  !
558  ! For a non-local second order term uncomment lines starting with '!2'
559  !-----------------------------------------------------------------------
560  ! remember: factor 2*PI because of the cylindrical coordinates
561 
562  ! dF_drho_att(1:ncomp) = 2.0*PI*parame(1,1)**2 * parame(1,3)/t * mu_int_z(1:ncomp)
563 
564 
565 
566  !-----------------------------------------------------------------------
567  ! make the dispersive attraction term consistent with PC-SAFT, by
568  ! adding the difference ( PC-SAFT - 1PT )_dispersion locally (LDA)
569  !-----------------------------------------------------------------------
570 
571  !*****************************************************************************
572  ! under construction
573  !*****************************************************************************
574  dense(1) = pi / 6.0 * sum( rhop(i,1:ncomp) * mseg(1:ncomp) * dhs(1:ncomp)**3 )
575  densta(1) = dense(1)
576  xi(1,1:ncomp) = rhop(i,1:ncomp) / sum( rhop(i,1:ncomp) )
577 
578  ensemble_flag = 'tv'
579 
580  call only_one_term_eos_numerical ( 'disp_term', 'PT_MIX ' )
581  call fugacity ( lnphi )
582  call restore_previous_eos_numerical
583  f_disp_1pt = fres * sum( rhop(i,1:ncomp) )
584  mu_disp_1pt(1:ncomp) = lnphi(1,1:ncomp)
585 
586  call only_one_term_eos_numerical ( 'disp_term', 'PC-SAFT ' )
587  call fugacity ( lnphi )
588  call restore_previous_eos_numerical
589  f_disp_pcsaft = fres * sum( rhop(i,1:ncomp) )
590  mu_disp_pcsaft(1:ncomp) = lnphi(1,1:ncomp)
591 
592  fres_polar = 0.0
593  mu_polar(:) = 0.0
594  IF ( sum( parame(1:ncomp,6) ) > 1.e-10 .OR. sum( parame(1:ncomp,7) ) > 1.e-10 ) THEN
595  eta = densta(1)
596  rho = sum(rhop(i,1:ncomp))
597  x(1:ncomp) = xi(1,1:ncomp)
598  call f_polar ( fdd, fqq, fdq )
599  fres_polar = ( fdd + fqq + fdq ) * sum( rhop(i,1:ncomp) )
600  DO ik = 1, ncomp
601  z3_rk = pi/6.0 * mseg(ik) * dhs(ik)**3
602  call phi_polar ( ik, z3_rk, fdd_rk, fqq_rk, fdq_rk )
603  mu_polar(ik) = fdd_rk + fqq_rk + fdq_rk
604  END DO
605  END IF
606 
607  f_assoc = 0.0
608  mu_assoc(:) = 0.0
609  IF ( sum( nint(parame(1:ncomp,12)) ) > 0) THEN
610  call only_one_term_eos_numerical ( 'hb_term ', 'TPT1_Chap' )
611  call fugacity ( lnphi )
612  call restore_previous_eos_numerical
613  f_assoc = fres * sum( rhop(i,1:ncomp) )
614  mu_assoc(1:ncomp) = lnphi(1,1:ncomp)
615  END IF
616 
617  ensemble_flag = 'tp'
618 
619  ! CALL mu_pert_theory_mix ( mu_dsp )
620  ! write (*,*) dF_drho_att(1) - mu_disp_1PT(1), dF_drho_att(1), mu_disp_1PT(1)
621  ! write (*,*) dF_drho_att(2) - mu_disp_1PT(2), dF_drho_att(2), mu_disp_1PT(2)
622  ! write (*,*) f_att - f_disp_1PT, f_att, f_disp_1PT
623  ! write (*,*) 'err dis a',dF_drho_att(1) - mu_dsp(1)
624  ! write (*,*) 'err dis b',dF_drho_att(2) - mu_dsp(2)
625  ! write (*,*) 'err dis 1',dF_drho_att(1) - mu_disp_1PT(1)
626  ! write (*,*) 'err dis 2',dF_drho_att(2) - mu_disp_1PT(2)
627  ! pause
628 
629  df_drho_att(1:ncomp)= df_drho_att(1:ncomp) + ( mu_disp_pcsaft(1:ncomp) - mu_disp_1pt(1:ncomp) ) &
630  + mu_assoc(1:ncomp) + mu_polar(1:ncomp)
631  f_att = f_att + ( f_disp_pcsaft - f_disp_1pt ) + f_assoc + fres_polar
632 
633  !*****************************************************************************
634  !*****************************************************************************
635 
636 
637  !-----------------------------------------------------------------------
638  ! collect the total Helmholtz energy density 'f_tot' and the functional
639  ! derivative to rhop(z) 'dF_drho_tot' (both are residual quantities)
640  !-----------------------------------------------------------------------
641 
642  df_drho_tot(i,1:ncomp) = df_drho_fmt(1:ncomp) + df_drho_ch(1:ncomp) + df_drho_att(1:ncomp)
643  f_tot(i) = f_fmt + f_ch + f_att
644 
645 
646  !-----------------------------------------------------------------------
647  ! on-the-fly report on local errors in the profile
648  !-----------------------------------------------------------------------
649  IF ( mod(i, outpt) == 0 )write(*,'(i5,4(G15.6))') i, rhop(i,1), rhop(i,2), &
650  my0(1)-df_drho_tot(i,1)+log( rhob(1,1)/rhop(i,1) ), &
651  my0(2)-df_drho_tot(i,2)+log( rhob(1,2)/rhop(i,2) )
652 
653 
654  END DO ! end of loop (i = 0, discret) along the profile
655 
656 
657  !--------------------------------------------------------------------------
658  ! correction for numerical inaccuracies in calculating the density profile
659  !--------------------------------------------------------------------------
660 
661  DO j = 1, ncomp
662  f_l = exp( my0(j) - df_drho_tot(10,j) +log( rhob(1,j)/rhop(10,j) ) )
663  f_r = exp( my0(j) - df_drho_tot(discret-10,j) +log( rhob(1,j)/rhop(discret-10,j) ))
664  DO i = 0, discret
665  r_corr(i,j) = ( tanh(-(zp(i)-zp(int(discret/2))) / parame(j,2) *tanhfac) + 1.0 ) * f_l/2.0 &
666  - ( tanh(-(zp(i)-zp(int(discret/2))) / parame(j,2) *tanhfac) - 1.0 ) * f_r/2.0
667  END DO
668  END DO
669 
670  !--------------------------------------------------------------------------
671  ! update the density profile using either a Picard-iteration scheme
672  ! (i.e. a direct substitution scheme, or a LDA - Inversion Procedure
673  ! similar to R.Evans (Bristol), which is done in function DENS_INVERS2.
674  !--------------------------------------------------------------------------
675  dev = 0.0
676  DO j = 1, ncomp
677  DO i = 0, discret
678  rhop_o(i,j) = rhop(i,j)
679  rhop(i,j) = rhob(1,j) * exp( my0(j) - df_drho_tot(i,j) ) /r_corr(i,j)
680  rhop(i,j) = rhop_o(i,j) + ( rhop(i,j) - rhop_o(i,j) )*damppic(j)
681  dev = dev + ( ( rhop(i,j) - rhop_o(i,j) )*parame(j,2)**3 / damppic(j) )**2
682  END DO
683  END DO
684  ! if ( kk == 4 ) damppic(1) = damppic(1) * 2.0
685  if ( kk == 8 ) damppic(:) = damppic(:) * 10.0
686  if ( kk == 16 ) damppic(:) = damppic(:) * 5.0
687  dev = dev / real(discret)
688 
689 
690  !--------------------------------------------------------------------------
691  ! correction for numerical inaccuracies in calculating the surface tension
692  !--------------------------------------------------------------------------
693 
694  f_l = f_tot(10) + sum( rhop(10,1:ncomp)*( log(rhop(10,1:ncomp)/rhob(1,1:ncomp))-1.0 ) ) &
695  - ( sum(rhop(10,1:ncomp)*my0(1:ncomp)) - pbulk)
696  f_r = f_tot(discret-10) + sum( rhop(discret-10,1:ncomp)*( log(rhop(discret-10,1:ncomp)/rhob(1,1:ncomp))-1.0 ) ) &
697  - ( sum(rhop(discret-10,1:ncomp)*my0(1:ncomp)) - pbulk)
698  j=1
699  DO i = 1, discret
700  f_corr(i) = ( tanh(-(zp(i)-zp(int(discret/2))) / parame(j,2) *tanhfac) + 1.0 ) * f_l/2.0 &
701  - ( tanh(-(zp(i)-zp(int(discret/2))) / parame(j,2) *tanhfac) - 1.0 ) * f_r/2.0
702  END DO
703 
704  !--------------------------------------------------------------------------
705  ! calculate surface tension
706  !--------------------------------------------------------------------------
707  free = 0.0
708  DO i = 1, discret
709 
710  !-----------------------------------------------------------------------
711  ! now add the ideal gas term
712  !-----------------------------------------------------------------------
713  f_tot(i) = f_tot(i) + sum( rhop(i,1:ncomp)*( log(rhop(i,1:ncomp)/rhob(1,1:ncomp))-1.0 ) )
714 
715  delta_f = f_tot(i) - ( sum(rhop(i,1:ncomp)*my0(1:ncomp)) - pbulk) - f_corr(i) ! all quantities .../(kT)
716  free = free + delta_f*dzp
717 
718  END DO
719  surftens(kk) = kbol * t *1.e20*1000.0 *free
720 
721 
722  !--------------------------------------------------------------------------
723  ! add an approximate capillary wave contribution to the surface tension
724  !--------------------------------------------------------------------------
725  m_average = sum( rhob(1,1:ncomp)*parame(1:ncomp,1) ) / rhob(1,0)
726  !m_average = 0.5 * ( SUM( rhob(1,1:ncomp)*parame(1:ncomp,1) ) / rhob(1,0) &
727  ! + SUM( rhob(2,1:ncomp)*parame(1:ncomp,1) ) / rhob(2,0) )
728  st_macro(kk) = surftens(kk) / ( 1.0 + 3.0/8.0/pi *t/tc &
729  * (1.0/2.55)**2 / (0.0674*m_average+0.0045) )
730 
731 
732  delta_st = 1.0
733  !IF ( kk > 1 ) delta_st = ABS( surftens(kk)-surftens(kk-1) ) / surftens(kk)
734 
735  WRITE (*,*) '-----------------------------------------------------------'
736  WRITE (*,*) ' # error-sum intrinsic ST total ST'
737  WRITE (*,'(i3,3F15.6)') kk, dev, surftens(kk), st_macro(kk)
738  WRITE (*,*) '-----------------------------------------------------------'
739  WRITE (71,'(i3,3G18.8)') kk, dev, surftens(kk), st_macro(kk)
740  kk = kk + 1
741 
742  !--------------------------------------------------------------------------
743  ! convergence criterion
744  !--------------------------------------------------------------------------
745  ! write (*,*) 'error measure',dev, delta_st
746  IF ( dev < maxdev .OR. delta_st < 1.e-8 ) THEN
747  EXIT dft_convergence_loop
748  ELSE
749  IF ( kk > maxiter ) THEN
750  WRITE(*,*)' no convergence in ',maxiter,' steps'
751  EXIT dft_convergence_loop
752  ENDIF
753  ENDIF
754 
755  ENDDO dft_convergence_loop
756 
757  !-----------------------------------------------------------------------------
758  ! write resulting density profile
759  !-----------------------------------------------------------------------------
760  WRITE(68,'(a,4(G18.8))') 't,p,x11,x21',t, val_conv(4)/1.e5, rhob(1,1)/sum(rhob(1,1:ncomp)), &
761  rhob(2,1)/sum(rhob(2,1:ncomp))
762  WRITE(68,'(a,2G18.8)') 'rho1,rho2',density(1), density(2)
763  WRITE(68,'(a,2G18.8)') 'gamma',surftens(kk-1), st_macro(kk-1)
764  WRITE(68,*) ' '
765  DO i = 0, discret
766  WRITE (char_len,'(I3)') 2*ncomp+1
767  WRITE (68,'(i6,'//char_len//'(G18.10))') i,zp(i)-zp(int(discret/2)), &
768  rhop(i,1:ncomp), -( df_drho_tot(i,1:ncomp)-my0(1:ncomp) )
769  ! IF ( MOD(i, outpt) == 0 ) WRITE(*,'(i6,3(f18.12))') i,zp(i)-zp(INT(discret/2)),rhop(i)
770  ! write (69,*) zp(i)*parame(1,2), pN_pT(i)/parame(1,2)**3 /1E-30*parame(1,3)*KBOL/1.E6 ! in MPa
771  END DO
772  WRITE (68,*) ' '
773 
774 
775 !!$ ! ---------------------------------------------------------------------
776 !!$ ! summarize results
777 !!$ ! ---------------------------------------------------------------------
778 !!$ WRITE (*,*) ' '
779 !!$ WRITE (*,*) 'SUMMARY FOR A SINGLE TEMPERATURE'
780 !!$ WRITE (*,*) ' '
781 !!$ WRITE (*,*) 'Temp. [K], Pressure [bar] ',t,val_conv(4)/1.E5
782 !!$ WRITE (*,*) 'Critical point Temp., Press. ',tc,pc/1.E5
783 !!$ WRITE (*,*) 'Density [kg/m**3] ',density(1),density(2)
784 !!$ WRITE (*,*) 'Dimensionless Density (rho*) ',rhob(1),rhob(2)
785 !!$ WRITE (*,*) 'Excess Grand Potential ',free
786 !!$ WRITE (*,*) 'Intrinsic Interf. Tension [mN/m] ',surftens(kk-1)
787 !!$ WRITE (*,*) 'Macroscop.Interf. Tension [mN/m] ',st_macro(kk-1)
788 !!$ WRITE (*,*) '============================================================'
789 !!$ WRITE (*,*) ' '
790 !!$ WRITE (69,'(9(f18.10))') t, val_conv(4)/1.E5, &
791 !!$ rhob(1),rhob(2),surftens(kk-1),st_macro(kk-1),free,dev
792 !!$
793 !!$ ! ---------------------------------------------------------------------
794 !!$ ! when calc. a phase diagram & diagram of surface tens., loop for higher T
795 !!$ ! ---------------------------------------------------------------------
796 !!$ ensemble_flag = 'tp' ! this flag is for 'regular calculations'
797 !!$ subtract1 = 'no' ! this flag is for 'regular calculations'
798 !!$ IF ( diagram ) THEN
799 !!$ IF ( (t+8.0) <= tc ) THEN
800 !!$ t = t + 5.0
801 !!$ IF ( (t+15.0) <= tc ) t = t + 5.0
802 !!$ IF ( (t+25.0) <= tc ) t = t + 10.0
803 !!$ IF ( (t+45.0) <= tc ) t = t + 20.0
804 !!$ nphas = 2
805 !!$ n_unkw = ncomp ! number of quantities to be iterated
806 !!$ it(1) = 'p' ! iteration of pressure
807 !!$ val_init(3) = t ! value of temperature
808 !!$ running = 't' ! T is running variable in PHASE_EQUILIB - here T=const.
809 !!$ end_x = t ! end temperature - here T=const.
810 !!$ steps = 1.0 ! number of steps of running var. - here 1, since T=const.
811 !!$ d_hs = parame(1,2) * ( 1.0 - 0.12*EXP( -3.0*parame(1,3)/t ) )
812 !!$ z3t_st = PI/6.0* parame(1,1) * d_hs**3
813 !!$ dhs_st = d_hs/parame(1,2)
814 !!$ CALL rdf_matrix_units
815 !!$ ELSE
816 !!$ EXIT Diagram_for_various_T_Loop
817 !!$ END IF
818 !!$ WRITE (69,'(7(f18.10))') tc, pc/1.E5, rhoc, rhoc, 0., 0., 0.
819 !!$ ELSE
820 !!$ EXIT Diagram_for_various_T_Loop
821 !!$ END IF
822 !!$
823 !!$ ENDDO Diagram_for_various_T_Loop
824 
825 END SUBROUTINE dft_nmft_mix
826 
827 
828 
829 
830 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
831 !
832 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
833 
834 SUBROUTINE aux_chain_mix ( discret, fa, dzp, zp, rhop, rhobar, lambda )
835 
836  USE basic_variables
837  USE eos_variables, ONLY: dhs
838  USE dft_module, ONLY: ndft
839  IMPLICIT NONE
840 
841  !-----------------------------------------------------------------------------
842  INTEGER, INTENT(IN) :: discret
843  INTEGER, INTENT(IN) :: fa(nc)
844  REAL, INTENT(IN) :: dzp
845  REAL, INTENT(IN) :: zp(-ndft:ndft)
846  REAL, INTENT(IN) :: rhop(-ndft:ndft,2)
847  ! REAL, INTENT(IN) :: rhob(2,0:nc)
848  REAL, INTENT(OUT) :: rhobar(-ndft:ndft,2)
849  REAL, INTENT(OUT) :: lambda(-ndft:ndft,2)
850 
851  !-----------------------------------------------------------------------------
852  INTEGER :: i, j, k, nn
853  REAL :: int1, int2, zz1, xl, xh, dz = 0.0
854  REAL, DIMENSION(100) :: y2, rx, ry1, ry2
855  !-----------------------------------------------------------------------------
856 
857  ! DO i=(-fa),(discret+fa)
858  ! lambda(i) = 0.0
859  ! rhobar(i) = 0.0
860  ! f01 = 0.0 ! = 0.75*rhop(i-fa)*( 1.**2 -1.**2 )
861  ! f02 = 0.5 *rhop(i-fa)
862  ! DO j=(i-fa+1),(i+fa) ! fa=NINT(1./dzp)
863  ! zz1 = ABS(zp(j)-zp(i)) ! distance z12 between 1 and 2
864  ! int_r= 0.75*rhop(j)*(1.0-zz1*zz1)
865  ! int_l= 0.5 *rhop(j)
866  ! rhobar(i) = rhobar(i) + dzp * ( int_r+f01 )/2.0
867  ! lambda(i) = lambda(i) + dzp * ( int_l+f02 )/2.0
868  ! f01=int_r
869  ! f02=int_l
870  ! ENDDO
871  ! write(69,*) i,rhobar(i)
872  ! ENDDO
873 
874 
875  DO k = 1, ncomp
876  ry1(:) = 0.0
877  ry2(:) = 0.0
878  rx(:) = 0.0
879  y2(:) = 0.0
880  DO i = (-fa(k)), (discret+fa(k))
881  nn = 1
882  DO j = (i-fa(k)), (i+fa(k)) ! fa(k) = NINT(sigma(k)/dzp)
883  IF ( zp(j+1) > (zp(i)-dhs(k)) .AND. (zp(i)-dhs(k)) >= zp(j) ) THEN
884  dz = zp(j+1) - (zp(i)-dhs(k))
885  rx(1) = 0.0
886  ry1(1) = 0.0 ! = 0.75*rhop(i-fa,k)*( d.**2 -d.**2 )
887  ry2(1) = 0.5/dhs(k)*rhop(j,k) + ((zp(i)-dhs(k))-zp(j))/dzp &
888  *( 0.5/dhs(k)*rhop(j+1,k) - 0.5/dhs(k)*rhop(j,k) )
889  ELSE IF ( zp(j) > (zp(i)-dhs(k)) .AND. zp(j) <= (zp(i)+dhs(k)) ) THEN
890  zz1 = abs( zp(j)-zp(i) ) ! distance z12 between 1 and 2
891  nn = nn + 1
892  ry1(nn) = 0.75/dhs(k)**3 * rhop(j,k) * (dhs(k)**2 -zz1*zz1)
893  ry2(nn) = 0.5/dhs(k) * rhop(j,k)
894  rx(nn) = rx(nn-1) + dz
895  dz = dzp
896  IF ( (zp(i)+dhs(k)) < zp(j+1) ) THEN
897  dz = (zp(i)+dhs(k)) - zp(j)
898  nn = nn + 1
899  rx(nn) = rx(nn-1) + dz
900  ry1(nn) = 0.0
901  ry2(nn) = 0.5/dhs(k)*rhop(j,k) + ((zp(i)+dhs(k))-zp(j))/dzp &
902  *( 0.5/dhs(k)*rhop(j+1,k) - 0.5/dhs(k)*rhop(j,k) )
903  END IF
904  END IF
905  END DO
906  xl = rx(1)
907  xh = rx(nn)
908  CALL spline ( rx(1:nn), ry1(1:nn), nn, 1.e30, 1.e30, y2(1:nn) )
909  CALL splint_integral ( rx(1:nn), ry1(1:nn), y2(1:nn), nn, xl, xh, int1 )
910  rhobar(i,k) = int1
911  if ( rhobar(i,k) < 0.0 ) then
912  rhobar(i,k) = 0.0
913  do j = 2, nn
914  rhobar(i,k) = rhobar(i,k) + (ry1(j)+ry1(j-1))/2.0 *(rx(j)-rx(j-1))
915  end do
916  end if
917  if ( rhobar(i,k) < 0.0 ) rhobar(i,k) = rhop(i,k)
918  CALL spline ( rx(1:nn), ry2(1:nn), nn, 1.e30, 1.e30, y2(1:nn) )
919  CALL splint_integral ( rx(1:nn), ry2(1:nn), y2(1:nn), nn, xl, xh, int2 )
920  lambda(i,k) = int2
921  if ( lambda(i,k) < 0.0 ) then
922  lambda(i,k) = 0.0
923  do j = 2, nn
924  lambda(i,k) = lambda(i,k) + (ry2(j)+ry2(j-1))/2.0 *(rx(j)-rx(j-1))
925  end do
926  end if
927  if ( lambda(i,k) < 0.0 ) lambda(i,k) = rhop(i,k)
928 
929  ! write (*,'(i5,i5,3G19.10)') k,i,rhob(1,k),rhob(2,k),lambda(i,k)
930 
931  ! if ( k == 1 .AND. lambda(i,k) < 0.5*rhob(2,k) ) write (*,*) 'warning: lambda too low',i,lambda(i,k),rhob(2,k)
932  ! if ( k == 1 .AND. lambda(i,k) < 0.5*rhob(2,k) ) lambda(i,k) = 0.5*rhob(2,k)
933  END DO
934  END DO
935 
936 END SUBROUTINE aux_chain_mix
937 
938 
939 
940 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
941 !
942 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
943 
944 SUBROUTINE fmt_dens_mix ( discret, fa, dzp, zp, rhop, ni, &
945  phi_dn0, phi_dn1, phi_dn2, phi_dn3, phi_dn4, phi_dn5 )
946 
947  USE parameters, ONLY: pi
948  USE basic_variables
949  USE eos_variables, ONLY: dhs
950  USE dft_module, ONLY: ndft
951  IMPLICIT NONE
952 
953  !-----------------------------------------------------------------------------
954  INTEGER, INTENT(IN) :: discret
955  INTEGER, INTENT(IN) :: fa(nc)
956  REAL, INTENT(IN) :: dzp
957  REAL, INTENT(IN) :: zp(-ndft:ndft)
958  REAL, INTENT(IN) :: rhop(-ndft:ndft,2)
959  REAL, INTENT(OUT) :: ni(0:5,-ndft:ndft)
960  REAL, INTENT(OUT) :: phi_dn0(-ndft:ndft)
961  REAL, INTENT(OUT) :: phi_dn1(-ndft:ndft)
962  REAL, INTENT(OUT) :: phi_dn2(-ndft:ndft)
963  REAL, INTENT(OUT) :: phi_dn3(-ndft:ndft)
964  REAL, INTENT(OUT) :: phi_dn4(-ndft:ndft)
965  REAL, INTENT(OUT) :: phi_dn5(-ndft:ndft)
966 
967  !-----------------------------------------------------------------------------
968  INTEGER :: i, j, k, nn, fa2
969  REAL :: int2, int3, int5
970  REAL :: zz1, zs, d2, xl, xh, dz = 0.0
971  REAL, DIMENSION(100) :: y2, hx, hy2, hy3, hy5
972  !-----------------------------------------------------------------------------
973 
974  ni(2,:) = 0.0 ! corresponds to z2*6
975  ni(3,:) = 0.0 ! corresponds to z3
976  ni(5,:) = 0.0
977  ni(0,:) = 0.0 ! corresponds to z0/PI*6.0
978  ni(1,:) = 0.0 ! corresponds to z1/PI*3.0
979  ni(4,:) = 0.0
980 
981  DO k = 1, ncomp
982  hx(:) = 0.0
983  hy2(:) = 0.0
984  hy3(:) = 0.0 ! = rhop(i-fa/2)*parame(1,1)*(0.25-(zp(i-fa/2)-zp(i))**2 )
985  hy5(:) = 0.0
986  DO i = (-fa(k)), (discret+fa(k))
987  nn = 1
988  d2 = dhs(k) / 2.0
989  fa2= ( fa(k) + 1 ) / 2
990  ! if (i==-fa(k)) write (*,*) 'even?',k,fa(k),fa2
991  ! if (i==-fa(k)) pause
992  DO j = (i-fa2), (i+fa2)
993  IF ( zp(j+1) > (zp(i)-d2) .AND. (zp(i)-d2) >= zp(j) ) THEN
994  dz = zp(j+1) - (zp(i)-d2)
995  hx(1) = 0.0
996  hy2(1) = rhop(j,k) *parame(k,1) + ((zp(i)-d2)-zp(j))/dzp &
997  * ( rhop(j+1,k)*parame(k,1) - rhop(j,k)*parame(k,1) )
998  hy3(1) = 0.0 ! = rhop(i-fa2)*parame(k,1)*(0.25-(zp(i-fa2)-zp(i))**2 )
999  hy5(1) = rhop(j,k)*parame(k,1)*(zp(j)-zp(i)) +((zp(i)-d2)-zp(j))/dzp &
1000  * ( rhop(j+1,k)*parame(k,1)*(zp(j+1)-zp(i)) - rhop(j,k) *parame(k,1)*(zp(j)-zp(i)) )
1001  ELSE IF ( zp(j) > (zp(i)-d2) .AND. zp(j) <= (zp(i)+d2) ) THEN
1002  zz1 = zp(j)-zp(i) ! distance z12 between 1 and 2
1003  nn = nn + 1
1004  hy2(nn) = rhop(j,k)*parame(k,1)
1005  hy3(nn) = rhop(j,k)*parame(k,1) * ( d2**2 - zz1**2 )
1006  hy5(nn) = rhop(j,k)*parame(k,1) * zz1
1007  hx(nn) = hx(nn-1) + dz
1008  dz = dzp
1009  IF ( (zp(i)+d2) < zp(j+1) ) THEN
1010  dz = (zp(i)+d2) - zp(j)
1011  nn = nn + 1
1012  hy2(nn) = rhop(j,k)*parame(k,1) + ((zp(i)+d2)-zp(j))/dzp &
1013  * ( rhop(j+1,k)*parame(k,1) - rhop(j,k)*parame(k,1) )
1014  hy3(nn) = 0.0
1015  hy5(nn) = rhop(j,k)*parame(k,1)*(zp(j)-zp(i)) +((zp(i)+d2)-zp(j))/dzp &
1016  * ( rhop(j+1,k)*parame(k,1)*(zp(j+1)-zp(i)) - rhop(j,k)*parame(k,1)*(zp(j) -zp(i)) )
1017  hx(nn) = hx(nn-1) + dz
1018  END IF
1019  END IF
1020  END DO
1021  xl = hx(1)
1022  xh = hx(nn)
1023  CALL spline ( hx(1:nn), hy2(1:nn), nn, 1.e30, 1.e30, y2(1:nn) )
1024  CALL splint_integral ( hx(1:nn), hy2(1:nn), y2(1:nn), nn, xl, xh, int2 )
1025  CALL spline ( hx(1:nn), hy3(1:nn), nn, 1.e30, 1.e30, y2(1:nn) )
1026  CALL splint_integral ( hx(1:nn), hy3(1:nn), y2(1:nn), nn, xl, xh, int3 )
1027  CALL spline ( hx(1:nn), hy5(1:nn), nn, 1.e30, 1.e30, y2(1:nn) )
1028  CALL splint_integral ( hx(1:nn), hy5(1:nn), y2(1:nn), nn, xl, xh, int5 )
1029  ni(2,i) = ni(2,i) + pi * dhs(k) *int2 ! corresponds to z2*6
1030  ni(3,i) = ni(3,i) + pi * int3 ! corresponds to z3
1031  ni(5,i) = ni(5,i) + 2.0* pi *int5
1032  ni(0,i) = ni(0,i) + (pi * dhs(k) *int2) / (pi*dhs(k)**2 ) ! corresponds to z0/PI*6.0
1033  ni(1,i) = ni(1,i) + (pi * dhs(k) *int2) / (2.0*pi*dhs(k)) ! corresponds to z1/PI*3.0
1034  ni(4,i) = ni(4,i) + (2.0* pi *int5) / (2.0*pi*dhs(k))
1035  END DO
1036  END DO
1037 
1038  ! derivatives of phi to ni
1039  DO i = (-maxval(fa(1:ncomp))), (discret+maxval(fa(1:ncomp)))
1040  zs = 1.0 - ni(3,i)
1041  phi_dn0(i) = - log(zs) !!!!
1042  phi_dn1(i) = ni(2,i)/zs
1043  phi_dn2(i) = ni(1,i)/zs + 3.0*(ni(2,i)*ni(2,i) - ni(5,i)*ni(5,i)) &
1044  * (ni(3,i)+zs*zs*log(zs))/(36.0*pi*ni(3,i)*ni(3,i)*zs*zs)
1045  phi_dn3(i) = ni(0,i)/zs + (ni(1,i)*ni(2,i)-ni(4,i)*ni(5,i))/zs/zs &
1046  - (ni(2,i)**3 - 3.0*ni(2,i)*ni(5,i)*ni(5,i)) &
1047  *( ni(3,i)*(ni(3,i)*ni(3,i)-5.0*ni(3,i)+2.0) &
1048  + 2.0*zs**3 *log(zs) ) / ( 36.0*pi*(ni(3,i)*zs)**3 )
1049  phi_dn4(i) = - ni(5,i)/zs
1050  phi_dn5(i) = - ni(4,i)/zs - 6.0*ni(2,i)*ni(5,i) &
1051  * (ni(3,i)+zs*zs*log(zs))/(36.0*pi*ni(3,i)*ni(3,i)*zs*zs)
1052  END DO
1053 
1054 END SUBROUTINE fmt_dens_mix
1055 
1056 
1057 
1058 
1059 
1060 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1061 ! SUBROUTINE DFT_RAD_INT_mix
1062 !
1063 ! This subroutine integrates the kernel of the perturbation theory in
1064 ! radial direction (more precisely in distance coordinate r^ hat). A
1065 ! non-mean field approach is taken, using a radial distribution function
1066 ! (currently at a Percus-Yevick level).
1067 !
1068 ! The first and second order contributions to the perturbation theory
1069 ! are calculated - although the second order contribution is rarely used.
1070 !
1071 ! ToDo: comment the subroutine and remove the GOTO constructs
1072 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1073 
1074 SUBROUTINE dft_rad_int_mix ( i, j, l, m, ih, zz_1, rhop, sigij, dij, f_int_r, my_int1_r, my_int2_r )
1075 
1076  USE dft_module
1077  USE eos_variables, ONLY: pi, ncomp, mseg, dhs
1078  IMPLICIT NONE
1079 
1080  !-----------------------------------------------------------------------------
1081  INTEGER, INTENT(IN) :: i
1082  INTEGER, INTENT(IN) :: j
1083  INTEGER, INTENT(IN) :: l
1084  INTEGER, INTENT(IN) :: m
1085  INTEGER, INTENT(IN OUT) :: ih
1086  REAL, INTENT(IN) :: zz_1
1087  REAL, INTENT(IN) :: rhop(-ndft:ndft,2)
1088  REAL, INTENT(IN) :: sigij
1089  REAL, INTENT(IN) :: dij
1090  REAL, INTENT(OUT) :: f_int_r
1091  REAL, INTENT(OUT) :: my_int1_r
1092  REAL, INTENT(OUT) :: my_int2_r
1093 
1094  !-----------------------------------------------------------------------------
1095  INTEGER :: k
1096 
1097  REAL :: zz1
1098  REAL :: dzr_local
1099  REAL :: fint0, fint1
1100  REAL :: myint1_0, myint1_1
1101  REAL :: myint2_0, myint2_1
1102  REAL :: dg_dz3, dg_dr
1103  REAL :: rad, xg, rdf, z3_bar, ua, rs
1104  REAL :: analytic1, tau_rs
1105  LOGICAL :: shortcut
1106  !-----------------------------------------------------------------------------
1107 
1108  shortcut = .true.
1109  fint0 = rc * tau_cut
1110  myint1_0 = rc * tau_cut
1111  myint2_0 = 0.0
1112 
1113  f_int_r = 0.0
1114  my_int1_r = 0.0
1115  my_int2_r = 0.0
1116 
1117  !-----------------------------------------------------------------------------
1118  ! for mixtures it is advantageous to write all distances in dimensionless
1119  ! form, e.g. r^hat = r^hat / sigma.
1120  ! For mixtures, sigij = ( sig_i + sig_j ) / 2.
1121  !-----------------------------------------------------------------------------
1122  zz1 = zz_1 / sigij ! here, dimensionless
1123 
1124  !-----------------------------------------------------------------------------
1125  ! this block only speeds up the integration
1126  !-----------------------------------------------------------------------------
1127  IF ( shortcut ) THEN
1128  rs = max( rg, zz1 ) ! +dzr
1129  IF ( rs > rc ) WRITE (*,*) 'error !!!!'
1130  analytic1 = 0.4*rs**-10 - 0.4*rc**-10 - rs**-4 + rc**-4
1131  f_int_r = f_int_r + analytic1
1132  my_int1_r = my_int1_r + analytic1
1133  IF ( rs == zz1 ) GO TO 10
1134  tau_rs = 4.0 * ( rs**-12 - rs**-6 )
1135  fint0 = rs * tau_rs
1136  myint1_0 = rs * tau_rs
1137  rad = rs ! the simple integration scheme: set to rc
1138  k = 0 + nint( (rc-rs)/dzr ) ! in simple scheme: set to 0
1139  ELSE
1140  rad = rc
1141  k = 0
1142  END IF
1143  !-----------------------------------------------------------------------------
1144  !
1145  !-----------------------------------------------------------------------------
1146  !rad = rc
1147  !k = 0
1148 
1149 
1150  DO WHILE ( rad /= max( 1.0, zz1 ) )
1151 
1152  dzr_local = dzr ! dzr is set in f_dft (dimensionless)
1153 
1154  IF ( rad - dzr_local <= max( 1.0, zz1 ) ) THEN
1155  dzr_local = rad - max( 1.0, zz1 )
1156  rad = max( 1.0, zz1 )
1157  ELSE
1158  rad = rad - dzr_local
1159  END IF
1160 
1161  k = k + 1
1162  ua = 4.0 * ( rad**-12 -rad**-6 )
1163  xg = rad / dij * sigij
1164  !rho_bar = 0.5 * ( rhop(i) + rhop(j) ) * sigij**3
1165  z3_bar = pi / 6.0 * sum( 0.5*( rhop(i,1:ncomp) + rhop(j,1:ncomp) ) * mseg(1:ncomp) * dhs(1:ncomp)**3 )
1166  rdf = 1.0
1167  dg_dz3 = 0.0
1168  IF ( rad <= rg ) THEN
1169  IF ( l == 1 .AND. m == 1 ) CALL bi_cub_spline (z3_bar,xg,ya_11,x1a_11,x2a_11,y1a_11,y2a_11,y12a_11, &
1170  c_bicub_11,rdf,dg_dz3,dg_dr,den_step,ih,k)
1171  IF ( l /= m ) CALL bi_cub_spline (z3_bar,xg,ya_12,x1a_12,x2a_12,y1a_12,y2a_12,y12a_12, &
1172  c_bicub_12,rdf,dg_dz3,dg_dr,den_step,ih,k)
1173  IF ( l == 2 .AND. m == 2 ) CALL bi_cub_spline (z3_bar,xg,ya_22,x1a_22,x2a_22,y1a_22,y2a_22,y12a_22, &
1174  c_bicub_22,rdf,dg_dz3,dg_dr,den_step,ih,k)
1175  END IF
1176 
1177  fint1 = rdf * rad * ua
1178  myint1_1 = rdf * rad * ua
1179  myint2_1 = dg_dz3 * rad * ua
1180 
1181  f_int_r = f_int_r + dzr_local * ( fint1 + fint0 ) / 2.0
1182  my_int1_r = my_int1_r + dzr_local * ( myint1_1 + myint1_0) / 2.0
1183  my_int2_r = my_int2_r + dzr_local * ( myint2_1 + myint2_0) / 2.0
1184 
1185  fint0 = fint1
1186  myint1_0 = myint1_1
1187  myint2_0 = myint2_1
1188 
1189  ENDDO
1190 
1191 
1192 10 CONTINUE
1193 
1194  ! analytic1 = 4.0/10.0*rc**-10 - rc**-4
1195  ! f_int_r = f_int_r + analytic1
1196  ! my_int1_r = my_int1_r + analytic1
1197 
1198  f_int_r = f_int_r * sigij*sigij
1199  my_int1_r = my_int1_r * sigij*sigij
1200  my_int2_r = my_int2_r * sigij*sigij
1201 
1202 END SUBROUTINE dft_rad_int_mix
1203 
1204 
1205 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1206 !
1207 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1208 
1209 SUBROUTINE f_pert_theory_mix ( fdsp )
1210 
1211  USE eos_variables, ONLY: nc, pi, ncomp, t, rho, eta, x, mseg, dhs, sig_ij, uij
1212  USE dft_module
1213  IMPLICIT NONE
1214 
1215  !-----------------------------------------------------------------------------
1216  REAL, INTENT(IN OUT) :: fdsp
1217 
1218  !-----------------------------------------------------------------------------
1219  INTEGER :: k, ih
1220  INTEGER :: l, m
1221  REAL :: z3
1222  REAL :: ua, ua_c, rm
1223  REAL, DIMENSION(nc,nc) :: i1
1224  REAL :: int10, int11
1225  REAL :: d_ij, dzr_local
1226  REAL :: rad, xg, rdf
1227  REAL :: dg_dz3, dg_dr
1228  ! REAL :: intgrid(0:5000),intgri2(0:5000), utri(5000),I1_spline
1229  !-----------------------------------------------------------------------------
1230 
1231  !-----------------------------------------------------------------------------
1232  ! constants
1233  !-----------------------------------------------------------------------------
1234  ua_c = 4.0 * ( rc**-12 - rc**-6 )
1235  rm = 2.0**(1.0/6.0)
1236 
1237  i1(:,:) = 0.0
1238 
1239  DO l = 1, ncomp
1240  DO m = 1, ncomp
1241 
1242  rad = rc
1243 
1244  int10 = rc * rc * ua_c
1245  ! intgrid(0)= int10
1246 
1247  k = 0
1248  ih = 85
1249 
1250  DO WHILE ( rad /= 1.0 )
1251 
1252  dzr_local = dzr
1253  IF ( rad - dzr_local <= 1.0 ) dzr_local = rad - 1.0
1254 
1255  rad = rad - dzr_local
1256 
1257  k = k + 1
1258 
1259  d_ij = 0.5*(dhs(l)+dhs(m)) / sig_ij(l,m) ! dimensionless effective hs-diameter d(T)/sig
1260  xg = rad / d_ij
1261  z3 = eta
1262  rdf = 1.0
1263  dg_dz3 = 0.0
1264  IF ( rad <= rg ) THEN
1265  IF ( l == 1 .AND. m == 1 ) CALL bi_cub_spline (z3,xg,ya_11,x1a_11,x2a_11,y1a_11,y2a_11,y12a_11, &
1266  c_bicub_11,rdf,dg_dz3,dg_dr,den_step,ih,k)
1267  IF ( l /= m ) CALL bi_cub_spline (z3,xg,ya_12,x1a_12,x2a_12,y1a_12,y2a_12,y12a_12, &
1268  c_bicub_12,rdf,dg_dz3,dg_dr,den_step,ih,k)
1269  IF ( l == 2 .AND. m == 2 ) CALL bi_cub_spline (z3,xg,ya_22,x1a_22,x2a_22,y1a_22,y2a_22,y12a_22, &
1270  c_bicub_22,rdf,dg_dz3,dg_dr,den_step,ih,k)
1271  END IF
1272 
1273  ua = 4.0 * ( rad**-12 - rad**-6 )
1274 
1275  int11 = rdf * rad * rad * ua
1276  i1(l,m) = i1(l,m) + dzr_local * ( int11 + int10 ) / 2.0
1277 
1278  int10 = int11
1279  ! intgrid(k)= int11
1280 
1281  END DO
1282 
1283  ! stepno = k
1284  ! CALL SPLINE_PARA (dzr,intgrid,utri,stepno)
1285  ! CALL SPLINE_INT (I1_spline,dzr,intgrid,utri,stepno)
1286 
1287 
1288  ! caution: 1st order integral is in F_EOS.f defined with negative sign
1289  !-----------------------------------------------------------------------
1290  ! cut-off corrections
1291  !-----------------------------------------------------------------------
1292  ! I1(l,m) = I1(l,m) + ( 4.0/9.0 * rc**-9 - 4.0/3.0 * rc**-3 )
1293  ! I2(l,m) = I2(l,m) + 16.0/21.0 * rc**-21 - 32.0/15.0 * rc**-15 + 16.0/9.0 * rc**-9
1294 
1295  END DO
1296  END DO
1297 
1298 
1299  fdsp = 0.0
1300  DO l = 1, ncomp
1301  DO m = 1, ncomp
1302  fdsp = fdsp + 2.0*pi*rho*x(l)*x(m)* mseg(l)*mseg(m)*sig_ij(l,m)**3 * uij(l,m)/t *i1(l,m)
1303  ! ( 2.0*PI*rho*I1*order1 - PI*rho*c1_con*m_mean*I2*order2 )
1304  END DO
1305  END DO
1306 
1307 
1308 !!$ IF (disp_term == 'PT1') THEN
1309 !!$ c1_con = 0.0
1310 !!$ I2 = 0.0
1311 !!$ ELSEIF (disp_term == 'PT2') THEN
1312 !!$ zms = 1.0 - z3
1313 !!$ c1_con = 1.0/ ( 1.0 + m_mean*(8.0*z3-2.0*z3**2 )/zms**4 &
1314 !!$ + (1.0 - m_mean)*( 20.0*z3 -27.0*z3**2 +12.0*z3**3 -2.0*z3**4 ) &
1315 !!$ /(zms*(2.0-z3))**2 )
1316 !!$ ELSE
1317 !!$ write (*,*) 'define the type of perturbation theory'
1318 !!$ stop
1319 !!$ END IF
1320 
1321 
1322 END SUBROUTINE f_pert_theory_mix
1323 
1324 
1325 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1326 !
1327 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1328 
1329 SUBROUTINE mu_pert_theory_mix ( mu_dsp )
1330 
1331  USE eos_variables, ONLY: nc, pi, ncomp, t, rho, eta, x, mseg, dhs, sig_ij, uij
1332  USE dft_module
1333  IMPLICIT NONE
1334 
1335  !-----------------------------------------------------------------------------
1336  REAL, INTENT(IN OUT) :: mu_dsp(nc)
1337 
1338  !-----------------------------------------------------------------------------
1339  INTEGER :: k, ih
1340  INTEGER :: l, m
1341  REAL :: z3
1342  REAL :: ua, ua_c, rm
1343  REAL, DIMENSION(nc,nc) :: i1, i2
1344  REAL :: int1_0, int1_1, int2_0, int2_1
1345  REAL :: d_ij, dzr_local
1346  REAL :: rad, xg, rdf
1347  REAL :: dg_dz3, dg_dr
1348  REAL :: term1(nc), term2
1349  ! REAL :: intgrid(0:5000),intgri2(0:5000), utri(5000),I1_spline
1350  !-----------------------------------------------------------------------------
1351 
1352  !-----------------------------------------------------------------------------
1353  ! constants
1354  !-----------------------------------------------------------------------------
1355  ua_c = 4.0 * ( rc**-12 - rc**-6 )
1356  rm = 2.0**(1.0/6.0)
1357 
1358  i1(:,:) = 0.0
1359  i2(:,:) = 0.0
1360 
1361  DO l = 1, ncomp
1362 
1363  term1(l) = 0.0
1364 
1365  DO m = 1, ncomp
1366 
1367  rad = rc
1368 
1369  int1_0 = rc * rc * ua_c
1370  int2_0 = 0.0
1371 
1372  k = 0
1373  ih = 85
1374 
1375  DO WHILE ( rad /= 1.0 )
1376 
1377  dzr_local = dzr
1378  IF ( rad - dzr_local <= 1.0 ) dzr_local = rad - 1.0
1379 
1380  rad = rad - dzr_local
1381  k = k + 1
1382 
1383  d_ij = 0.5*(dhs(l)+dhs(m)) / sig_ij(l,m) ! dimensionless effective hs-diameter d(T)/sig
1384  xg = rad / d_ij
1385  z3 = eta
1386  rdf = 1.0
1387  dg_dz3 = 0.0
1388  IF ( rad <= rg ) THEN
1389  IF ( l == 1 .AND. m == 1 ) CALL bi_cub_spline (z3,xg,ya_11,x1a_11,x2a_11,y1a_11,y2a_11,y12a_11, &
1390  c_bicub_11,rdf,dg_dz3,dg_dr,den_step,ih,k)
1391  IF ( l /= m ) CALL bi_cub_spline (z3,xg,ya_12,x1a_12,x2a_12,y1a_12,y2a_12,y12a_12, &
1392  c_bicub_12,rdf,dg_dz3,dg_dr,den_step,ih,k)
1393  IF ( l == 2 .AND. m == 2 ) CALL bi_cub_spline (z3,xg,ya_22,x1a_22,x2a_22,y1a_22,y2a_22,y12a_22, &
1394  c_bicub_22,rdf,dg_dz3,dg_dr,den_step,ih,k)
1395  END IF
1396 
1397  ua = 4.0 * ( rad**-12 - rad**-6 )
1398 
1399  int1_1 = rdf * rad * rad * ua
1400  int2_1 = dg_dz3 * rad * rad * ua
1401  i1(l,m) = i1(l,m) + dzr_local * ( int1_1 + int1_0 ) / 2.0
1402  i2(l,m) = i2(l,m) + dzr_local * ( int2_1 + int2_0 ) / 2.0
1403 
1404  int1_0 = int1_1
1405  int2_0 = int2_1
1406 
1407  term1(l) = term1(l) +4.0*pi*rho*x(m)* mseg(l)*mseg(m) *sig_ij(l,m)**3 *uij(l,m)/t* dzr_local*(int1_1+int1_0)/2.0
1408 
1409  END DO
1410 
1411  END DO
1412  END DO
1413 
1414 
1415  ! DO l = 1, ncomp
1416  ! term1(l) = 0.0
1417  ! DO m = 1, ncomp
1418  ! term1(l) = term1(l) + 4.0*PI*rho*x(m)* mseg(l)*mseg(m) * sig_ij(l,m)**3 * uij(l,m)/t *I1(l,m)
1419  ! END DO
1420  ! END DO
1421 
1422  term2 = 0.0
1423  DO l = 1, ncomp
1424  DO m = 1, ncomp
1425  term2 = term2 + 2.0*pi*rho*x(l) * rho*x(m)* mseg(l)*mseg(m) * sig_ij(l,m)**3 * uij(l,m)/t *i2(l,m)
1426  END DO
1427  END DO
1428 
1429  DO l = 1, ncomp
1430  mu_dsp(l) = term1(l) + term2 * pi/ 6.0 * mseg(l)*dhs(l)**3
1431  END DO
1432 
1433 END SUBROUTINE mu_pert_theory_mix
1434 
1435 
1436 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1437 !
1438 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1439 
1440 SUBROUTINE rdf_matrix_mix
1441 
1442  USE parameters, ONLY: pi
1443  USE eos_variables, ONLY: parame, dhs
1444  USE dft_module
1445  IMPLICIT NONE
1446 
1447  !-----------------------------------------------------------------------------
1448  INTEGER :: i, j, k = 0
1449  REAL :: rdf, rad, xg, rho_rdf, z3, mseg_ij, d_ij = 0.0
1450  !-----------------------------------------------------------------------------
1451 
1452 
1453  DO j = 1,3
1454 
1455  IF( j == 1) mseg_ij = 0.5 * ( parame(1,1) + parame(1,1) )
1456  IF( j == 2) mseg_ij = 0.5 * ( parame(1,1) + parame(2,1) )
1457  IF( j == 3) mseg_ij = 0.5 * ( parame(2,1) + parame(2,1) )
1458 
1459  DO i = 1,den_step
1460  ! rho_rdf= rhob(1) + (rhob(2)-rhob(1))*REAL(i)/den_step
1461  rho_rdf = 1.e-5 + (1.0)*REAL(i-1)/REAL(den_step-1) ! segment density, rho_s*sigma**3
1462  rad = rc
1463  k = 0
1464  !write (*,*) 'eta=',rho_rdf*mseg_ij * PI/6.0 * dhs_st**3
1465  DO WHILE ( rad - 1.e-8 > 0.95 )
1466  rad = rad - dzr
1467  k = k + 1
1468  IF( j == 1) d_ij = dhs(1) / parame(1,2)
1469  IF( j == 2) d_ij = 0.5*(dhs(1)+dhs(2)) / ( 0.5*(parame(1,2)+parame(2,2)) )
1470  IF( j == 3) d_ij = dhs(2) / parame(2,2)
1471 
1472  xg = rad / d_ij
1473  z3 = rho_rdf * pi/6.0 * d_ij**3 ! d_ij is the dimensionless effective diameter d*(T)=d(T)/sigma
1474 
1475  IF ( j == 1) ya_11(i,k) = 1.0
1476  IF ( j == 2) ya_12(i,k) = 1.0
1477  IF ( j == 3) ya_22(i,k) = 1.0
1478 
1479  IF ( xg <= rg .AND. z3 > 0.0 ) CALL rdf_int ( z3, mseg_ij, xg, rdf )
1480 
1481  IF ( xg <= rg .AND. z3 > 0.0 .AND. j == 1 ) ya_11(i,k) = rdf
1482  IF ( xg <= rg .AND. z3 > 0.0 .AND. j == 2 ) ya_12(i,k) = rdf
1483  IF ( xg <= rg .AND. z3 > 0.0 .AND. j == 3 ) ya_22(i,k) = rdf
1484  ! ya(i,k) = y(x1a(i), x2a(k)) with x1a: density-vector, x2a: r-vector
1485  IF ( j == 1) x1a_11(i) = z3
1486  IF ( j == 1) x2a_11(k) = xg
1487  IF ( j == 2) x1a_12(i) = z3
1488  IF ( j == 2) x2a_12(k) = xg
1489  IF ( j == 3) x1a_22(i) = z3
1490  IF ( j == 3) x2a_22(k) = xg
1491  END DO
1492  END DO
1493 
1494  if ( xg > 1.0 ) stop 'rdf_matrix_mix: 0.95*sigma is too high for lower bound'
1495 
1496  WRITE (*,*) ' done with calculating g(r)',dhs_st
1497 
1498  kmax = k
1499  IF( j == 1) CALL bicub_derivative ( ya_11, x1a_11, x2a_11, y1a_11, y2a_11, y12a_11, den_step, kmax )
1500  IF( j == 1) CALL bicub_c ( ya_11, x1a_11, x2a_11, y1a_11, y2a_11, y12a_11, c_bicub_11, den_step, kmax )
1501  IF( j == 2) CALL bicub_derivative ( ya_12, x1a_12, x2a_12, y1a_12, y2a_12, y12a_12, den_step, kmax )
1502  IF( j == 2) CALL bicub_c ( ya_12, x1a_12, x2a_12, y1a_12, y2a_12, y12a_12, c_bicub_12, den_step, kmax )
1503  IF( j == 3) CALL bicub_derivative ( ya_22, x1a_22, x2a_22, y1a_22, y2a_22, y12a_22, den_step, kmax )
1504  IF( j == 3) CALL bicub_c ( ya_22, x1a_22, x2a_22, y1a_22, y2a_22, y12a_22, c_bicub_22, den_step, kmax )
1505 
1506  END DO
1507 
1508 END SUBROUTINE rdf_matrix_mix
1509 
1510 
1511 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1512 !
1513 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1514 
1515 SUBROUTINE df_fmt_drho_mix ( i, fa, dzp, zp, phi_dn0, &
1516  phi_dn1, phi_dn2, phi_dn3, phi_dn4, phi_dn5, dF_drho_fmt )
1517 
1518  USE parameters, ONLY: pi
1519  USE basic_variables
1520  USE eos_variables, ONLY: dhs
1521  USE dft_module, ONLY: ndft
1522  IMPLICIT NONE
1523 
1524  !-----------------------------------------------------------------------------
1525  INTEGER, INTENT(IN) :: i
1526  INTEGER, INTENT(IN) :: fa(nc)
1527  REAL, INTENT(IN) :: dzp
1528  REAL, INTENT(IN) :: zp(-ndft:ndft)
1529  REAL, INTENT(IN) :: phi_dn0(-ndft:ndft)
1530  REAL, INTENT(IN) :: phi_dn1(-ndft:ndft)
1531  REAL, INTENT(IN) :: phi_dn2(-ndft:ndft)
1532  REAL, INTENT(IN) :: phi_dn3(-ndft:ndft)
1533  REAL, INTENT(IN) :: phi_dn4(-ndft:ndft)
1534  REAL, INTENT(IN) :: phi_dn5(-ndft:ndft)
1535  REAL, INTENT(OUT) :: df_drho_fmt(nc)
1536 
1537  !-----------------------------------------------------------------------------
1538  INTEGER :: j, k
1539  INTEGER :: nn, fa2
1540  REAL :: int0, int1, int2, int3, int4, int5
1541  REAL :: zz1, d2, xl, xh, dz = 0.0
1542  REAL, DIMENSION(100) :: y2, hx, hy0, hy1, hy2, hy3, hy4, hy5
1543  !-----------------------------------------------------------------------------
1544 
1545 
1546  DO k = 1, ncomp
1547 
1548  hx(:) = 0.0
1549  hy0(:) = 0.0
1550  hy1(:) = 0.0
1551  hy2(:) = 0.0
1552  hy3(:) = 0.0
1553  hy4(:) = 0.0
1554  hy5(:) = 0.0
1555  y2(:) = 0.0
1556 
1557  nn = 1
1558  d2 = dhs(k) / 2.0
1559  fa2 = ( fa(k) + 1 ) / 2
1560  DO j = (i-fa2), (i+fa2)
1561  IF ( zp(j+1) > (zp(i)-d2) .AND. (zp(i)-d2) >= zp(j) ) THEN
1562  hx(1) = 0.0
1563  hy0(1) = phi_dn0(j) + ((zp(i)-d2)-zp(j))/dzp *( phi_dn0(j+1) - phi_dn0(j) )
1564  hy1(1) = phi_dn1(j) + ((zp(i)-d2)-zp(j))/dzp *( phi_dn1(j+1) - phi_dn1(j) )
1565  hy2(1) = phi_dn2(j) + ((zp(i)-d2)-zp(j))/dzp *( phi_dn2(j+1) - phi_dn2(j) )
1566  hy3(1) = 0.0 ! = rhop(i-fa/2)*parame(1,1)*(0.25-(zp(i-fa/2)-zp(i))**2)
1567  hy4(1) = phi_dn4(j)*(zp(j)-zp(i)) +((zp(i)-d2)-zp(j))/dzp &
1568  * ( phi_dn4(j+1)*(zp(j+1)-zp(i)) - phi_dn4(j)*(zp(j)-zp(i)) )
1569  hy5(1) = phi_dn5(j)*(zp(j)-zp(i)) +((zp(i)-d2)-zp(j))/dzp &
1570  * ( phi_dn5(j+1)*(zp(j+1)-zp(i)) - phi_dn5(j)*(zp(j)-zp(i)) )
1571  dz = zp(j+1) - (zp(i)-d2)
1572  ELSE IF ( zp(j) > (zp(i)-d2).AND.zp(j) <= (zp(i)+d2) ) THEN
1573  zz1 = zp(j) - zp(i) ! distance z12 between 1 and 2
1574  nn = nn + 1
1575  hy0(nn) = phi_dn0(j)
1576  hy1(nn) = phi_dn1(j)
1577  hy2(nn) = phi_dn2(j)
1578  hy3(nn) = phi_dn3(j) * ( d2**2 - zz1**2 )
1579  hy4(nn) = phi_dn4(j) * zz1
1580  hy5(nn) = phi_dn5(j) * zz1
1581  hx(nn) = hx(nn-1) + dz
1582  dz = dzp
1583  IF ( (zp(i)+d2) < zp(j+1) ) THEN
1584  dz = (zp(i)+d2) - zp(j)
1585  nn = nn + 1
1586  hx(nn) = hx(nn-1) + dz
1587  hy0(nn) = phi_dn0(j) + ( (zp(i)+d2) - zp(j) ) / dzp * ( phi_dn0(j+1) - phi_dn0(j) )
1588  hy1(nn) = phi_dn1(j) + ( (zp(i)+d2) - zp(j) ) / dzp * ( phi_dn1(j+1) - phi_dn1(j) )
1589  hy2(nn) = phi_dn2(j) + ( (zp(i)+d2) - zp(j) ) / dzp * ( phi_dn2(j+1) - phi_dn2(j) )
1590  hy3(nn) = 0.0
1591  hy4(nn) = phi_dn4(j) * (zp(j)-zp(i)) + ((zp(i)+d2)-zp(j)) / dzp &
1592  * ( phi_dn4(j+1)*(zp(j+1)-zp(i)) -phi_dn4(j)*(zp(j)-zp(i)) )
1593  hy5(nn) = phi_dn5(j) * (zp(j)-zp(i)) + ((zp(i)+d2)-zp(j))/dzp &
1594  * ( phi_dn5(j+1)*(zp(j+1)-zp(i)) -phi_dn5(j)*(zp(j)-zp(i)) )
1595  END IF
1596  END IF
1597  END DO
1598  xl = hx(1)
1599  xh = hx(nn)
1600  CALL spline ( hx(1:nn), hy0(1:nn), nn, 1.e30, 1.e30, y2(1:nn) )
1601  CALL splint_integral( hx(1:nn), hy0(1:nn), y2(1:nn), nn, xl, xh, int0 )
1602  CALL spline ( hx(1:nn), hy1(1:nn), nn, 1.e30, 1.e30, y2(1:nn) )
1603  CALL splint_integral( hx(1:nn), hy1(1:nn), y2(1:nn), nn, xl, xh, int1 )
1604  CALL spline ( hx(1:nn), hy2(1:nn), nn, 1.e30, 1.e30, y2(1:nn) )
1605  CALL splint_integral( hx(1:nn), hy2(1:nn), y2(1:nn), nn ,xl, xh, int2 )
1606  CALL spline ( hx(1:nn), hy3(1:nn), nn, 1.e30, 1.e30, y2(1:nn) )
1607  CALL splint_integral( hx(1:nn), hy3(1:nn), y2(1:nn), nn, xl, xh, int3 )
1608  CALL spline ( hx(1:nn), hy4(1:nn), nn, 1.e30, 1.e30, y2(1:nn) )
1609  CALL splint_integral( hx(1:nn), hy4(1:nn), y2(1:nn), nn, xl, xh, int4 )
1610  CALL spline ( hx(1:nn), hy5(1:nn), nn, 1.e30, 1.e30, y2(1:nn) )
1611  CALL splint_integral( hx(1:nn), hy5(1:nn), y2(1:nn), nn, xl, xh, int5 )
1612  ! write (*,*) int3,int4,int5
1613  ! pause
1614  ! int0 = int0/d_hs
1615  ! int1 = 0.5*parame(1,2)*int1
1616  ! int2 = PI*d_hs*parame(1,2)**2 *int2
1617  ! int3 = PI*parame(1,2)**3 *int3
1618  ! int4 = parame(1,2)*int4
1619  ! int5 = 2.0*PI*parame(1,2)**2 *int5
1620  int0 = int0 / dhs(k)
1621  int1 = 0.5 * int1
1622  int2 = pi * dhs(k) * int2
1623  int3 = pi * int3
1624  int4 = int4 / dhs(k)
1625  int5 = 2.0 * pi * int5
1626  ! dF_drho_fmt=dF_hs/drho = dF_hs/drho_segment *m_mean
1627  df_drho_fmt(k) = (int0+int1+int2+int3+int4+int5)*parame(k,1)
1628  END DO
1629 
1630 END SUBROUTINE df_fmt_drho_mix
1631 
1632 
1633 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1634 ! SUBROUTINE dF_chain_drho_mix
1635 !
1636 ! dF/kT[rho] / d(rho_k(r)) of the chain term is here calculated as
1637 ! = + (m_k - 1) * ln(rho_k) - (m_k - 1)* [ ln( y_kk(rho^bar)*lambda_k ) - 1 ]
1638 ! - SUM_j[ (m_j - 1) * INTEGRAL[ d(ln(y_jj))/r(rho^bar_k) *rho_j * 0.75/d_k^3*(d_k^2-zp^2) dz ] ]
1639 ! - (m_k - 1) * INTEGRAL[ 1/lambda_k * rho_k * 0.5/d_k dz ]
1640 ! The second last line is abbreviated as I_ch1. The last line is
1641 ! abbreviated as - (m_k - 1)*I_ch1.
1642 !
1643 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1644 
1645 SUBROUTINE df_chain_drho_mix ( i, fa, dzp, zp, rhop, lambda, rhobar, f_ch, dF_drho_ch )
1646 
1647  USE parameters, ONLY: pi
1648  USE basic_variables
1649  USE eos_variables, ONLY: dhs, mseg
1650  USE dft_module, ONLY: ndft
1651  IMPLICIT NONE
1652 
1653  !-----------------------------------------------------------------------------
1654  INTEGER, INTENT(IN) :: i
1655  INTEGER :: fa(nc)
1656  REAL, INTENT(IN) :: dzp
1657  REAL, INTENT(IN) :: zp(-ndft:ndft)
1658  REAL, INTENT(IN) :: rhop(-ndft:ndft,2)
1659  REAL, INTENT(IN) :: lambda(-ndft:ndft,2)
1660  REAL, INTENT(IN) :: rhobar(-ndft:ndft,2)
1661  REAL, INTENT(OUT) :: f_ch
1662  REAL, INTENT(OUT) :: df_drho_ch(nc)
1663 
1664  !-----------------------------------------------------------------------------
1665  INTEGER :: j, k, nn
1666  REAL :: zz1, xl, xh, dz = 0.0
1667  REAL, DIMENSION(100) :: y2, hx, hy0, hy1
1668  REAL, DIMENSION(2) :: ycorr
1669  REAL, DIMENSION(2,2) :: dlnydr
1670  REAL :: i_ch1, i_ch2
1671  !-----------------------------------------------------------------------------
1672 
1673  CALL cavity_mix ( rhobar(i,:), ycorr, dlnydr ) ! this is: d( ln( yij ) ) / d( rho(k) ) used with i=j
1674 
1675  f_ch = 0.0
1676 
1677  DO k = 1, ncomp
1678  hx(:) = 0.0
1679  hy0(:) = 0.0
1680  hy1(:) = 0.0
1681  nn = 1
1682  DO j = (i-fa(k)), (i+fa(k)) ! fa=NINT(1./dzp)
1683  IF ( zp(j+1) > (zp(i)-dhs(k)) .AND. zp(j)-1e-11 <= (zp(i)-dhs(k)) ) THEN
1684  dz = zp(j+1) - (zp(i)-dhs(k))
1685  hx(1) = 0.0
1686  hy0(1) = 0.0 ! = 0.75*rhop(i-fa)*DLNYDR(d_hs,rhobar(i-fa))*(d.**2 -d.**2 )
1687  hy1(1) = 0.5/dhs(k)*rhop(j,k)/lambda(j,k) + ((zp(i)-dhs(k))-zp(j))/dzp &
1688  * ( 0.5/dhs(k)*rhop(j+1,k)/lambda(j+1,k) - 0.5/dhs(k)*rhop(j,k)/lambda(j,k) )
1689  ELSE IF ( zp(j) > (zp(i)-dhs(k)) .AND. zp(j)+1.e-11 < (zp(i)+dhs(k)) ) THEN
1690  zz1 = zp(j) - zp(i) ! distance z12 between 1 and 2
1691  nn = nn + 1
1692 
1693  hy0(nn) = - sum( ( mseg(1:ncomp)-1.0 ) * rhop(j,1:ncomp)*dlnydr(1:ncomp,k) ) &
1694  * 0.75/dhs(k)**3 * (dhs(k)*dhs(k)-zz1*zz1)
1695  hy1(nn) = 0.5/dhs(k)*rhop(j,k)/lambda(j,k)
1696  hx(nn) = hx(nn-1) + dz
1697  dz = dzp
1698  IF ( zp(j+1)+1.e-11 >= (zp(i)+dhs(k)) ) THEN
1699  dz = (zp(i)+dhs(k)) - zp(j)
1700  nn = nn + 1
1701  hy0(nn) = 0.0
1702  hy1(nn) = 0.5/dhs(k)*rhop(j,k)/lambda(j,k) +((zp(i)+dhs(k))-zp(j))/dzp &
1703  * ( 0.5/dhs(k)*rhop(j+1,k)/lambda(j+1,k)-0.5/dhs(k)*rhop(j,k)/lambda(j,k) )
1704  hx(nn) = hx(nn-1) + dz
1705  END IF
1706  END IF
1707  END DO
1708  xl = hx(1)
1709  xh = hx(nn)
1710  CALL spline ( hx(1:nn), hy0(1:nn), nn, 1.e30, 1.e30, y2(1:nn) )
1711  CALL splint_integral( hx(1:nn), hy0(1:nn), y2(1:nn), nn, xl, xh, i_ch1 )
1712  CALL spline ( hx(1:nn), hy1(1:nn), nn, 1.e30, 1.e30, y2(1:nn) )
1713  CALL splint_integral( hx(1:nn), hy1(1:nn), y2(1:nn), nn, xl, xh, i_ch2 )
1714 
1715  f_ch = f_ch + ( parame(k,1) - 1.0 ) * rhop(i,k) * ( log(rhop(i,k)) - 1.0 ) &
1716  - ( parame(k,1) - 1.0 ) * rhop(i,k) * ( log(ycorr(k)*lambda(i,k)) - 1.0 )
1717  df_drho_ch(k) = + ( parame(k,1) - 1.0 ) * log(rhop(i,k)) &
1718  - ( parame(k,1) - 1.0 ) * ( log( ycorr(k)*lambda(i,k) ) - 1.0 + i_ch2 ) &
1719  + i_ch1
1720  END DO
1721 
1722 END SUBROUTINE df_chain_drho_mix
1723 
1724 
1725 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1726 !
1727 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1728 
1729 SUBROUTINE cavity_mix ( rhoi, ycorr, dlnydr )
1730 
1731  USE parameters, ONLY: pi
1732  USE basic_variables
1733  USE eos_variables, ONLY: dhs, mseg
1734  IMPLICIT NONE
1735 
1736  !-----------------------------------------------------------------------------
1737  REAL, INTENT(IN) :: rhoi(2)
1738  REAL, INTENT(OUT) :: ycorr(2)
1739  REAL, INTENT(OUT) :: dlnydr(2,2) ! this is: d( ln( yij ) ) / d( rho(k) ) used with i=j
1740 
1741  !-----------------------------------------------------------------------------
1742  INTEGER :: i, j, k
1743  REAL :: z0, z1, z2, z3, zms, z1_rk, z2_rk, z3_rk
1744  REAL, DIMENSION(nc,nc) :: dij_ab, gij, gij_rk
1745  !-----------------------------------------------------------------------------
1746 
1747  z0 = pi / 6.0 * sum( rhoi(1:ncomp) * mseg(1:ncomp) )
1748  z1 = pi / 6.0 * sum( rhoi(1:ncomp) * mseg(1:ncomp) * dhs(1:ncomp) )
1749  z2 = pi / 6.0 * sum( rhoi(1:ncomp) * mseg(1:ncomp) * dhs(1:ncomp)**2 )
1750  z3 = pi / 6.0 * sum( rhoi(1:ncomp) * mseg(1:ncomp) * dhs(1:ncomp)**3 )
1751 
1752  zms = 1.0 - z3
1753 
1754  DO i = 1,ncomp
1755  DO j=1,ncomp
1756  dij_ab(i,j)=dhs(i)*dhs(j)/(dhs(i)+dhs(j))
1757  ENDDO
1758  END DO
1759 
1760  DO k = 1, ncomp
1761  DO i = 1, ncomp
1762  z1_rk = pi/6.0 * mseg(k) * dhs(k)
1763  z2_rk = pi/6.0 * mseg(k) * dhs(k)*dhs(k)
1764  z3_rk = pi/6.0 * mseg(k) * dhs(k)**3
1765  !DO j = 1, ncomp
1766  j = i
1767  gij(i,j) = 1.0/zms + 3.0*dij_ab(i,j)*z2/zms/zms &
1768  + 2.0*(dij_ab(i,j)*z2)**2 /zms**3
1769  !dgijdz(i,j)= 1.0/zms/zms + 3.0*dij_ab(i,j)*z2*(1.0+z3)/z3/zms**3 &
1770  ! + (dij_ab(i,j)*z2/zms/zms)**2 *(4.0+2.0*z3)/z3
1771  gij_rk(i,j) = z3_rk/zms/zms &
1772  + 3.0*dij_ab(i,j)*(z2_rk+2.0*z2*z3_rk/zms)/zms/zms &
1773  + dij_ab(i,j)**2 *z2/zms**3 *(4.0*z2_rk+6.0*z2*z3_rk/zms)
1774  !END DO
1775 
1776  ycorr(i) = gij(i,i)
1777  dlnydr(i,k) = gij_rk(i,i) / gij(i,i)
1778 
1779  END DO
1780  END DO
1781 
1782 END SUBROUTINE cavity_mix
1783 
1784 
1785 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1786 !
1787 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1788 
1789 SUBROUTINE cutoff_corrections_mix ( i, irc, rhop, rhob, f_att, dF_drho_att )
1790 
1791  USE eos_variables, ONLY: nc, ncomp, pi, mseg, sig_ij, uij, t
1792  USE dft_module, ONLY: rc, ndft
1793  IMPLICIT NONE
1794 
1795  !-----------------------------------------------------------------------------
1796  INTEGER, INTENT(IN) :: i
1797  INTEGER, INTENT(IN) :: irc(nc)
1798  REAL, DIMENSION(-NDFT:NDFT,2), INTENT(IN) :: rhop
1799  REAL, DIMENSION(2,0:nc), INTENT(IN) :: rhob
1800  REAL, INTENT(IN OUT) :: f_att
1801  REAL, DIMENSION(2), INTENT(IN OUT) :: df_drho_att
1802 
1803  !-----------------------------------------------------------------------------
1804  INTEGER :: l, m
1805  INTEGER :: irc_j
1806  REAL :: cutoffz, rho_l, rho_m
1807  !-----------------------------------------------------------------------------
1808 
1809  cutoffz = 4.0/90.0 * rc**-9 - 1.0/3.0 * rc**-3
1810 
1811  !-----------------------------------------------------------------------------
1812  ! cutoff-corrections to the left
1813  !-----------------------------------------------------------------------------
1814  DO l = 1, ncomp
1815  DO m = 1, ncomp
1816  irc_j = ( irc(l) + irc(m) + 1 ) / 2
1817  rho_l = rhob(1,l)
1818  rho_m = rhob(1,m)
1819  IF ( ( abs(rhop(i-irc_j,l)-rhob(1,l))/rhob(1,l) ) > 1.e-3 ) rho_l = rhop(i-irc_j,l) ! a crude approximation !
1820  IF ( ( abs(rhop(i-irc_j,m)-rhob(1,m))/rhob(1,m) ) > 1.e-3 ) rho_m = rhop(i-irc_j,m) ! a crude approximation !
1821  write (*,*) 'le',i, rho_l,rhop(i,l)
1822  write (*,*) 'le',i, rho_m,rhop(i,m)
1823 
1824  f_att = f_att + pi*mseg(l)*mseg(m) *rho_l* rho_m *sig_ij(l,m) *uij(l,m)/t *cutoffz
1825  write(*,*) 'AAA', l,m,f_att,pi*mseg(l)*mseg(m) *rho_l* rho_m *sig_ij(l,m) *uij(l,m)/t *cutoffz
1826  df_drho_att(l) = df_drho_att(l) + 2.0*pi*mseg(l)*mseg(m) *rho_m *sig_ij(l,m) *uij(l,m)/t *cutoffz
1827  END DO
1828  END DO
1829 
1830  !-----------------------------------------------------------------------------
1831  ! cutoff-corrections to the right
1832  !-----------------------------------------------------------------------------
1833  DO l = 1, ncomp
1834  DO m = 1, ncomp
1835  irc_j = ( irc(l) + irc(m) + 1 ) / 2
1836  rho_l = rhob(2,l)
1837  rho_m = rhob(2,m)
1838  IF ( ( abs(rhop(i+irc_j,l)-rhob(2,l))/rhob(2,l) ) > 1.e-3 ) rho_l = rhop(i+irc_j,l) ! a crude approximation !
1839  IF ( ( abs(rhop(i+irc_j,m)-rhob(2,m))/rhob(2,m) ) > 1.e-3 ) rho_m = rhop(i+irc_j,m) ! a crude approximation !
1840 
1841  f_att = f_att + pi*mseg(l)*mseg(m) *rho_l* rho_m *sig_ij(l,m)**3 *uij(l,m)/t *cutoffz
1842  df_drho_att(l) = df_drho_att(l) + 2.0*pi*mseg(l)*mseg(m) *rho_m *sig_ij(l,m)**3 *uij(l,m)/t *cutoffz
1843  write(*,*) 'BBB', l,m,f_att,pi*mseg(l)*mseg(m) *rho_l* rho_m *sig_ij(l,m) *uij(l,m)/t *cutoffz
1844  END DO
1845  END DO
1846 
1847 END SUBROUTINE cutoff_corrections_mix
parameters
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains constant...
Definition: modules.f90:6
lambda
double lambda
Latent heat of blowing agent, J/kg.
Definition: experimentalInputs.h:118
dft_module
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW Module DFT_MODULE This module...
Definition: modules.f90:272
eos_variables
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:120
starting_values
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW Module STARTING_VALUES This m...
Definition: modules.f90:251
basic_variables
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:29
eos_numerical_derivatives
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:220
starting_values::start_var
subroutine, public start_var(converg)
subroutine start_var
Definition: starting_value.f90:65