MoDeNa  1.0
Software framework facilitating sequential multi-scale modelling
eos_f_contributions.f90
1 MODULE eos_f_contributions
2 
3  IMPLICIT NONE
4 
5  PRIVATE
6  PUBLIC :: f_ideal_gas, f_hard_sphere, f_chain_tpt1, f_chain_tpt_d, &
7  f_chain_hu_liu, f_chain_hu_liu_rc, f_spt, f_disp_pcsaft, &
8  f_association, f_ion_dipole_tbh, f_ion_ion_primmsa, &
9  f_ion_ion_nonprimmsa, f_lc_mayersaupe, f_pert_theory
10 
11 CONTAINS
12 
13 
14 
15 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
16 !
17 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
18 
19 SUBROUTINE f_ideal_gas ( fid )
20 
21  USE eos_variables, ONLY: nc, ncomp, x, rho, pi, kbol, nav
22 
23  !-----------------------------------------------------------------------------
24  REAL, INTENT(IN OUT) :: fid
25  !-----------------------------------------------------------------------------
26  INTEGER :: i
27  REAL, DIMENSION(nc) :: rhoi
28  !-----------------------------------------------------------------------------
29 
30  !h_Planck = 6.62606896E-34 ! Js
31  DO i = 1, ncomp
32  rhoi(i) = x(i) * rho
33  ! debroglie(i) = h_Planck *1d10 & ! in units Angstrom
34  ! *SQRT( 1.0 / (2.0*PI *1.0 / NAv / 1000.0 * KBOL*T) )
35  ! ! *SQRT( 1.0 / (2.0*PI *mm(i) /NAv/1000.0 * KBOL*T) )
36  ! fid = fid + x(i) * ( LOG(rhoi(i)*debroglie(i)**3) - 1.0 )
37  fid = fid + x(i) * ( log(rhoi(i)) - 1.0 )
38  END DO
39 
40 END SUBROUTINE f_ideal_gas
41 
42 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
43 !
44 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
45 
46 SUBROUTINE f_hard_sphere ( m_mean2, fhs )
47 
48  USE eos_variables, ONLY: z0t, z1t, z2t, z3t, eta, rho
49 
50  !-----------------------------------------------------------------------------
51  REAL, INTENT(IN) :: m_mean2
52  REAL, INTENT(IN OUT) :: fhs
53  !-----------------------------------------------------------------------------
54  REAL :: z0, z1, z2, z3, zms
55  !-----------------------------------------------------------------------------
56 
57  rho = eta / z3t
58  z0 = z0t * rho
59  z1 = z1t * rho
60  z2 = z2t * rho
61  z3 = z3t * rho
62  zms = 1.0 - z3
63 
64  fhs= m_mean2*( 3.0*z1*z2/zms + z2**3 /z3/zms/zms + (z2**3 /z3/z3-z0)*log(zms) )/z0
65 
66 
67 END SUBROUTINE f_hard_sphere
68 
69 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
70 !
71 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
72 
73 SUBROUTINE f_chain_tpt1 ( fhc )
74 
75  USE eos_variables, ONLY: nc, ncomp, mseg, x, z0t, z1t, z2t, z3t, &
76  rho, eta, dij_ab, gij
77 
78  !-----------------------------------------------------------------------------
79  REAL, INTENT(OUT) :: fhc
80  !-----------------------------------------------------------------------------
81  INTEGER :: i, j
82  REAL :: z0, z1, z2, z3, zms
83  !-----------------------------------------------------------------------------
84 
85  rho = eta / z3t
86  z0 = z0t * rho
87  z1 = z1t * rho
88  z2 = z2t * rho
89  z3 = z3t * rho
90  zms = 1.0 - z3
91 
92  DO i = 1, ncomp
93  DO j = 1, ncomp
94  gij(i,j) = 1.0/zms + 3.0*dij_ab(i,j)*z2/zms/zms + 2.0*(dij_ab(i,j)*z2)**2 / zms**3
95  END DO
96  END DO
97 
98  fhc = 0.0
99  DO i = 1, ncomp
100  fhc = fhc + x(i) *(1.0- mseg(i)) *log(gij(i,i))
101  END DO
102 
103 END SUBROUTINE f_chain_tpt1
104 
105 
106 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
107 !
108 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
109 
110 SUBROUTINE f_chain_tpt_d ( fhc )
111 
112  USE eos_variables, ONLY: nc, ncomp, mseg, x, z0t, z1t, z2t, z3t, rho, eta, &
113  mseg, dij_ab, gij
114 
115  !-----------------------------------------------------------------------------
116  REAL, INTENT(OUT) :: fhc
117  !-----------------------------------------------------------------------------
118  INTEGER :: i, j
119  REAL, DIMENSION(nc) :: gij_hd
120  REAL :: z0, z1, z2, z3, zms
121  !-----------------------------------------------------------------------------
122 
123  rho = eta / z3t
124  z0 = z0t * rho
125  z1 = z1t * rho
126  z2 = z2t * rho
127  z3 = z3t * rho
128  zms = 1.0 - z3
129 
130  DO i = 1, ncomp
131  DO j = 1, ncomp
132  gij(i,j) = 1.0/zms + 3.0*dij_ab(i,j)*z2/zms/zms + 2.0*(dij_ab(i,j)*z2)**2 / zms**3
133  END DO
134  END DO
135 
136  DO i = 1, ncomp
137  gij_hd(i) = 1.0/(2.0*zms) + 3.0*dij_ab(i,i)*z2 / zms**2
138  END DO
139 
140  fhc = 0.0
141  DO i = 1, ncomp
142  IF ( mseg(i) >= 2.0 ) THEN
143  fhc = fhc - x(i) * ( mseg(i)/2.0 * log( gij(i,i) ) + ( mseg(i)/2.0 - 1.0 ) * log( gij_hd(i)) )
144  ELSE
145  fhc = fhc + x(i) * ( 1.0 - mseg(i) ) * log( gij(i,i) )
146  END IF
147  END DO
148 
149 END SUBROUTINE f_chain_tpt_d
150 
151 
152 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
153 !
154 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
155 
156 SUBROUTINE f_chain_hu_liu ( fhc )
157 
158  USE eos_variables, ONLY: nc, ncomp, mseg, x, eta
159 
160  ! This subroutine calculates the hard chain contribution of the TPT-Liu-Hu Eos.
161  !-----------------------------------------------------------------------------
162  REAL, INTENT(OUT) :: fhc
163  !-----------------------------------------------------------------------------
164  REAL :: a2, b2, c2, a3, b3, c3
165  REAL :: a20, b20, c20, a30, b30, c30
166  REAL :: sum1, sum2, am, bm, cm
167  REAL :: zms
168  !-----------------------------------------------------------------------------
169 
170  zms = 1.0 - eta
171 
172  sum1 = sum( x(1:ncomp)*(mseg(1:ncomp)-1.0) )
173  sum2 = sum( x(1:ncomp)/mseg(1:ncomp)*(mseg(1:ncomp)-1.0)*(mseg(1:ncomp)-2.0) )
174 
175  a2 = 0.45696
176  a3 = -0.74745
177  b2 = 2.10386
178  b3 = 3.49695
179  c2 = 1.75503
180  c3 = 4.83207
181  a20 = - a2 + b2 - 3.0*c2
182  b20 = - a2 - b2 + c2
183  c20 = c2
184  a30 = - a3 + b3 - 3.0*c3
185  b30 = - a3 - b3 + c3
186  c30 = c3
187  am = (3.0 + a20) * sum1 + a30 * sum2
188  bm = (1.0 + b20) * sum1 + b30 * sum2
189  cm = (1.0 + c20) * sum1 + c30 * sum2
190 
191  fhc = - ( (am*eta - bm) / (2.0*zms) + bm/2.0/zms**2 - cm *log(zms) )
192 
193 
194 END SUBROUTINE f_chain_hu_liu
195 
196 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
197 !
198 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
199 
200 SUBROUTINE f_chain_hu_liu_rc ( fhs, fhc )
201 
202  USE eos_variables, ONLY: mseg, chir, eta
203 
204  ! This subroutine calculates the hard chain contribution of the TPT-Liu-Hu Eos.
205  !-----------------------------------------------------------------------------
206  REAL, INTENT(IN) :: fhs
207  REAL, INTENT(OUT) :: fhc
208  !-----------------------------------------------------------------------------
209  REAL :: a2, b2, c2, a3, b3, c3
210  REAL :: para1,para2,para3,para4
211  REAL :: alh,blh,clh
212  !-----------------------------------------------------------------------------
213 
214  ! This routine is only for pure components
215 
216  a2 = 0.45696
217  b2 = 2.10386
218  c2 = 1.75503
219 
220  para1 = -0.74745
221  para2 = 0.299154629727814
222  para3 = 1.087271036653154
223  para4 = -0.708979110326831
224  a3 = para1 + para2*chir(1) + para3*chir(1)**2 + para4*chir(1)**3
225  b3 = 3.49695 - (3.49695 + 0.317719074806190)*chir(1)
226  c3 = 4.83207 - (4.83207 - 3.480163780334421)*chir(1)
227 
228  alh = mseg(1)*(1.0 + ((mseg(1)-1.0)/mseg(1))*a2 + ((mseg(1)-1.0)/mseg(1))*((mseg(1)-2.0)/mseg(1))*a3 )
229  blh = mseg(1)*(1.0 + ((mseg(1)-1.0)/mseg(1))*b2 + ((mseg(1)-1.0)/mseg(1))*((mseg(1)-2.0)/mseg(1))*b3 )
230  clh = mseg(1)*(1.0 + ((mseg(1)-1.0)/mseg(1))*c2 + ((mseg(1)-1.0)/mseg(1))*((mseg(1)-2.0)/mseg(1))*c3 )
231 
232  fhc = ((3.0 + alh - blh + 3.0*clh)*eta - (1.0 + alh + blh - clh)) / (2.0*(1.0-eta)) + &
233  (1.0 + alh + blh - clh) / ( 2.0*(1.0-eta)**2 ) + (clh - 1.0)*log(1.0-eta)
234 
235  fhc = fhc - fhs
236 
237 END SUBROUTINE f_chain_hu_liu_rc
238 
239 
240 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
241 ! SUBROUTINE F_SPT
242 !
243 ! This routine calculates the residual Helmholtz energy contribution
244 ! obtained from Scaled Particle Theory (SPT).
245 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
246 
247 SUBROUTINE f_spt ( fhs, fhc )
248 
249  USE eos_variables, ONLY: nc, ncomp, pi, mseg, chir, eta
250 
251 
252  !-----------------------------------------------------------------------------
253  REAL, INTENT(OUT) :: fhc
254  REAL, INTENT(IN) :: fhs
255  !-----------------------------------------------------------------------------
256  real, dimension(nc,nc) :: c0, c1, c2
257 
258  INTEGER :: i,j
259  REAL :: msegav,aexv_1,bexv_1,chirav,aexv_2,aexv_3
260  REAL :: bexv_2,vm,b2,a_spt,sm,psi_spt
261  REAL, DIMENSION(nc,nc) :: b2iso
262  LOGICAL :: vex_williamson, vex_correlation
263  !-----------------------------------------------------------------------------
264 
265  !-----------------------------------------------------------------------------
266  ! Define which model to use for calculating the excluded volume
267  !-----------------------------------------------------------------------------
268  vex_williamson = .false.
269  vex_correlation = .true.
270 
271  ! ---------------------------------------------------------------------
272  ! Calculate isotropic 2nd virial coefficient B2iso
273  ! ---------------------------------------------------------------------
274 
275  IF ( vex_williamson ) THEN
276  DO i = 1, ncomp
277  DO j = 1, ncomp
278  msegav = ( mseg(i) + mseg(j) ) / 2.0
279  c0(i,j) = (pi/6.0)*(11.0*msegav - 3.0)
280  c1(i,j) = (pi/6.0)*3.53390 *(mseg(i)-1.0) * (mseg(j)-1.0)
281  c2(i,j) = 0.0
282  b2iso(i,j) = c0(i,j)/2.0 + (pi/8.0)*c1(i,j) + (1.0/3.0)*c2(i,j)
283  END DO
284  END DO
285  ELSEIF ( vex_correlation ) THEN
286  aexv_1 = 4.63
287  bexv_1 = 0.305
288  DO i = 1, ncomp
289  DO j = 1, ncomp
290  msegav = ( mseg(i) + mseg(j) ) / 2.0
291  chirav = ( chir(i) + chir(j) ) / 2.0
292  aexv_2 = -4.70 + 7.84/msegav
293  aexv_3 = 1.31 - 6.18/msegav
294  bexv_2 = -0.171 + 3.32/msegav
295 
296  c0(i,j) = (pi/6.0)* ( (11.0*msegav - 3.0) + (mseg(i) - 1.0)* &
297  (mseg(j) -1.0)*( aexv_1*(1.0-chirav) + aexv_2*(1.0-chirav)**2 &
298  + aexv_3*(1.0-chirav)**3 ) )
299  c1(i,j) = (pi/6.0)*3.53390*(mseg(i) - 1.0)*(mseg(j) -1.0)*chirav**2
300  c2(i,j) = (pi/6.0)*(mseg(i) - 1.0)*(mseg(j) -1.0)*( bexv_1*(1.0-chirav) &
301  + bexv_2*(1.0-chirav)**2 )
302  b2iso(i,j) = c0(i,j)/2.0 + (pi/8.0)*c1(i,j) + (1.0/3.0)*c2(i,j)
303  END DO
304  END DO
305  ELSE
306  write(*,*) 'no model for excluded volume defined'
307  stop
308  END IF
309 
310  !write(*,*)'f_spt'
311  !write(*,*)'C0, C1, C2 :',C0(1,1),C1(1,1),C2(1,1)
312 
313  ! This only holds for pure components!!!!!
314  vm = (pi/6.0) * mseg(1)
315  b2 = b2iso(1,1) / vm;
316  a_spt = (b2 - 1.0) / 3.0
317  sm = pi*mseg(1)
318  psi_spt = (sm/(9.0*vm)) * (3.0*a_spt - (sm/(4.0*vm))*(1.0 - (mseg(1)-1.0)*(1.0/4.0) ) )
319  fhc = (psi_spt -1.0)*log(1.0-eta) + 3.0*a_spt*eta/(1.0-eta) + psi_spt*eta/( (1.0-eta)**2 )
320 
321  fhc = fhc - fhs
322 
323 END SUBROUTINE f_spt
324 
325 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
326 !
327 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
328 
329 SUBROUTINE f_disp_pcsaft ( fdsp )
330 
331  USE eos_variables, ONLY: pi, rho, eta, z0t, apar, bpar, order1, order2
332 
333  !-----------------------------------------------------------------------------
334  REAL, INTENT(IN OUT) :: fdsp
335  !-----------------------------------------------------------------------------
336  INTEGER :: m
337  REAL :: i1, i2, c1_con, z3, zms, m_mean
338  !-----------------------------------------------------------------------------
339 
340  z3 = eta
341  zms = 1.0 - eta
342  m_mean = z0t / ( pi / 6.0 )
343 
344  i1 = 0.0
345  i2 = 0.0
346  DO m = 0, 6
347  i1 = i1 + apar(m) * z3**m
348  i2 = i2 + bpar(m) * z3**m
349  END DO
350 
351  c1_con= 1.0/ ( 1.0 + m_mean*(8.0*z3-2.0*z3**2 )/zms**4 &
352  + (1.0-m_mean)*( 20.0*z3-27.0*z3**2 +12.0*z3**3 -2.0*z3**4 )/(zms*(2.0-z3))**2 )
353 
354  fdsp = -2.0*pi*rho*i1*order1 - pi*rho*c1_con*m_mean*i2*order2
355 
356 END SUBROUTINE f_disp_pcsaft
357 
358 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
359 !
360 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
361 
362 !!$SUBROUTINE F_DISP_CKSAFT ( fdsp )
363 !!$
364 !!$ USE EOS_VARIABLES, ONLY: nc, ncomp, PI, TAU, t, rho, eta, x, z0t, mseg, vij, uij, parame, um
365 !!$ USE EOS_CONSTANTS, ONLY: DNM
366 !!$ IMPLICIT NONE
367 !!$
368 !!$ !-----------------------------------------------------------------------------
369 !!$ REAL, INTENT(IN OUT) :: fdsp
370 !!$ !-----------------------------------------------------------------------------
371 !!$ INTEGER :: i, j, n, m
372 !!$ REAL :: zmr, nmr, m_mean
373 !!$ !-----------------------------------------------------------------------------
374 !!$
375 !!$ m_mean = z0t / ( PI / 6.0 )
376 !!$
377 !!$ DO i = 1, ncomp
378 !!$ DO j = 1, ncomp
379 !!$ vij(i,j)=(0.5*((parame(i,2)*(1.0-0.12 *EXP(-3.0*parame(i,3)/t))**3 )**(1.0/3.0) &
380 !!$ +(parame(j,2)*(1.0-0.12 *EXP(-3.0*parame(j,3)/t))**3 )**(1.0/3.0)))**3
381 !!$ END DO
382 !!$ END DO
383 !!$ zmr = 0.0
384 !!$ nmr = 0.0
385 !!$ DO i = 1, ncomp
386 !!$ DO j = 1, ncomp
387 !!$ zmr = zmr + x(i)*x(j)*mseg(i)*mseg(j)*vij(i,j)*uij(i,j)
388 !!$ nmr = nmr + x(i)*x(j)*mseg(i)*mseg(j)*vij(i,j)
389 !!$ END DO
390 !!$ END DO
391 !!$ um = zmr / nmr
392 !!$ fdsp = 0.0
393 !!$ DO n = 1, 4
394 !!$ DO m = 1, 9
395 !!$ fdsp = fdsp + DNM(n,m) * (um/t)**n *(eta/TAU)**m
396 !!$ END DO
397 !!$ END DO
398 !!$ fdsp = m_mean * fdsp
399 !!$
400 !!$
401 !!$END SUBROUTINE F_DISP_CKSAFT
402 
403 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
404 !
405 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
406 
407 SUBROUTINE f_association ( eps_kij, k_kij, fhb )
408 
409  USE eos_variables, ONLY: nc, nsite, ncomp, t, z0t, z1t, z2t, z3t, rho, eta, x, &
410  parame, sig_ij, dij_ab, gij, nhb_typ, mx, nhb_no
411  USE utilities
412 
413  !-----------------------------------------------------------------------------
414  REAL, INTENT(IN) :: eps_kij, k_kij
415  REAL, INTENT(IN OUT) :: fhb
416  !-----------------------------------------------------------------------------
417  LOGICAL :: assoc
418  INTEGER :: i, j, k, l, no, ass_cnt, max_eval
419  REAL, DIMENSION(nc,nc) :: kap_hb
420  REAL, DIMENSION(nc,nc,nsite,nsite) :: eps_hb
421  REAL, DIMENSION(nc,nc,nsite,nsite) :: delta
422  REAL, DIMENSION(nc,nsite) :: mx_itr
423  REAL :: err_sum, sum0, attenu, tol, ass_s1, ass_s2
424  REAL :: z0, z1, z2, z3, zms
425  !-----------------------------------------------------------------------------
426 
427  assoc = .false.
428  DO i = 1,ncomp
429  IF (nint(parame(i,12)) /= 0) assoc = .true.
430  END DO
431  IF (assoc) THEN
432 
433  rho = eta / z3t
434  z0 = z0t * rho
435  z1 = z1t * rho
436  z2 = z2t * rho
437  z3 = z3t * rho
438  zms = 1.0 - z3
439 
440  DO i = 1, ncomp
441  DO j = 1, ncomp
442  gij(i,j) = 1.0/zms + 3.0*dij_ab(i,j)*z2/zms/zms + 2.0*(dij_ab(i,j)*z2)**2 / zms**3
443  END DO
444  END DO
445 
446 
447  DO i = 1, ncomp
448  IF ( nint(parame(i,12)) /= 0 ) THEN
449  nhb_typ(i) = nint( parame(i,12) )
450  kap_hb(i,i) = parame(i,13)
451  no = 0
452  DO k = 1,nhb_typ(i)
453  DO l = 1,nhb_typ(i)
454  eps_hb(i,i,k,l) = parame(i,(14+no))
455  no = no + 1
456  END DO
457  END DO
458  DO k = 1,nhb_typ(i)
459  nhb_no(i,k) = parame(i,(14+no))
460  no = no + 1
461  END DO
462  ELSE
463  nhb_typ(i) = 0
464  kap_hb(i,i)= 0.0
465  DO k = 1,nsite
466  DO l = 1,nsite
467  eps_hb(i,i,k,l) = 0.0
468  END DO
469  END DO
470  END IF
471  END DO
472 
473  DO i = 1,ncomp
474  DO j = 1,ncomp
475  IF ( i /= j .AND. (nhb_typ(i) /= 0.AND.nhb_typ(j) /= 0) ) THEN
476  ! kap_hb(i,j)= (kap_hb(i,i)+kap_hb(j,j))/2.0
477  ! kap_hb(i,j)= ( ( kap_hb(i,i)**(1.0/3.0) + kap_hb(j,j)**(1.0/3.0) )/2.0 )**3
478  kap_hb(i,j) = (kap_hb(i,i)*kap_hb(j,j))**0.5 &
479  *((parame(i,2)*parame(j,2))**3 )**0.5 &
480  / (0.5*(parame(i,2)+parame(j,2)))**3
481  kap_hb(i,j)= kap_hb(i,j)*(1.0-k_kij)
482  DO k = 1,nhb_typ(i)
483  DO l = 1,nhb_typ(j)
484  IF ( k /= l .AND. nhb_typ(i) >= 2 .AND. nhb_typ(j) >= 2 ) THEN
485  eps_hb(i,j,k,l) = (eps_hb(i,i,k,l)+eps_hb(j,j,l,k))/2.0
486  ! eps_hb(i,j,k,l) = (eps_hb(i,i,k,l)*eps_hb(j,j,l,k))**0.5
487  eps_hb(i,j,k,l) = eps_hb(i,j,k,l)*(1.0-eps_kij)
488  ELSE IF ( nhb_typ(i) == 1 .AND. l > k ) THEN
489  eps_hb(i,j,k,l) = (eps_hb(i,i,k,k)+eps_hb(j,j,l,k))/2.0
490  eps_hb(j,i,l,k) = (eps_hb(i,i,k,k)+eps_hb(j,j,l,k))/2.0
491  eps_hb(i,j,k,l) = eps_hb(i,j,k,l)*(1.0-eps_kij)
492  eps_hb(j,i,l,k) = eps_hb(j,i,l,k)*(1.0-eps_kij)
493  END IF
494  END DO
495  END DO
496  END IF
497  END DO
498  END DO
499 
500  !--------------------------------------------------------------------------
501  ! setting the self-association to zero for ionic compounds
502  !--------------------------------------------------------------------------
503  DO i = 1,ncomp
504  IF ( parame(i,10) /= 0) kap_hb(i,i)=0.0
505  DO j = 1,ncomp
506  IF ( parame(i,10) /= 0 .AND. parame(j,10) /= 0 ) kap_hb(i,j) = 0.0
507  END DO
508  END DO
509  ! kap_hb(1,2)=0.050
510  ! kap_hb(2,1)=0.050
511  ! eps_hb(2,1,1,1)=465.0
512  ! eps_hb(1,2,1,1)=465.0
513  ! nhb_typ(1) = 1
514  ! nhb_typ(2) = 1
515  ! nhb_no(1,1)= 1.0
516  ! nhb_no(2,1)= 1.0
517 
518 
519  DO i = 1, ncomp
520  DO k = 1, nhb_typ(i)
521  DO j = 1, ncomp
522  DO l = 1, nhb_typ(j)
523  delta(i,j,k,l) =gij(i,j) *kap_hb(i,j) *(exp(eps_hb(i,j,k,l)/t) - 1.0) *sig_ij(i,j)**3
524 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
525  ! IF ((i+j).EQ.3) delta(i,j,k,l)=94.0
526 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
527  END DO
528  END DO
529  IF ( mx(i,k) == 0.0 ) mx(i,k) = 1.0
530  END DO
531  END DO
532 
533  !--------------------------------------------------------------------------
534  ! constants for Iteration
535  !--------------------------------------------------------------------------
536  attenu = 0.7
537  tol = 1.e-10
538  IF (eta < 0.2) tol = 1.e-11
539  IF (eta < 0.01) tol = 1.e-14
540  max_eval = 200
541  err_sum = 2.0 * tol
542 
543  !--------------------------------------------------------------------------
544  ! Iterate over all components and all sites
545  !--------------------------------------------------------------------------
546  ass_cnt = 0
547  DO WHILE ( err_sum > tol .AND. ass_cnt <= max_eval )
548 
549  ass_cnt = ass_cnt + 1
550 
551  DO i = 1, ncomp
552  DO k = 1, nhb_typ(i)
553  sum0 = 0.0
554  DO j = 1, ncomp
555  DO l = 1, nhb_typ(j)
556  sum0 = sum0 + x(j) * mx(j,l) * nhb_no(j,l) * delta(i,j,k,l)
557  END DO
558  END DO
559  mx_itr(i,k) = 1.0 / (1.0 + sum0 * rho)
560  END DO
561  END DO
562 
563  err_sum = 0.0
564  DO i = 1, ncomp
565  DO k = 1, nhb_typ(i)
566  err_sum = err_sum + abs( mx_itr(i,k) - mx(i,k) ) ! / ABS(mx_itr(i,k))
567  mx(i,k) = mx_itr(i,k) * attenu + mx(i,k) * (1.0 - attenu)
568  IF ( mx(i,k) <= 0.0 ) mx(i,k) = 1.e-50
569  IF ( mx(i,k) > 1.0 ) mx(i,k) = 1.0
570  END DO
571  END DO
572 
573  END DO
574 
575  IF ( err_sum /= err_sum ) write (*,*) 'F_EOS: association NaN',ass_cnt, rho, sum0
576  IF ( err_sum /= err_sum ) read (*,*)
577  IF ( ass_cnt >= max_eval .AND. err_sum > sqrt(tol) ) THEN
578  WRITE(*,'(a,2G15.7)') 'F_EOS: Max_eval violated (mx) Err_Sum = ',err_sum,tol
579  END IF
580 
581  DO i = 1, ncomp
582  ass_s1 = 0.0
583  ass_s2 = 0.0
584  DO k = 1, nhb_typ(i)
585  ass_s1 = ass_s1 + nhb_no(i,k) * ( 1.0 - mx(i,k) )
586  ass_s2 = ass_s2 + nhb_no(i,k) * log(mx(i,k))
587  END DO
588  fhb = fhb + x(i) * ( ass_s2 + ass_s1 / 2.0 )
589  END DO
590 
591  END IF
592 
593 END SUBROUTINE f_association
594 
595 
596 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
597 !
598 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
599 
600 SUBROUTINE f_ion_dipole_tbh ( fhend )
601 
602  USE eos_variables, ONLY: nc, pi, kbol, nav, ncomp, t, rho, eta, x, z0t, parame, uij, sig_ij
603 
604  !-----------------------------------------------------------------------------
605  REAL, INTENT(IN OUT) :: fhend
606  !-----------------------------------------------------------------------------
607  INTEGER :: i, dipole, ions
608  REAL :: m_mean
609  REAL :: fh32, fh2, fh52, fh3
610  REAL :: e_elem, eps_cc0, rho_sol, dielec
611  REAL :: polabil, ydd, kappa, x_dipol, x_ions
612  REAL, DIMENSION(nc) :: my2dd, z_ii, e_cd, x_dd, x_ii
613  REAL :: sig_c, sig_d, sig_cd, r_s
614  REAL :: i0cc, i1cc, i2cc, icd, idd
615  REAL :: iccc, iccd, icdd, iddd
616  !-----------------------------------------------------------------------------
617 
618  m_mean = z0t / ( pi / 6.0 )
619 
620  !-----------------------------------------------------------------------------
621  ! Dieletric Constant of Water
622  !-----------------------------------------------------------------------------
623  e_elem = 1.602189246e-19 ! in Unit [Coulomb]
624  eps_cc0 = 8.854187818e-22 ! in Unit [Coulomb**2 /(J*Angstrom)]
625  ! Correlation of M. Uematsu and E. U. Frank
626  ! (Static Dieletric Constant of Water and Steam)
627  ! valid range of conditions 273,15 K <=T<= 823,15 K
628  ! and density <= 1150 kg/m3 (i.e. 0 <= p <= 500 MPa)
629  rho_sol = rho * 18.015 * 1.e27/ nav
630  rho_sol = rho_sol/1000.0
631  dielec = 1.0+(7.62571/(t/293.15))*rho_sol +(2.44e2/(t/293.15)-1.41e2 &
632  + 2.78e1*(t/293.15))*rho_sol**2 &
633  + (-9.63e1/(t/293.15)+4.18e1*(t/293.15) &
634  - 1.02e1*(t/293.15)**2 )*rho_sol**3 +(-4.52e1/(t/293.15)**2 &
635  + 8.46e1/(t/293.15)-3.59e1)*rho_sol**4
636 
637 
638  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
639 
640  dielec = 1.0
641 
642  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
643 
644 
645  !-----------------------------------------------------------------------------
646  ! Dipole-Ion Term
647  !-----------------------------------------------------------------------------
648  dipole = 0
649  ions = 0
650  fhend = 0.0
651  DO i = 1, ncomp
652  IF ( parame(i,6) /= 0.0 .AND. uij(i,i) /= 0.0 .AND. x(i) > 0.0 ) THEN
653  my2dd(i) = (parame(i,6))**2 *1.e-49 / (uij(i,i)*kbol*sig_ij(i,i)**3 *1.e-30)
654  dipole = 1
655  ELSE
656  my2dd(i) = 0.0
657  END IF
658 
659  z_ii(i) = parame(i,10)
660  IF ( z_ii(i) /= 0.0 .AND. uij(i,i) /= 0.0 .AND. x(i) > 0.0 ) THEN
661  e_cd(i) = ( parame(i,10)*e_elem* 1.e5 / sqrt(1.11265005) )**2 &
662  / ( uij(i,i)*kbol*sig_ij(i,i)*1.e-10 )
663  ions = 1
664  ELSE
665  e_cd(i) = 0.0
666  END IF
667  END DO
668 
669 
670  IF ( dipole == 1 .AND. ions == 1 ) THEN
671 
672  ydd = 0.0
673  kappa = 0.0
674  x_dipol = 0.0
675  x_ions = 0.0
676  polabil = 0.0
677  DO i = 1, ncomp
678  ydd = ydd + x(i)*(parame(i,6))**2 *1.e-49/ (kbol*t*1.e-30)
679  kappa = kappa + x(i) &
680  *(parame(i,10)*e_elem* 1.e5/sqrt(1.11265005))**2 /(kbol*t*1.e-10)
681  IF (parame(i,10) /= 0.0) THEN
682  x_ions = x_ions + x(i)
683  ELSE
684  polabil = polabil + 4.0*pi*x(i)*rho*1.4573 *1.e-30 &
685  / (sig_ij(3,3)**3 *1.e-30)
686  END IF
687  IF (parame(i,6) /= 0.0) x_dipol= x_dipol+ x(i)
688  END DO
689  ydd = ydd * 4.0/9.0 * pi * rho
690  kappa = sqrt( 4.0 * pi * rho * kappa )
691 
692  fh2 = 0.0
693  sig_c = 0.0
694  sig_d = 0.0
695  DO i=1,ncomp
696  x_ii(i) = 0.0
697  x_dd(i) = 0.0
698  IF(parame(i,10) /= 0.0 .AND. x_ions /= 0.0) x_ii(i) = x(i)/x_ions
699  IF(parame(i,6) /= 0.0 .AND. x_dipol /= 0.0) x_dd(i) = x(i)/x_dipol
700  sig_c = sig_c + x_ii(i)*parame(i,2)
701  sig_d = sig_d + x_dd(i)*parame(i,2)
702  END DO
703  sig_cd = 0.5 * (sig_c + sig_d)
704 
705  r_s = 0.0
706  ! DO i=1,ncomp
707  ! r_s=r_s + rho * x(i) * dhs(i)**3
708  ! END DO
709  r_s = eta*6.0 / pi / m_mean
710 
711  i0cc = - (1.0 + 0.97743 * r_s + 0.05257*r_s*r_s) &
712  /(1.0 + 1.43613 * r_s + 0.41580*r_s*r_s)
713  ! I1cc = - (10.0 - 2.0*z3 + z3*z3) /20.0/(1.0 + 2.0*z3)
714  i1cc = - (10.0 - 2.0*r_s*pi/6.0 + r_s*r_s*pi/6.0*pi/6.0) &
715  /20.0/(1.0 + 2.0*r_s*pi/6.0)
716 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
717  ! I2cc = + (z3-4.0)*(z3*z3+2.0) /24.0/(1.0+2.0*z3)
718  ! relation of Stell and Lebowitz
719  i2cc = -0.33331+0.7418*r_s - 1.2047*r_s*r_s &
720  + 1.6139*r_s**3 - 1.5487*r_s**4 + 0.6626*r_s**5
721 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
722  icd = (1.0 + 0.79576 *r_s + 0.104556 *r_s*r_s) &
723  /(1.0 + 0.486704*r_s - 0.0222903*r_s*r_s)
724  idd = (1.0 + 0.18158*r_s - 0.11467*r_s*r_s) &
725  /3.0/(1.0 - 0.49303*r_s + 0.06293*r_s*r_s)
726  iccc= 3.0*(1.0 - 1.05560*r_s + 0.26591*r_s*r_s) &
727  /2.0/(1.0 + 0.53892*r_s - 0.94236*r_s*r_s)
728  iccd= 11.0*(1.0 + 2.25642 *r_s + 0.05679 *r_s*r_s) &
729  /6.0/(1.0 + 2.64178 *r_s + 0.79783 *r_s*r_s)
730  icdd= 0.94685*(1.0 + 2.97323 *r_s + 3.11931 *r_s*r_s) &
731  /(1.0 + 2.70186 *r_s + 1.22989 *r_s*r_s)
732  iddd= 5.0*(1.0 + 1.12754*r_s + 0.56192*r_s*r_s) &
733  /24.0/(1.0 - 0.05495*r_s + 0.13332*r_s*r_s)
734 
735  IF ( sig_c <= 0.0 ) WRITE (*,*) 'error in Henderson ion term'
736 
737  fh32= - kappa**3 /(12.0*pi*rho)
738  fh2 = - 3.0* kappa**2 * ydd*icd /(8.0*pi*rho) / sig_cd &
739  - kappa**4 *sig_c/(16.0*pi*rho)*i0cc
740  IF (sig_d /= 0.0) fh2 = fh2 - 27.0* ydd * ydd*idd &
741  /(8.0*pi*rho) / sig_d**3
742  fh52= (3.0*kappa**3 * ydd + kappa**5 *sig_c**2 *i1cc) &
743  /(8.0*pi*rho)
744  fh3 = - kappa**6 * sig_c**3 /(8.0*pi*rho) *(i2cc-iccc/6.0) &
745  + kappa**4 * ydd *sig_c/(16.0*pi*rho) &
746  *( (6.0+5.0/3.0*sig_d/sig_c)*i0cc + 3.0*sig_d/sig_c*iccd ) &
747  + 3.0*kappa**2 * ydd*ydd /(8.0*pi*rho) / sig_cd &
748  *( (2.0-3.21555*sig_d/sig_cd)*icd +3.0*sig_d/sig_cd*icdd )
749  IF (sig_d /= 0.0) fh3 = fh3 + 27.0*ydd**3 &
750  /(16.0*pi*rho)/sig_d**3 *iddd
751 
752  fhend = ( fh32 + (fh32*fh32*fh3-2.0*fh32*fh2*fh52+fh2**3 ) &
753  /(fh2*fh2-fh32*fh52) ) &
754  / ( 1.0 + (fh32*fh3-fh2*fh52) /(fh2*fh2-fh32*fh52) &
755  - (fh2*fh3-fh52*fh52) /(fh2*fh2-fh32*fh52) )
756  !----------
757  ! fH32= - kappa**3 /(12.0*PI*rho)
758  ! fH2 = - 3.0* kappa**2 * ydd*Icd /(8.0*PI*rho) / sig_cd
759  ! fH52= (3.0*kappa**3 * ydd)/(8.0*PI*rho)
760  ! fH3 = + kappa**4 * ydd *sig_c/(16.0*PI*rho) &
761  ! *( (6.0+5.0/3.0*sig_d/sig_c)*0.0*I0cc + 3.0*sig_d/sig_c*Iccd) &
762  ! + 3.0*kappa**2 * ydd*ydd /(8.0*PI*rho) / sig_cd &
763  ! *( (2.0-3.215550*sig_d/sig_cd)*Icd +3.0*sig_d/sig_cd*Icdd )
764 
765  ! fHcd = ( + (fH32*fH32*fH3-2.0*fH32*fH2*fH52+fH2**3 ) &
766  ! /(fH2*fH2-fH32*fH52) ) &
767  ! / ( 1.0 + (fH32*fH3-fH2*fH52) /(fH2*fH2-fH32*fH52) &
768  ! - (fH2*fH3-fH52*fH52) /(fH2*fH2-fH32*fH52) )
769 
770  END IF
771 
772 END SUBROUTINE f_ion_dipole_tbh
773 
774 
775 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
776 !
777 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
778 
779 SUBROUTINE f_ion_ion_primmsa ( fcc )
780 
781  USE eos_variables, ONLY: nc, pi, kbol, nav, ncomp, t, rho, x, parame, mx
782 
783  !-----------------------------------------------------------------------------
784  REAL, INTENT(IN OUT) :: fcc
785  !-----------------------------------------------------------------------------
786  INTEGER :: i, j, cc_it, ions
787  REAL :: e_elem, eps_cc0, rho_sol, dielec
788  REAL :: x_ions
789  REAL :: cc_sig1, cc_sig2, cc_sig3
790  REAL, DIMENSION(nc) :: z_ii, x_ii, sigm_i, my2dd
791  REAL :: alpha_2, kappa, ii_par
792  REAL :: cc_omeg, p_n, q2_i, cc_q2, cc_gam
793  REAL :: cc_error(2), cc_delt
794  REAL :: rhs, lambda, lam_s
795  !-----------------------------------------------------------------------------
796 
797  !----------------Dieletric Constant of Water--------------------------
798  e_elem = 1.602189246e-19 ! in Unit [Coulomb]
799  eps_cc0 = 8.854187818e-22 ! in Unit [Coulomb**2 /(J*Angstrom)]
800  ! Correlation of M. Uematsu and E. U. Frank
801  ! (Static Dieletric Constant of Water and Steam)
802  ! valid range of conditions 273,15 K <=T<= 823,15 K
803  ! and density <= 1150 kg/m3 (i.e. 0 <= p <= 500 MPa)
804  rho_sol = rho * 18.015 * 1.e27/ nav
805  rho_sol = rho_sol/1000.0
806  dielec = 1.0+(7.62571/(t/293.15))*rho_sol +(2.44e2/(t/293.15)-1.41e2 &
807  +2.78e1*(t/293.15))*rho_sol**2 &
808  +(-9.63e1/(t/293.15)+4.18e1*(t/293.15) &
809  -1.02e1*(t/293.15)**2 )*rho_sol**3 +(-4.52e1/(t/293.15)**2 &
810  +8.46e1/(t/293.15)-3.59e1)*rho_sol**4
811 
812 
813  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
814 
815  dielec = 1.0
816 
817  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
818 
819 
820  !----------------Ion-Ion: primitive MSA -------------------------------
821  ! the (dipole moment)**2 [my**2] corresponds to an attraction from
822  ! point charges of [ SUM(xi * zi**2 * e_elem**2) * 3 * di**2 ]
823 
824  ! parame(ion,6))**2 * 1.E-49 / (kbol*T)
825  ! = (e_elem* 1.E5/SQRT(1.112650050))**2
826  ! *x(i)*zi**2 *3.0*sig_ij(1,1)**2 *1.E-20
827 
828  ! parame(ion,6))**2 = (e_elem* 1.E5/SQRT(1.112650050))**2 /1.E-49
829  ! *x(i)*zi**2 *3.0*sig_ij(i,i)**2 *1.E-20
830 
831  ! with the units
832  ! my**2 [=] D**2 = 1.E-49 J*m3
833  ! e_elem **2 [=] C**2 = 1.E5 / SQRT(1.112650050) J*m
834 
835 
836  ions = 0
837  x_ions = 0.0
838  fcc = 0.0
839  DO i = 1, ncomp
840  z_ii(i) = parame(i,10)
841  IF (z_ii(i) /= 0.0) THEN
842  sigm_i(i) = parame(i,2)
843  ELSE
844  sigm_i(i) = 0.0
845  END IF
846  IF (z_ii(i) /= 0.0) ions = 1
847  IF (z_ii(i) /= 0.0) x_ions = x_ions + x(i)
848  END DO
849 
850  IF (ions == 1 .AND. x_ions > 0.0) THEN
851 
852  cc_sig1 = 0.0
853  cc_sig2 = 0.0
854  cc_sig3 = 0.0
855  DO i=1,ncomp
856  IF (z_ii(i) /= 0.0) THEN
857  x_ii(i) = x(i)/x_ions
858  ELSE
859  x_ii(i) =0.0
860  END IF
861  cc_sig1 = cc_sig1 +x_ii(i)*sigm_i(i)
862  cc_sig2 = cc_sig2 +x_ii(i)*sigm_i(i)**2
863  cc_sig3 = cc_sig3 +x_ii(i)*sigm_i(i)**3
864  END DO
865 
866 
867  ! alpha_2 = 4.0*PI*e_elem**2 /eps_cc0/dielec/kbol/T
868  alpha_2 = e_elem**2 /eps_cc0 / dielec / kbol/t
869  kappa = 0.0
870  DO i = 1, ncomp
871  kappa = kappa + x(i)*z_ii(i)*z_ii(i)*mx(i,1)
872  END DO
873  kappa = sqrt( rho * alpha_2 * kappa )
874  ii_par= kappa * cc_sig1
875 
876  ! Temporaer: nach der Arbeit von Krienke verifiziert
877  ! noch nicht fuer Mischungen mit unterschiedl. Ladung erweitert
878  ! ii_par = SQRT( e_elem**2 /eps_cc0/dielec/kbol/T &
879  ! *rho*(x(1)*Z_ii(1)**2 + x(2)*Z_ii(2)**2 ) )*cc_sig1
880 
881 
882  cc_gam = kappa/2.0
883 
884  ! noch offen !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
885  cc_delt = 0.0
886  DO i = 1, ncomp
887  cc_delt = cc_delt + x(i)*mx(i,1)*rho*sigm_i(i)**3
888  END DO
889  cc_delt= 1.0 - pi/6.0*cc_delt
890 
891  cc_it = 0
892 13 CONTINUE
893  j = 0
894  cc_it = cc_it + 1
895 131 CONTINUE
896  j = j + 1
897  cc_omeg = 0.0
898  DO i = 1, ncomp
899  cc_omeg = cc_omeg +x(i)*mx(i,1)*sigm_i(i)**3 /(1.0+cc_gam*sigm_i(i))
900  END DO
901  cc_omeg = 1.0 + pi/2.0 / cc_delt * rho * cc_omeg
902  p_n = 0.0
903  DO i = 1, ncomp
904  p_n = p_n + x(i)*mx(i,1)*rho / cc_omeg*sigm_i(i)*z_ii(i) / (1.0+cc_gam*sigm_i(i))
905  END DO
906  q2_i = 0.0
907  cc_q2= 0.0
908  DO i = 1, ncomp
909  q2_i = q2_i + rho*x(i)*mx(i,1)*( (z_ii(i)-pi/2.0/cc_delt*sigm_i(i)**2 *p_n) &
910  /(1.0+cc_gam*sigm_i(i)) )**2
911  cc_q2 = cc_q2 + x(i)*mx(i,1)*rho*z_ii(i)**2 / (1.0+cc_gam*sigm_i(i))
912  END DO
913  q2_i = q2_i*alpha_2 / 4.0
914 
915  cc_error(j) = cc_gam - sqrt(q2_i)
916  IF (j == 1) cc_gam = cc_gam*1.000001
917  IF (j == 2) cc_gam = cc_gam - cc_error(2)* (cc_gam-cc_gam/1.000001)/(cc_error(2)-cc_error(1))
918 
919  IF ( j == 1 .AND. abs(cc_error(1)) > 1.e-15 ) GO TO 131
920  IF ( cc_it >= 10 ) THEN
921  WRITE (*,*) ' cc error'
922  stop
923  END IF
924  IF ( j /= 1 ) GO TO 13
925 
926  fcc= - alpha_2 / pi/4.0 /rho* (cc_gam*cc_q2 &
927  + pi/2.0/cc_delt *cc_omeg*p_n**2 ) + cc_gam**3 /pi/3.0/rho
928  ! Restricted Primitive Model
929  ! fcc=-(3.0*ii_par*ii_par+6.0*ii_par+2.0 &
930  ! -2.0*(1.0+2.0*ii_par)**1.50) &
931  ! /(12.0*PI*rho *cc_sig1**3 )
932 
933  ! fcc = x_ions * fcc
934 
935  my2dd(3) = (parame(3,6))**2 *1.e-19 /(kbol*t)
936  my2dd(3) = (1.84)**2 *1.e-19 /(kbol*t)
937 
938  rhs = 12.0 * pi * rho * x(3) * my2dd(3)
939  lam_s = 1.0
940 12 CONTINUE
941  lambda = (rhs/((lam_s+2.0)**2 ) + 16.0/((1.0+lam_s)**4 ) )**0.5
942  IF ( abs(lam_s-lambda) > 1.e-10 )THEN
943  lam_s = ( lambda + lam_s ) / 2.0
944  GO TO 12
945  END IF
946 
947  ! f_cd = -(ii_par*ii_par)/(4.0*PI*rho*m_mean *cc_sig1**3 ) &
948  ! *(dielec-1.0)/(1.0 + parame(3,2)/cc_sig1/lambda)
949  ! write (*,*) ' ',f_cd,fcc,x_ions
950  ! f_cd = f_cd/(1.0 - fcc/f_cd)
951  ! fcc = 0.0
952 
953  END IF
954 
955 
956 END SUBROUTINE f_ion_ion_primmsa
957 
958 
959 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
960 !
961 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
962 
963 SUBROUTINE f_ion_ion_nonprimmsa ( fdd, fqq, fdq, fcc )
964 
965  USE eos_variables, ONLY: nc, ncomp, x, parame, mseg
966 
967  !-----------------------------------------------------------------------------
968  REAL, INTENT(IN OUT) :: fdd, fqq, fdq, fcc
969  !-----------------------------------------------------------------------------
970  INTEGER :: dipole
971  !REAL :: A_MSA !, A_CC, A_CD, A_DD, U_MSA, chempot
972  REAL, DIMENSION(nc) :: x_export, msegm
973  !-----------------------------------------------------------------------------
974 
975  dipole = 0
976  IF ( sum( parame(1:ncomp,6) ) > 1.e-5 ) dipole = 1
977 
978  IF ( dipole /= 0 ) THEN ! alternatively ions and dipoles = 1
979  fdd = 0.0
980  fqq = 0.0
981  fdq = 0.0
982  fcc = 0.0
983  msegm(:) = mseg(:) ! the entries of the vector mseg and x are changed
984  x_export(:) = x(:) ! in SEMIRESTRICTED because the ions should be positioned first
985  ! that is why dummy vectors msegm and x_export are defined
986  !CALL SEMIRESTRICTED (A_MSA,A_CC,A_CD,A_DD,U_MSA, &
987  ! chempot,ncomp,parame,t,eta,x_export,msegm,0)
988  !fdd = A_MSA
989  write (*,*) 'why are individual contrib. A_CC,A_CD,A_DD not used'
990  stop
991  END IF
992 
993 END SUBROUTINE f_ion_ion_nonprimmsa
994 
995 
996 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
997 !
998 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
999 
1000 SUBROUTINE f_lc_mayersaupe ( flc )
1001 
1002  USE eos_variables, ONLY: nc, pi, kbol, nav, ncomp, phas, t, rho, eta, &
1003  x, mseg, e_lc, s_lc, dhs
1004 
1005  !-----------------------------------------------------------------------------
1006  REAL, INTENT(IN OUT) :: flc
1007  !-----------------------------------------------------------------------------
1008  INTEGER :: i, j, k
1009  INTEGER :: liq_crystal, count_lc, steps_lc
1010  REAL :: alpha_lc, tolerance, deltay
1011  REAL :: integrand1, integrand2, accel_lc
1012  REAL :: error_lc, u_term, sphase
1013  REAL, DIMENSION(nc) :: z_lc, s_lc1, s_lc2, sumu
1014  REAL, DIMENSION(nc,nc) :: u_lc, klc
1015  !-----------------------------------------------------------------------------
1016  INTEGER :: stabil
1017  COMMON /stabil / stabil
1018  !-----------------------------------------------------------------------------
1019 
1020 
1021  klc(1,2) = 0.0
1022  klc(2,1) = klc(1,2)
1023 
1024  alpha_lc = 1.0
1025  accel_lc = 4.0
1026  IF ( eta < 0.35 ) accel_lc = 1.3
1027  IF ( eta < 0.15 ) accel_lc = 1.0
1028 
1029  liq_crystal = 0
1030  DO i = 1, ncomp
1031  DO j = 1, ncomp
1032  e_lc(i,j) = (e_lc(i,i)*e_lc(j,j))**0.5 *(1.0-klc(i,j)) !combining rule
1033  ! E_LC(i,j)= ( E_LC(i,i)+E_LC(j,j) ) * 0.5 !combining rule
1034  ! S_LC(i,j)= ( S_LC(i,i)+S_LC(j,j) ) * 0.5 !combining rule
1035  IF (e_lc(i,j) /= 0.0) liq_crystal = 1
1036  END DO
1037  END DO
1038  ! S_LC(1,2) = 0.0
1039  ! S_LC(2,1) = S_LC(1,2)
1040  ! E_LC(1,2) = 60.0
1041  ! E_LC(2,1) = E_LC(1,2)
1042 
1043  IF ( liq_crystal == 1 .AND. phas == 1 .AND. stabil == 0 ) THEN
1044 
1045  count_lc = 0
1046  tolerance = 1.e-6
1047 
1048  steps_lc = 200
1049  deltay = 1.0 / REAL(steps_lc)
1050 
1051  ! --- dimensionless function U_LC repres. anisotr. intermolecular interactions in l.c.
1052 
1053  DO i = 1, ncomp
1054  DO j = 1, ncomp
1055  u_lc(i,j) = 2.0/3.0*pi*mseg(i)*mseg(j) *(0.5*(dhs(i)+dhs(j)))**3 & ! sig_ij(i,j)**3
1056  *(e_lc(i,j)/t+s_lc(i,j))*rho
1057  END DO
1058  END DO
1059 
1060 
1061  DO i=1,ncomp
1062  ! S_lc2(i) = 0.0 !for isotropic
1063  s_lc2(i) = 0.5 !for nematic
1064  s_lc1(i) = s_lc2(i)
1065  END DO
1066 
1067 1 CONTINUE
1068 
1069  DO i = 1, ncomp
1070  IF (s_lc2(i) <= 0.3) s_lc1(i) = s_lc2(i)
1071  IF (s_lc2(i) > 0.3) s_lc1(i) = s_lc1(i) + (s_lc2(i)-s_lc1(i))*accel_lc
1072  END DO
1073 
1074  count_lc = count_lc + 1
1075 
1076  !--------------------------------------------------------------------------
1077  ! single-particle orientation partition function Z_LC in liquid crystals
1078  !--------------------------------------------------------------------------
1079 
1080  DO i = 1, ncomp
1081  sumu(i) = 0.0
1082  DO j = 1, ncomp
1083  sumu(i) = sumu(i) + x(j)*u_lc(i,j)*s_lc1(j)
1084  END DO
1085  END DO
1086 
1087  DO i = 1, ncomp
1088  z_lc(i) = 0.0
1089  integrand1 = exp(-0.5*sumu(i)) !eq. for Z_LC with y=0
1090  DO k = 1, steps_lc
1091  integrand2 = exp(0.5*sumu(i)*(3.0*(deltay*REAL(k)) **2 -1.0))
1092  z_lc(i) = z_lc(i) + (integrand1 + integrand2)/2.0*deltay
1093  integrand1 = integrand2
1094  END DO !k-index integration
1095  END DO !i-index Z_LC(i) calculation
1096 
1097  !--------------------------------------------------------------------------
1098  ! order parameter S_lc in liquid crystals
1099  !--------------------------------------------------------------------------
1100 
1101  error_lc = 0.0
1102  DO i = 1, ncomp
1103  s_lc2(i) = 0.0
1104  integrand1 = -1.0/z_lc(i)*0.5*exp(-0.5*sumu(i)) !for S_lc with y=0
1105  DO k = 1, steps_lc
1106  integrand2 = 1.0/z_lc(i)*0.5*(3.0*(deltay*REAL(k)) &
1107  **2 -1.0)*exp(0.5*sumu(i)*(3.0 *(deltay*REAL(k))**2 -1.0))
1108  s_lc2(i) = s_lc2(i) + (integrand1 + integrand2)/2.0*deltay
1109  integrand1 = integrand2
1110  END DO !k-index integration
1111  error_lc = error_lc + abs(s_lc2(i)-s_lc1(i))
1112  END DO !i-index Z_LC(i) calculation
1113 
1114  sphase = 0.0
1115  DO i = 1, ncomp
1116  sphase = sphase + s_lc2(i)
1117  END DO
1118  IF (sphase < 1.e-4) THEN
1119  error_lc = 0.0
1120  DO i = 1, ncomp
1121  s_lc2(i)= 0.0
1122  z_lc(i) = 1.0
1123  END DO
1124  END IF
1125 
1126 
1127  ! write (*,*) count_LC,S_lc2(1)-S_lc1(1),S_lc2(2)-S_lc1(2)
1128  IF (error_lc > tolerance .AND. count_lc < 400) GO TO 1
1129  ! write (*,*) 'done',eta,S_lc2(1),S_lc2(2)
1130 
1131  IF (count_lc == 400) WRITE (*,*) 'LC iteration not converg.'
1132 
1133  !--------------------------------------------------------------------------
1134  ! the anisotropic contribution to the Helmholtz energy
1135  !--------------------------------------------------------------------------
1136 
1137  u_term = 0.0
1138  DO i = 1, ncomp
1139  DO j = 1, ncomp
1140  u_term = u_term + 0.5*x(i)*x(j)*s_lc2(i) *s_lc2(j)*u_lc(i,j)
1141  END DO
1142  END DO
1143 
1144  flc = 0.0
1145  DO i = 1, ncomp
1146  IF (z_lc(i) /= 0.0) flc = flc - x(i) * log(z_lc(i))
1147  END DO
1148  flc = flc + u_term
1149 
1150  END IF
1151  ! write (*,'(i2,i2,4(f15.8))') phas,stabil,flc,eta,S_lc2(1),x(1)
1152 
1153 
1154 END SUBROUTINE f_lc_mayersaupe
1155 
1156 
1157 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1158 !
1159 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1160 
1161 SUBROUTINE f_pert_theory ( fdsp )
1162 
1163  USE eos_variables, ONLY: nc, pi, rho, eta, z0t, order1, order2
1164  USE eos_numerical_derivatives, ONLY: disp_term
1165  USE dft_module
1166 
1167  !-----------------------------------------------------------------------------
1168  REAL, INTENT(IN OUT) :: fdsp
1169  !-----------------------------------------------------------------------------
1170  REAL :: i1, i2
1171  REAL :: z3, zms, c1_con, m_mean
1172  !-----------------------------------------------------------------------------
1173 
1174  ! caution: positive sign of correlation integral is used here !
1175  ! (the Helmholtz energy terms are written with a negative sign, while I1 and I2 are positive)
1176 
1177  IF (disp_term == 'PT1') THEN
1178 
1179  CALL f_dft ( i1, i2)
1180  c1_con = 0.0
1181  i2 = 0.0
1182  fdsp = + ( - 2.0*pi*rho*i1*order1 )
1183 
1184  ELSEIF (disp_term == 'PT2') THEN
1185 
1186  CALL f_dft ( i1, i2)
1187  z3 = eta
1188  zms = 1.0 - z3
1189  m_mean = z0t / ( pi / 6.0 )
1190  c1_con = 1.0/ ( 1.0 + m_mean*(8.0*z3-2.0*z3**2 )/zms**4 &
1191  + (1.0 - m_mean)*( 20.0*z3 -27.0*z3**2 +12.0*z3**3 -2.0*z3**4 ) &
1192  /(zms*(2.0-z3))**2 )
1193  fdsp = + ( - 2.0*pi*rho*i1*order1 - pi*rho*c1_con*m_mean*i2*order2 )
1194 
1195  ELSEIF (disp_term == 'PT_MIX') THEN
1196 
1197  CALL f_pert_theory_mix ( fdsp )
1198 
1199  ELSEIF (disp_term == 'PT_MF') THEN
1200 
1201  ! mean field theory
1202  i1 = - ( - 8.0/9.0 - 4.0/9.0*(rc**-9 -3.0*rc**-3 ) - tau_cut/3.0*(rc**3 -1.0) )
1203  fdsp = + ( - 2.0*pi*rho*i1*order1 )
1204  write (*,*) 'caution: not thoroughly checked and tested'
1205 
1206  ELSE
1207  write (*,*) 'define the type of perturbation theory'
1208  stop
1209  END IF
1210 
1211  ! I1 = I1 + 4.0/9.0*(2.5**-9 -3.0*2.5**-3 )
1212  ! fdsp = + ( - 2.0*PI*rho*I1*order1 - PI*rho*c1_con*m_mean*I2*order2 )
1213 
1214 END SUBROUTINE f_pert_theory
1215 
1216 END MODULE eos_f_contributions
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
eos_numerical_derivatives
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:220