35 SUBROUTINE start_var(converg)
42 INTEGER,
INTENT(IN OUT) :: converg
46 INTEGER :: ncompsav, n_unkwsav, ph_split
47 LOGICAL :: lle_check, flashcase, renormalize
48 REAL :: den1, den2, x_1, x_2
49 CHARACTER (LEN=50) :: filename
55 IF (num == 2) renormalize = .true.
59 IF (xif(1) /= 0.0) flashcase = .true.
65 xi(1,i) = 1.0 /
REAL(ncomp)
66 xi(2,i) = 1.0 /
REAL(ncomp)
73 IF( ncomp == 2 .AND. .NOT.flashcase )
THEN 74 CALL vle_min( lle_check )
82 CALL phase_stability ( .false., flashcase, ph_split )
89 CALL select_sum_rel (1,0,1)
92 CALL select_sum_rel (2,0,2)
98 IF (ph_split == 1)
THEN 106 CALL determine_flash_it2
108 CALL select_sum_rel (1,0,1)
109 CALL select_sum_rel (2,0,2)
114 CALL objective_ctrl (converg)
124 IF (lle_check)
CALL phase_stability (lle_check,flashcase,ph_split)
131 IF (ph_split == 1)
THEN 136 IF (flashcase)
CALL select_sum_rel (1,0,1)
137 IF (flashcase)
CALL select_sum_rel (2,0,2)
143 CALL determine_flash_it2
145 CALL select_sum_rel (1,0,1)
146 CALL select_sum_rel (2,0,2)
152 CALL objective_ctrl (converg)
160 IF (converg == 1)
THEN 164 xi(ph,i) = exp( val_conv(4+i+(ph-1)*ncomp) )
167 dense(1:2) = val_conv(1:2)
174 END SUBROUTINE start_var
213 SUBROUTINE objective_ctrl (converg)
220 INTEGER,
INTENT(OUT) :: converg
224 SUBROUTINE objec_fct ( iter_no, y, residu, dummy )
225 INTEGER,
INTENT(IN) :: iter_no
226 REAL,
INTENT(IN) :: y(iter_no)
227 REAL,
INTENT(OUT) :: residu(iter_no)
228 INTEGER,
INTENT(IN OUT) :: dummy
229 END SUBROUTINE objec_fct
232 INTEGER :: info,k,posn,i
233 INTEGER,
PARAMETER :: mxr = nc*(nc+1)/2
234 REAL,
ALLOCATABLE :: y(:),diag(:),residu(:)
235 REAL :: x_init, x_solut, r_diff1, r_diff2, totres
236 REAL :: r_thrash, x_thrash
237 CHARACTER (LEN=2) :: compon
238 LOGICAL :: convergence
243 ALLOCATE( y(n_unkw), diag(n_unkw), residu(n_unkw) )
245 IF (num == 0) acc_a = 1.e-7
246 IF (num == 0) step_a = 2.e-8
247 IF (num == 1) acc_a = 1.e-7
248 IF (num == 1) step_a = 2.e-8
249 IF (num == 2) acc_a = 5.e-7
250 IF (num == 2) step_a = 1.e-7
255 IF (it(i) ==
't') y(posn) = val_init(3)
256 IF (it(i) ==
'p') y(posn) = val_init(4)
257 IF (it(i) ==
'lnp') y(posn) = log( val_init(4) )
258 IF (it(i) ==
'fls') y(posn) = alpha
259 IF (it(i)(1:1) ==
'x' .AND. it(i)(2:2) ==
'1') compon = it(i)(3:3)
260 IF (it(i)(1:1) ==
'x' .AND. it(i)(2:2) ==
'1')
READ(compon,*) k
261 IF (it(i)(1:1) ==
'x' .AND. it(i)(2:2) ==
'1') y(posn) = val_init(4+k)
262 IF (it(i)(1:1) ==
'x' .AND. it(i)(2:2) ==
'2') compon = it(i)(3:3)
263 IF (it(i)(1:1) ==
'x' .AND. it(i)(2:2) ==
'2')
READ(compon,*) k
264 IF (it(i)(1:1) ==
'x' .AND. it(i)(2:2) ==
'2') y(posn) = val_init(4+ncomp+k)
265 IF (it(i)(1:1) ==
'x' .AND. it(i)(2:2) ==
'3') compon = it(i)(3:3)
266 IF (it(i)(1:1) ==
'x' .AND. it(i)(2:2) ==
'3')
READ(compon,*) k
267 IF (it(i)(1:1) ==
'x' .AND. it(i)(2:2) ==
'3') y(posn) = val_init(4+ncomp+ncomp+k)
274 IF (lnx(1,i) /= 0.0 .AND. lnx(2,i) /= 0.0)
THEN 275 x_init = x_init + abs( 1.0 - lnx(2,i)/lnx(1,i) )
277 x_init = x_init + abs( 1.0 - exp(lnx(2,i))/exp(lnx(1,i)) )
281 CALL hbrd (objec_fct, n_unkw, y, residu, step_a, acc_a, info, diag)
285 IF ( lnx(1,i) /= 0.0 .AND. lnx(2,i) /= 0.0 )
THEN 286 x_solut = x_solut + abs( 1.0 - lnx(2,i)/lnx(1,i) )
288 IF (lnx(1,i) < 1e300 .AND. lnx(1,i) > -1.e300 ) &
289 x_solut = x_solut + abs( 1.0 - exp(lnx(2,i))/exp(lnx(1,i)) )
292 r_diff1 = abs( 1.0 - dense(1)/dense(2) )
293 IF ( val_conv(2) > 0.0 )
THEN 294 r_diff2 = abs( 1.0 - val_conv(1)/val_conv(2) )
299 totres = sum( abs( residu(1:n_unkw) ) )
303 if (num > 0 ) r_thrash = r_thrash * 10.0
304 if (num > 0 ) x_thrash = x_thrash * 100.0
308 IF ( info >= 2 ) convergence = .false.
309 IF ( abs( 1.0- dense(1)/dense(2) ) < r_thrash .AND. x_solut < x_thrash )
THEN 310 IF ( x_init > 0.050 ) convergence = .false.
311 IF ( ( abs( 1.0- dense(1)/dense(2) ) + x_solut ) < 1.e-7 ) convergence = .false.
313 IF ( r_diff2 /= 0.0 .AND. r_diff2 > (4.0*r_diff1) .AND. bindiag == 1 ) convergence = .false.
314 IF ( ncomp == 1 .AND. totres > 100.0*acc_a ) convergence = .false.
315 IF ( totres > 1000.0*acc_a ) convergence = .false.
316 IF ( ncomp == 1 .AND. r_diff1 < 1.d-5 ) convergence = .false.
318 IF ( convergence )
THEN 322 IF (num <= 1)
CALL enthalpy_etc
327 DEALLOCATE( y, diag, residu )
329 END SUBROUTINE objective_ctrl
346 SUBROUTINE objec_fct ( iter_no, y, residu, dummy )
353 INTEGER,
INTENT(IN) :: iter_no
354 REAL,
INTENT(IN) :: y(iter_no)
355 REAL,
INTENT(OUT) :: residu(iter_no)
356 INTEGER,
INTENT(IN OUT) :: dummy
359 INTEGER :: i, ph,k,posn, skip,phase
360 REAL :: lnphi(np,nc),isofugacity
361 CHARACTER (LEN=2) :: compon
368 IF (it(i) ==
't') t = y(posn)
369 IF (it(i) ==
'p') p = y(posn)
370 IF (it(i) ==
'lnp') p = exp( y(posn) )
371 IF (it(i) ==
'fls') alpha = y(posn)
372 IF (it(i)(1:1) ==
'x' .AND. it(i)(2:2) ==
'1') compon = it(i)(3:3)
373 IF (it(i)(1:1) ==
'x' .AND. it(i)(2:2) ==
'1')
READ(compon,*) k
374 IF (it(i)(1:1) ==
'x' .AND. it(i)(2:2) ==
'1') lnx(1,k) = y(posn)
375 IF (it(i)(1:1) ==
'x' .AND. it(i)(2:2) ==
'2') compon = it(i)(3:3)
376 IF (it(i)(1:1) ==
'x' .AND. it(i)(2:2) ==
'2')
READ(compon,*) k
377 IF (it(i)(1:1) ==
'x' .AND. it(i)(2:2) ==
'2') lnx(2,k) = y(posn)
378 IF (it(i)(1:1) ==
'x' .AND. it(i)(2:2) ==
'3') compon = it(i)(3:3)
379 IF (it(i)(1:1) ==
'x' .AND. it(i)(2:2) ==
'3')
READ(compon,*) k
380 IF (it(i)(1:1) ==
'x' .AND. it(i)(2:2) ==
'3') lnx(3,k) = y(posn)
384 IF (lnx(1,k) > 0.0) lnx(1,k) = 0.0
385 IF (lnx(2,k) > 0.0) lnx(2,k) = 0.0
388 IF (p < 1.e-100) p = 1.e-12
392 IF ( p /= p ) p = 1000.0
393 IF ( t /= t ) t = 300.0
394 IF ( alpha /= alpha ) alpha = 0.5
399 IF ( lnx(ph,i) < -300.0 )
THEN 402 xi(ph,i) = exp( lnx(ph,i) )
407 IF (ncomp > 1)
CALL x_summation
409 CALL fugacity (lnphi)
414 IF (n_unkw < (ncomp*(nphas-1))) skip = ncomp*(nphas-1) - n_unkw
415 IF ((i+skip-ncomp*(phase-2)) > ncomp) phase = phase + 1
416 residu(i) = isofugacity((i+skip-ncomp*(phase-2)),phase,lnphi)
417 if ( density_error(phase) /= 0.0 ) residu(i) = residu(i) + sign( density_error(phase), residu(i) ) * 0.001
420 END SUBROUTINE objec_fct
431 REAL FUNCTION isofugacity (i,phase,lnphi)
437 INTEGER,
INTENT(IN) :: i
438 INTEGER,
INTENT(IN) :: phase
439 REAL,
INTENT(IN) :: lnphi(np,nc)
450 isofugacity = scaling(i) *( lnphi(p2,i)+lnx(p2,i)-lnx(p1,i)-lnphi(p1,i) )
458 END FUNCTION isofugacity
480 SUBROUTINE vle_min(lle_check)
488 LOGICAL,
INTENT(OUT) :: lle_check
490 INTEGER :: i,j,k,phasen(0:40),steps
492 REAL :: vlemin(0:40),llemin(0:40),xval(0:40)
493 REAL :: start_xv(0:40),start_xl(0:40),x_sav,dg_dx2
515 xi(1,1) = 1.0 -
REAL(i) /
REAL(steps)
516 IF ( xi(1,1) <= 1.e-50 ) xi(1,1) = 1.e-50
518 lnx(1,1) = log(xi(1,1))
519 lnx(2,1) = log(xi(2,1))
522 CALL fugacity (lnphi)
529 llemin(i)= gibbs(1) +(xi(1,1)*lnx(1,1)+xi(1,2)*lnx(1,2))*rgas*t
531 IF ( abs(1.0-dense(1)/dense(2)) > 0.0001 )
THEN 532 vlemin(i)= gibbs(1) +(xi(1,1)*lnx(1,1)+xi(1,2)*lnx(1,2))*rgas*t &
533 - ( gibbs(2) +(xi(2,1)*lnx(2,1)+xi(2,2)*lnx(2,2))*rgas*t )
539 IF (i > 0 .AND. phasen(i) == 2)
THEN 540 IF (phasen(i-1) == 2 .AND. abs(vlemin(i)+vlemin(i-1)) < &
541 abs(vlemin(i))+abs(vlemin(i-1)))
THEN 543 start_xv(j)=xval(i-1) + (xval(i)-xval(i-1)) &
544 * abs(vlemin(i-1))/abs(vlemin(i)-vlemin(i-1))
552 dg_dx2 = (-llemin(i-2)+16.0*llemin(i-1)-30.0*llemin(i) &
553 +16.0*llemin(i+1)-llemin(i+2)) / (12.0*((xval(i)-xval(i-1))**2))
554 IF (dg_dx2 < 0.0)
THEN 560 IF (start_xl(1) == 0.0 .AND. start_xv(1) /= 0.0)
THEN 562 xi(1,1) = start_xv(1)
563 xi(1,2) = 1.0-xi(1,1)
566 ELSE IF (start_xl(1) /= 0.0 .AND. start_xv(1) == 0.0)
THEN 568 xi(1,1) = start_xl(1)
569 xi(1,2) = 1.0-xi(1,1)
572 ELSE IF (start_xl(1) /= 0.0 .AND. start_xv(1) /= 0.0)
THEN 574 xi(1,1) = start_xv(1)
575 xi(1,2) = 1.0-xi(1,1)
581 xi(1,2) = 1.0 - xi(1,1)
587 END SUBROUTINE vle_min
599 SUBROUTINE phase_stability ( lle_check, flashcase, ph_split )
603 USE eos_variables, ONLY: dhs, pi, x, eta, eta_start, z3t, fres
609 LOGICAL,
INTENT(IN OUT) :: flashcase
610 INTEGER,
INTENT(OUT) :: ph_split
614 REAL FUNCTION f_stability ( optpara, n )
615 INTEGER,
INTENT(IN) :: n
616 REAL,
INTENT(IN OUT) :: optpara(n)
643 REAL :: fmin, t0, h0, machep, praxis
644 REAL,
ALLOCATABLE :: optpara(:)
646 INTEGER :: i, feedphases, trial
647 REAL :: rhoi(nc),rho_start
648 REAL :: feeddens, rho_phas(np)
653 REAL :: w(np,nc), mean_mass
657 ALLOCATE( optpara(n) )
662 IF (.NOT.flashcase) xif(i) = xi(1,i)
663 IF (flashcase) xi(1,i) = xif(i)
671 CALL dens_calc(rho_phas)
672 IF ( abs(1.0-dense(1)/dense(2)) > 0.0005 )
THEN 691 w(2,i) = 1.0 /
REAL(ncomp)
693 mean_mass = 1.0 / sum( w(2,1:ncomp)/mm(1:ncomp) )
694 xi(2,1:ncomp) = w(2,1:ncomp)/mm(1:ncomp) * mean_mass
699 rhoif(i) = rho_phas(1) * xif(i)
710 rhot = sum( rhoi(1:ncomp) )
711 x(1:ncomp) = rhoi(1:ncomp) / rhot
712 CALL perturbation_parameter
713 xi(1,1:ncomp) = x(1:ncomp)
716 densta(1) = eta_start
718 CALL fugacity (lnphi)
721 call fden_calc ( fden, rhoi )
724 grad_fd(1:ncomp) = lnphi(1,1:ncomp) + log( rhoi(1:ncomp) )
733 CALL dens_calc(rho_phas)
734 rho_start = rho_phas(2)*0.45/dense(2)
737 rhoi(i) = xi(2,i)*rho_start
738 optpara(i) = log( rhoi(i) )
749 fmin = praxis( t0, machep, h0, n, prin, optpara, f_stability, fmin )
757 fmin = f_stability( optpara, n )
774 IF (fmin < -1.e-7 .AND. &
775 abs( 1.0 - maxval(exp(optpara),mask=optpara /= 0.0) /maxval(rhoif) ) > 0.0005)
THEN 779 IF (ph_split == 1)
THEN 786 rhoi2(1:ncomp) = exp( optpara(1:ncomp) )
787 dens = pi/6.0 * sum( rhoi2(1:ncomp) * parame(1:ncomp,1) * dhs(1:ncomp)**3 )
788 rhot = sum( rhoi2(1:ncomp) )
789 xi(2,1:ncomp) = rhoi2(1:ncomp) / rhot
793 IF (trial <= ncomp + ncomp)
THEN 797 IF (trial <= ncomp)
THEN 799 w(2,i) = 1.0 /
REAL(ncomp-1) * 0.05
804 w(2,i) = 1.0 /
REAL(ncomp-1) * 0.00001
806 w(2,trial-ncomp) = 0.99999
808 mean_mass = 1.0 / sum( w(2,1:ncomp)/mm(1:ncomp) )
809 xi(2,1:ncomp) = w(2,1:ncomp)/mm(1:ncomp) * mean_mass
815 IF (feedphases > 1 .AND. .NOT.lle_check .AND. densta(1) > 0.2)
THEN 817 CALL dens_calc(rho_phas)
824 DEALLOCATE( optpara )
826 END SUBROUTINE phase_stability
851 SUBROUTINE select_sum_rel (ph,excl,startindex)
857 INTEGER,
INTENT(IN) :: ph
858 INTEGER,
INTENT(IN) :: excl
859 INTEGER,
INTENT(IN) :: startindex
861 INTEGER :: i,j, sum_index
869 IF ( xi(ph,i) > xmax(ph) )
THEN 873 IF (ph == 1 .AND. i == 1) sum_rel(1) =
'x11' 874 IF (ph == 1 .AND. i == 2) sum_rel(1) =
'x12' 875 IF (ph == 1 .AND. i == 3) sum_rel(1) =
'x13' 876 IF (ph == 1 .AND. i == 4) sum_rel(1) =
'x14' 877 IF (ph == 1 .AND. i == 5) sum_rel(1) =
'x15' 879 IF (ph == 2 .AND. i == 1) sum_rel(2) =
'x21' 880 IF (ph == 2 .AND. i == 2) sum_rel(2) =
'x22' 881 IF (ph == 2 .AND. i == 3) sum_rel(2) =
'x23' 882 IF (ph == 2 .AND. i == 4) sum_rel(2) =
'x24' 883 IF (ph == 2 .AND. i == 5) sum_rel(2) =
'x25' 885 IF (ph == 3 .AND. i == 1) sum_rel(3) =
'x31' 886 IF (ph == 3 .AND. i == 2) sum_rel(3) =
'x32' 887 IF (ph == 3 .AND. i == 3) sum_rel(3) =
'x33' 888 IF (ph == 3 .AND. i == 4) sum_rel(3) =
'x34' 889 IF (ph == 3 .AND. i == 5) sum_rel(3) =
'x35' 898 IF ( i /= sum_index .AND. i /= excl )
THEN 899 IF (ph == 1 .AND. i == 1) it(startindex+j) =
'x11' 900 IF (ph == 1 .AND. i == 2) it(startindex+j) =
'x12' 901 IF (ph == 1 .AND. i == 3) it(startindex+j) =
'x13' 902 IF (ph == 1 .AND. i == 4) it(startindex+j) =
'x14' 903 IF (ph == 1 .AND. i == 5) it(startindex+j) =
'x15' 905 IF (ph == 2 .AND. i == 1) it(startindex+j) =
'x21' 906 IF (ph == 2 .AND. i == 2) it(startindex+j) =
'x22' 907 IF (ph == 2 .AND. i == 3) it(startindex+j) =
'x23' 908 IF (ph == 2 .AND. i == 4) it(startindex+j) =
'x24' 909 IF (ph == 2 .AND. i == 5) it(startindex+j) =
'x25' 911 IF (ph == 3 .AND. i == 1) it(startindex+j) =
'x31' 912 IF (ph == 3 .AND. i == 2) it(startindex+j) =
'x32' 913 IF (ph == 3 .AND. i == 3) it(startindex+j) =
'x33' 914 IF (ph == 3 .AND. i == 4) it(startindex+j) =
'x34' 915 IF (ph == 3 .AND. i == 5) it(startindex+j) =
'x35' 922 END SUBROUTINE select_sum_rel
931 SUBROUTINE tangent_plane
955 REAL FUNCTION praxis( t0, MACHEP, h0, n, PRIN, optpara, TANGENT_VALUE, fmin )
956 REAL,
INTENT(IN OUT) :: t0
957 REAL,
INTENT(IN) :: machep
958 REAL,
INTENT(IN) :: h0
960 INTEGER,
INTENT(IN OUT) :: prin
961 REAL,
INTENT(IN OUT) :: optpara(n)
962 REAL,
EXTERNAL :: tangent_value
963 REAL,
INTENT(IN OUT) :: fmin
966 REAL FUNCTION tangent_value2 ( optpara, n )
967 INTEGER,
INTENT(IN) :: n
968 REAL,
INTENT(IN) :: optpara(n)
975 INTEGER :: small_i, min_ph, other_ph
977 REAL :: fmin , t0, h0, machep
979 REAL,
ALLOCATABLE :: optpara(:)
992 ALLOCATE( optpara(n) )
1000 lnx(1,i) = log(xi(1,i))
1001 lnx(2,i) = log(xi(2,i))
1005 optpara(i) = log( xi(2,i) * 0.001 )
1017 fmin = praxis( t0, machep, h0, n, prin, optpara, tangent_value2, fmin )
1021 fmin = tangent_value2( optpara, n )
1030 IF ( minval( lnx(1,1:ncomp) ) < minval( lnx(2,1:ncomp) ) )
THEN 1037 small_i = minloc( lnx(min_ph,1:ncomp), 1 )
1039 IF ( minval( lnx(min_ph,1:ncomp) ) < -20.0 )
THEN 1040 CALL fugacity ( lnphi )
1041 lnx(min_ph,small_i) = lnx(other_ph,small_i)+lnphi(other_ph,small_i) - lnphi(min_ph,small_i)
1042 optpara(small_i) = lnx(2,small_i) + log( sum( exp( optpara(1:ncomp) ) ) )
1050 val_init(1) = dense(1)
1051 val_init(2) = dense(2)
1056 val_init(4+k+(ph-1)*ncomp) = lnx(ph,k)
1063 DEALLOCATE( optpara )
1065 END SUBROUTINE tangent_plane
1072 SUBROUTINE determine_flash_it2
1079 REAL :: n_phase1, n_phase2, max_x_diff
1082 IF ( minval( lnx(1,1:ncomp) ) < minval( lnx(2,1:ncomp) ) )
THEN 1085 IF (ncomp >= 3) it(3) =
'x13' 1086 IF (ncomp >= 4) it(4) =
'x14' 1087 IF (ncomp >= 5) it(5) =
'x15' 1092 IF (ncomp >= 3) it(3) =
'x23' 1093 IF (ncomp >= 4) it(4) =
'x24' 1094 IF (ncomp >= 5) it(5) =
'x25' 1099 IF ( abs( exp( lnx(1,i) ) - exp( lnx(2,i) ) ) > max_x_diff )
THEN 1100 max_x_diff = abs( exp( lnx(1,i) ) - exp( lnx(2,i) ) )
1101 n_phase1 = ( xif(i) - exp( lnx(2,i) ) ) / ( exp( lnx(1,i) ) - exp( lnx(2,i) ) )
1102 n_phase2 = 1.0 - n_phase1
1105 lnx(1,1:ncomp) = lnx(1,1:ncomp) + log( n_phase1 )
1106 lnx(2,1:ncomp) = lnx(2,1:ncomp) + log( n_phase2 )
1109 val_init(1) = dense(1)
1110 val_init(2) = dense(2)
1115 val_init(4+k+(ph-1)*ncomp) = lnx(ph,k)
1120 END SUBROUTINE determine_flash_it2
1139 SUBROUTINE poly_sta_var (converg)
1145 INTEGER,
INTENT(IN OUT) :: converg
1149 REAL :: p_spec,solution(10,4+nc*np)
1154 find_equilibrium:
DO 1156 CALL polymer_free(p_spec,sol,solution)
1159 WRITE (*,*)
' GENERATING STARTING VALUES' 1161 val_init(1) = solution(1,1)
1162 val_init(2) = solution(1,2)
1163 val_init(3) = solution(1,3)
1164 val_init(4) = solution(1,4)
1167 val_init(4+k+(ph-1)*ncomp) = solution(1,4+k+(ph-1)*ncomp)
1170 val_init(7) = -10000.0
1173 WRITE (*,*)
' INITIAL EQUILIBRIUM CALC. FAILD. NEXT STEP STARTS' 1175 IF (p == p_spec)
THEN 1179 CALL objective_ctrl (converg)
1183 CALL phase_equilib(p_spec,5.0,converg)
1186 IF (converg == 1)
EXIT find_equilibrium
1188 IF ( p < (0.7*p_spec) )
WRITE (*,*)
' NO SOLUTION FOUND' 1189 IF ( p < (0.7*p_spec) ) stop
1191 END DO find_equilibrium
1193 WRITE (*,*)
' FINISHED: POLY_STA_VAR' 1195 END SUBROUTINE poly_sta_var
1213 SUBROUTINE x_summation
1219 INTEGER :: i, j, comp_i, ph_i
1221 CHARACTER (LEN=2) :: phasno
1222 CHARACTER (LEN=2) :: compno
1223 LOGICAL :: flashcase2
1227 IF (sum_rel(j)(1:3) ==
'nfl')
THEN 1235 flashcase2 = .false.
1239 IF (sum_rel(j)(1:1) ==
'x')
THEN 1241 phasno = sum_rel(j)(2:2)
1243 compno = sum_rel(j)(3:3)
1244 READ(compno,*) comp_i
1245 IF ( sum_rel(nphas+j)(1:1) ==
'e' )
CALL neutr_charge(nphas+j)
1249 IF ( i /= comp_i ) sum_x = sum_x + xi(ph_i,i)
1251 xi(ph_i,comp_i) = 1.0 - sum_x
1252 IF ( xi(ph_i,comp_i ) < 0.0 ) xi(ph_i,comp_i) = 0.0
1253 IF ( xi(ph_i,comp_i ) /= 0.0 )
THEN 1254 lnx(ph_i,comp_i) = log( xi(ph_i,comp_i) )
1256 lnx(ph_i,comp_i) = -100000.0
1260 ELSE IF ( sum_rel(j)(1:2) ==
'fl' )
THEN 1276 WRITE (*,*)
'summation relation not defined' 1282 IF ( it(1) ==
'fls' )
CALL flash_sum
1283 IF ( flashcase2 )
CALL flash_alpha
1285 END SUBROUTINE x_summation
1305 SUBROUTINE fugacity (ln_phi)
1308 USE eos_variables, ONLY: phas, x, eta, eta_start, lnphi, fres, rho, pges, kbol
1312 REAL,
INTENT(OUT) :: ln_phi(np,nc)
1323 eta_start = densta(ph)
1324 x(1:ncomp) = xi(ph,1:ncomp)
1326 IF (p < 1.e-100)
THEN 1327 WRITE(*,*)
' FUGACITY: PRESSURE TOO LOW', p
1331 IF (num == 0)
CALL phi_eos
1332 IF (num == 1)
CALL phi_numerical
1334 IF (num == 2)
write(*,*)
'CRITICAL RENORM NOT INCLUDED YET' 1337 ln_phi(ph,1:ncomp) = lnphi(1:ncomp)
1363 END SUBROUTINE fugacity
1377 SUBROUTINE enthalpy_etc
1393 x(1:ncomp) = xi(ph,1:ncomp)
1398 IF(num == 1)
CALL h_numerical
1399 IF(num == 2)
write (*,*)
'enthalpy_etc: incorporate H_EOS_RN' 1409 IF (nphas == 2) h_lv = enthal(2)-enthal(1)
1413 END SUBROUTINE enthalpy_etc
1424 SUBROUTINE dens_calc(rho_phas)
1432 REAL,
INTENT(OUT) :: rho_phas(np)
1444 eta_start = densta(ph)
1445 x(1:ncomp) = xi(ph,1:ncomp)
1447 CALL perturbation_parameter
1448 CALL density_iteration
1451 rho_phas(ph) = eta/z3t
1454 write (*,*)
' SUBROUTINE DENS_CALC not available for cubic EOS' 1460 END SUBROUTINE dens_calc
1467 SUBROUTINE fden_calc (fden, rhoi)
1474 REAL,
INTENT(OUT) :: fden
1475 REAL,
INTENT(IN OUT) :: rhoi(nc)
1477 REAL :: rhot, fden_id
1483 rhot = sum( rhoi(1:ncomp) )
1484 x(1:ncomp) = rhoi(1:ncomp) / rhot
1486 CALL perturbation_parameter
1492 ELSE IF(num == 1)
THEN 1495 write (*,*)
'deactivated this line when making a transition to f90' 1500 fden_id = sum( rhoi(1:ncomp) * ( log( rhoi(1:ncomp) ) - 1.0 ) )
1502 fden = fres * rhot + fden_id
1505 write (*,*)
' SUBROUTINE FDEN_CALC not available for cubic EOS' 1509 END SUBROUTINE fden_calc
1527 SUBROUTINE polymer_free (p_spec,sol,solution)
1533 REAL,
INTENT(IN OUT) :: p_spec
1534 INTEGER,
INTENT(OUT) :: sol
1535 REAL,
INTENT(OUT) :: solution(10,4+nc*np)
1538 INTEGER :: k,j,ph, converg
1553 DO WHILE ( sol == 0 )
1559 xi(1,1) = grid(j) / ( (1.0-grid(j)) * mm(1) / mm(2) )
1560 IF ( mm(1) < 5000.0 ) xi(1,1) = xi(1,1) * 0.8
1561 xi(1,2) = 1.0 - xi(1,1)
1562 lnx(1,1) = log(xi(1,1))
1563 lnx(1,2) = log(xi(1,2))
1570 val_init(2) = 0.0001
1576 val_init(4+k+(ph-1)*ncomp) = lnx(ph,k)
1589 CALL objective_ctrl (converg)
1591 IF (converg == 1 .AND. abs(dense(1)/dense(2)-1.0) > 1.d-3 .AND. dense(1) > 0.1)
THEN 1594 DO k = 1,4+ncomp*nphas
1595 solution(sol,k) = val_conv(k)
1597 ELSE IF (abs(solution(sol,5)/lnx(1,1)-1.0) > 1.d-2)
THEN 1599 DO k = 1,4+ncomp*nphas
1600 solution(sol,k) = val_conv(k)
1612 WRITE (*,*)
' no initial solution found' 1614 IF (p < (0.7*p_spec))
WRITE (*,*)
' NO SOLUTION FOUND' 1615 IF (p < (0.7*p_spec)) stop
1616 ELSE IF (sol > 1)
THEN 1633 END SUBROUTINE polymer_free
1670 SUBROUTINE phase_equilib (end_x,steps,converg)
1676 REAL,
INTENT(IN) :: end_x
1677 REAL,
INTENT(IN) :: steps
1678 INTEGER,
INTENT(OUT) :: converg
1681 INTEGER :: k, count1,count2,runindex,maxiter
1682 REAL :: delta_x,delta_org,val_org,runvar
1683 CHARACTER (LEN=2) :: compon
1684 LOGICAL :: continue_cycle
1687 IF (running(1:2) ==
'd1') runindex = 1
1688 IF (running(1:2) ==
'd2') runindex = 2
1689 IF (running(1:1) ==
't') runindex = 3
1690 IF (running(1:1) ==
'p') runindex = 4
1691 IF (running(1:2) ==
'x1') compon = running(3:3)
1692 IF (running(1:2) ==
'x1')
READ(compon,*) k
1693 IF (running(1:2) ==
'x1') runindex = 4+k
1694 IF (running(1:2) ==
'x2') compon = running(3:3)
1695 IF (running(1:2) ==
'x2')
READ(compon,*) k
1696 IF (running(1:2) ==
'x2') runindex = 4+ncomp+k
1697 IF (running(1:2) ==
'l1') compon = running(3:3)
1698 IF (running(1:2) ==
'l1')
READ(compon,*) k
1699 IF (running(1:2) ==
'l1') runindex = 4+k
1700 IF (running(1:2) ==
'l2') compon = running(3:3)
1701 IF (running(1:2) ==
'l2')
READ(compon,*) k
1702 IF (running(1:2) ==
'l2') runindex = 4+ncomp+k
1705 IF ( ncomp >= 3 ) maxiter = 1000
1708 delta_x = ( end_x - val_init(runindex) ) / steps
1709 delta_org = ( end_x - val_init(runindex) ) / steps
1710 val_org = val_init(runindex)
1711 IF ( running(1:1) ==
'x' )
THEN 1712 delta_x = ( end_x - exp(val_init(runindex)) ) / steps
1713 delta_org = ( end_x - exp(val_init(runindex)) ) / steps
1714 val_org = exp(val_init(runindex))
1717 continue_cycle = .true.
1719 DO WHILE ( continue_cycle )
1726 CALL objective_ctrl (converg)
1728 IF (converg == 1)
THEN 1729 val_init( 1:(4+ncomp*nphas) ) = val_conv( 1:(4+ncomp*nphas) )
1730 IF (outp == 1 .AND. (abs(delta_x) > 0.1*abs(delta_org) .OR. count2 == 2))
CALL output
1732 delta_x = delta_x / 2.0
1733 IF (num == 2) delta_x = delta_x / 2.0
1734 val_init(runindex) = val_org
1735 IF (running(1:1) ==
'x') val_init(runindex) = log(val_org)
1736 continue_cycle = .true.
1739 runvar = val_init(runindex)
1740 IF (running(1:1) ==
'x') runvar = exp(val_init(runindex))
1742 IF ( end_x == 0.0 .AND. running(1:1) /=
'x' )
THEN 1743 IF ( abs(runvar-end_x) < 1.e-8 ) continue_cycle = .false.
1744 ELSE IF ( abs((runvar-end_x)/end_x) < 1.e-8 )
THEN 1746 continue_cycle = .false.
1747 ELSE IF ( count1 == maxiter )
THEN 1748 WRITE (*,*)
' MAX. NO OF ITERATIONS',count1
1750 continue_cycle = .false.
1751 ELSE IF ( abs(delta_x) < 1.e-5*abs(delta_org) )
THEN 1754 continue_cycle = .false.
1756 continue_cycle = .true.
1758 IF (abs(runvar+delta_x-end_x) > abs(runvar-end_x)) delta_x = end_x - runvar
1759 val_init(runindex) = runvar + delta_x
1760 IF (running(1:1) ==
'x') val_init(runindex) = log(runvar+delta_x)
1763 IF (abs(delta_x) < abs(delta_org) .AND. count2 >= 5)
THEN 1764 delta_x = delta_x * 2.0
1770 END SUBROUTINE phase_equilib
1776 SUBROUTINE new_flash (ph_it)
1782 INTEGER,
INTENT(IN) :: ph_it
1785 INTEGER :: i, ph_cal
1786 REAL,
DIMENSION(nc) :: ni_1, ni_2
1792 IF ( lnx(ph_it,i) < -300.0 )
THEN 1795 ni_2(i) = exp( lnx(ph_it,i) )
1800 ni_1(i) = xif(i)-ni_2(i)
1801 IF ( ni_2(i) > xif(i) )
THEN 1803 ni_1(i) = xif(i) * 1.e-20
1807 xi(ph_it,1:ncomp) = ni_2(1:ncomp) / sum( ni_2(1:ncomp) )
1809 IF ( xi(ph_it,i) >= 1.e-300 ) lnx(ph_it,i) = log( xi(ph_it,i) )
1811 xi(ph_cal,1:ncomp) = ni_1(1:ncomp) / sum( ni_1(1:ncomp) )
1812 lnx(ph_cal,1:ncomp) = log( xi(ph_cal,1:ncomp) )
1814 END SUBROUTINE new_flash
1840 INTEGER :: i, j, k, ki, l, m
1841 REAL :: z0, z1, z2, z3, z0_rk, z1_rk, z2_rk, z3_rk
1843 REAL,
DIMENSION(nc) :: mhs, mdsp, mhc, myres
1844 REAL,
DIMENSION(nc) :: m_rk
1845 REAL :: gij_rk(nc,nc)
1849 REAL :: i1, i2, i1_rk, i2_rk
1850 REAL :: ord1_rk, ord2_rk
1851 REAL :: c1_con, c2_con, c1_rk
1852 REAL :: zmr, nmr, zmr_rk, nmr_rk, um_rk
1853 REAL,
DIMENSION(nc,0:6) :: ap_rk, bp_rk
1856 REAL :: ass_s2, m_hbon(nc)
1858 REAL :: fdd_rk, fqq_rk, fdq_rk
1859 REAL,
DIMENSION(nc) :: my_dd, my_qq, my_dq
1865 CALL perturbation_parameter
1871 IF (ensemble_flag ==
'tp')
THEN 1872 CALL density_iteration
1873 ELSEIF (ensemble_flag ==
'tv')
THEN 1877 write (*,*)
'PHI_EOS: define ensemble, ensemble_flag == (pv) or (pt)' 1889 m_mean = z0t / (pi/6.0)
1895 zges = (p * 1.d-30) / (kbol*t*rho)
1896 IF ( ensemble_flag ==
'tv' ) zges = (pges * 1.d-30) / (kbol*t*rho)
1907 z0_rk = pi/6.0 * mseg(k)
1908 z1_rk = pi/6.0 * mseg(k) * dhs(k)
1909 z2_rk = pi/6.0 * mseg(k) * dhs(k)*dhs(k)
1910 z3_rk = pi/6.0 * mseg(k) * dhs(k)**3
1913 m_rk(k) = ( mseg(k) - m_mean ) / rho
1928 mhs(k) = 6.0/pi* ( 3.0*(z1_rk*z2+z1*z2_rk)/zms + 3.0*z1*z2*z3_rk/zms/zms &
1929 + 3.0*z2*z2*z2_rk/z3/zms/zms + z2**3 *z3_rk*(3.0*z3-1.0)/z3/z3/zms**3 &
1930 + ((3.0*z2*z2*z2_rk*z3-2.0*z2**3 *z3_rk)/z3**3 -z0_rk) *log(zms) &
1931 + (z0-z2**3 /z3/z3)*z3_rk/zms )
1938 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
1939 gij_rk(i,j) = z3_rk/zms/zms &
1940 + 3.0*dij_ab(i,j)*(z2_rk+2.0*z2*z3_rk/zms)/zms/zms &
1941 + dij_ab(i,j)**2 *z2/zms**3 *(4.0*z2_rk+6.0*z2*z3_rk/zms)
1947 mhc(k) = mhc(k) + x(i)*rho * (1.0-mseg(i)) / gij(i,i) * gij_rk(i,i)
1949 mhc(k) = mhc(k) + ( 1.0-mseg(k)) * log( gij(k,k) )
1959 ap_rk(k,m) = m_rk(k)/m_mean**2 * ( ap(m,2) + (3.0 -4.0/m_mean) *ap(m,3) )
1960 bp_rk(k,m) = m_rk(k)/m_mean**2 * ( bp(m,2) + (3.0 -4.0/m_mean) *bp(m,3) )
1968 i1 = i1 + apar(m)*eta**
REAL(m)
1969 i2 = i2 + bpar(m)*eta**
REAL(m)
1970 i1_rk = i1_rk + apar(m)*
REAL(m)*eta**
REAL(m-1)*z3_rk + ap_rk(k,m)*eta**
REAL(m)
1971 i2_rk = i2_rk + bpar(m)*
REAL(m)*eta**
REAL(m-1)*z3_rk + bp_rk(k,m)*eta**
REAL(m)
1977 ord1_rk = ord1_rk + 2.0*mseg(k)*rho*x(i)*mseg(i)*sig_ij(i,k)**3 *uij(i,k)/t
1978 ord2_rk = ord2_rk + 2.0*mseg(k)*rho*x(i)*mseg(i)*sig_ij(i,k)**3 *(uij(i,k)/t)**2
1981 c1_con= 1.0/ ( 1.0 + m_mean*(8.0*z3-2.0*z3*z3)/zms**4 &
1982 + (1.0 - m_mean)*(20.0*z3-27.0*z3*z3 +12.0*z3**3 -2.0*z3**4 ) &
1983 /(zms*(2.0-z3))**2 )
1984 c2_con= - c1_con*c1_con *( m_mean*(-4.0*z3*z3+20.0*z3+8.0)/zms**5 &
1985 + (1.0 - m_mean) *(2.0*z3**3 +12.0*z3*z3-48.0*z3+40.0) &
1986 /(zms*(2.0-z3))**3 )
1987 c1_rk= c2_con*z3_rk - c1_con*c1_con*m_rk(k) * ( (8.0*z3-2.0*z3*z3)/zms**4 &
1988 - (-2.0*z3**4 +12.0*z3**3 -27.0*z3*z3+20.0*z3) / (zms*(2.0-z3))**2 )
1990 mdsp(k) = -2.0*pi* ( order1*rho*rho*i1_rk + ord1_rk*i1 ) &
1991 - pi* c1_con*m_mean * ( order2*rho*rho*i2_rk + ord2_rk*i2 ) &
1992 - pi* ( c1_con*m_rk(k) + c1_rk*m_mean ) * order2*rho*rho*i2
2005 zmr = zmr + x(i)*x(j)*mseg(i)*mseg(j)*vij(i,j)*uij(i,j)
2006 nmr = nmr + x(i)*x(j)*mseg(i)*mseg(j)*vij(i,j)
2008 zmr_rk = zmr_rk + 2.0*mseg(k) * x(i)*mseg(i)*vij(k,i)*uij(k,i)
2009 nmr_rk = nmr_rk + 2.0*mseg(k) * x(i)*mseg(i)*vij(k,i)
2012 um_rk = 1.0/nmr**2 * ( nmr*zmr_rk - zmr*nmr_rk )
2017 mdsp(k) = mdsp(k) + dnm(i,j)*(um/t)**
REAL(i)*(eta/tau)**
REAL(j) &
2018 * ( 1.0 + z3_rk*rho/eta*
REAL(j) + um_rk*rho/um*
REAL(i) )
2032 IF (nhb_typ(i) /= 0) assoc = .true.
2037 DO l = 1, nhb_typ(k)
2038 ass_s2 = ass_s2 + nhb_no(k,l) * log(mx(k,l))
2043 DO ki = 1, nhb_typ(i)
2045 DO l = 1, nhb_typ(j)
2046 m_hbon(k)= m_hbon(k) - rho*rho/2.0*x(i)*x(j) *mx(i,ki)*mx(j,l) *nhb_no(i,ki)*nhb_no(j,l) &
2047 * gij_rk(i,j) * ass_d(i,j,ki,l)
2060 CALL phi_polar ( k, z3_rk, fdd_rk, fqq_rk, fdq_rk )
2069 myres(k) = mhs(k) +mhc(k) +mdsp(k) +m_hbon(k) +my_dd(k) +my_qq(k) +my_dq(k)
2082 IF (ensemble_flag ==
'tp' ) lnphi(k) = myres(k) - log(zges)
2083 IF (ensemble_flag ==
'tv' ) lnphi(k) = myres(k)
2091 END SUBROUTINE phi_eos
2098 SUBROUTINE phi_numerical
2108 REAL :: fres1, fres2, fres3, fres4, fres5
2110 REAL,
DIMENSION(nc) :: myres
2111 REAL,
DIMENSION(nc) :: rhoi, rhoi_0
2112 REAL :: tfr_1, tfr_2, tfr_3, tfr_4, tfr_5
2120 IF (ensemble_flag ==
'tp')
THEN 2121 CALL density_iteration
2122 ELSEIF (ensemble_flag ==
'tv')
THEN 2126 write (*,*)
'PHI_EOS: define ensemble, ensemble_flag == (tv) or (tp)' 2134 zges = (p * 1.e-30) / (kbol*t*eta/z3t)
2135 IF ( ensemble_flag ==
'tv' ) zges = (pges * 1.e-30) / (kbol*t*eta/z3t)
2139 rhoi_0(1:ncomp) = x(1:ncomp) * eta/z3t
2140 rhoi(1:ncomp) = rhoi_0(1:ncomp)
2149 IF ( rhoi_0(k) > 1.d-9 )
THEN 2150 delta_rho = 1.e-13 * 10.0**(0.5*(15.0+log10(rhoi_0(k))))
2155 rhoi(k) = rhoi_0(k) + delta_rho
2156 eta = pi/6.0 * sum( rhoi(1:ncomp)*mseg(1:ncomp)*dhs(1:ncomp)**3 )
2157 x(1:ncomp) = rhoi(1:ncomp) / sum( rhoi(1:ncomp) )
2158 rho = sum( rhoi(1:ncomp) )
2163 rhoi(k) = rhoi_0(k) + 0.5 * delta_rho
2164 eta = pi/6.0 * sum( rhoi(1:ncomp)*mseg(1:ncomp)*dhs(1:ncomp)**3 )
2165 x(1:ncomp) = rhoi(1:ncomp) / sum( rhoi(1:ncomp) )
2166 rho = sum( rhoi(1:ncomp) )
2171 IF ( rhoi_0(k) > 1.e-9 )
THEN 2172 rhoi(k) = rhoi_0(k) - 0.5 * delta_rho
2173 eta = pi/6.0 * sum( rhoi(1:ncomp)*mseg(1:ncomp)*dhs(1:ncomp)**3 )
2174 x(1:ncomp) = rhoi(1:ncomp) / sum( rhoi(1:ncomp) )
2175 rho = sum( rhoi(1:ncomp) )
2180 rhoi(k) = rhoi_0(k) - delta_rho
2181 eta = pi/6.0 * sum( rhoi(1:ncomp)*mseg(1:ncomp)*dhs(1:ncomp)**3 )
2182 x(1:ncomp) = rhoi(1:ncomp) / sum( rhoi(1:ncomp) )
2183 rho = sum( rhoi(1:ncomp) )
2190 eta = pi/6.0 * sum( rhoi(1:ncomp)*mseg(1:ncomp)*dhs(1:ncomp)**3 )
2191 x(1:ncomp) = rhoi(1:ncomp) / sum( rhoi(1:ncomp) )
2192 rho = sum( rhoi(1:ncomp) )
2197 IF ( rhoi_0(k) > 1.e-9 )
THEN 2198 myres(k) = ( fres5 - 8.0*fres4 + 8.0*fres2 - fres1 ) / ( 6.0*delta_rho )
2200 myres(k) = ( -3.0*fres3 + 4.0*fres2 - fres1 ) / delta_rho
2217 IF (ensemble_flag ==
'tp') lnphi(k) = myres(k) - log(zges)
2218 IF (ensemble_flag ==
'tv' .AND. eta >= 0.0) lnphi(k) = myres(k)
2226 END SUBROUTINE phi_numerical
2269 REAL :: zges, df_dt, dfdr, ddfdrdr
2270 REAL :: cv_res, df_dt2, df_drdt
2271 REAL :: fact, dist, t_tmp, rho_0
2272 REAL :: fdr1, fdr2, fdr3, fdr4
2275 REAL :: dhsdt(nc), dhsdt2(nc)
2276 REAL :: z0, z1, z2, z3, z1tdt, z2tdt, z3tdt
2277 REAL :: z1dt, z2dt, z3dt, zms, gii
2278 REAL :: fhsdt, fhsdt2
2279 REAL :: fchdt, fchdt2
2280 REAL :: fdspdt, fdspdt2
2281 REAL :: fhbdt, fhbdt2
2282 REAL :: sumseg, i1, i2, i1dt, i2dt, i1dt2, i2dt2
2283 REAL :: c1_con, c2_con, c3_con, c1_dt, c1_dt2
2284 REAL :: z1tdt2, z2tdt2, z3tdt2
2285 REAL :: z1dt2, z2dt2, z3dt2
2287 INTEGER :: j, k, l, no, ass_cnt, max_eval
2289 REAL :: dij, dijdt, dijdt2
2290 REAL :: gij1dt, gij2dt, gij3dt
2291 REAL :: gij1dt2, gij2dt2, gij3dt2
2292 REAL,
DIMENSION(nc,nc) :: gijdt, gijdt2, kap_hb
2293 REAL,
DIMENSION(nc,nc,nsite,nsite) :: ass_d_dt, ass_d_dt2, eps_hb, delta, deltadt, deltadt2
2294 REAL,
DIMENSION(nc,nsite) :: mxdt, mxdt2, mx_itr, mx_itrdt, mx_itrdt2
2295 REAL :: attenu, tol, suma, sumdt, sumdt2, err_sum
2298 REAL :: fdddt, fdddt2
2299 REAL,
DIMENSION(nc) :: my2dd, my0
2300 REAL,
DIMENSION(nc,nc) :: idd2, idd2dt, idd2dt2, idd4, idd4dt, idd4dt2
2301 REAL,
DIMENSION(nc,nc,nc) :: idd3, idd3dt, idd3dt2
2302 REAL :: factor2, factor3
2303 REAL :: fdd2, fdd3, fdd2dt, fdd3dt, fdd2dt2, fdd3dt2
2304 REAL :: eij, xijmt, xijkmt
2307 REAL :: fqqdt, fqqdt2
2308 REAL,
DIMENSION(nc) :: qq2
2309 REAL,
DIMENSION(nc,nc) :: iqq2, iqq2dt, iqq2dt2, iqq4, iqq4dt, iqq4dt2
2310 REAL,
DIMENSION(nc,nc,nc) :: iqq3, iqq3dt, iqq3dt2
2311 REAL :: fqq2, fqq2dt, fqq2dt2, fqq3, fqq3dt, fqq3dt2
2314 REAL :: fdqdt, fdqdt2
2315 REAL,
DIMENSION(nc) :: myfac, q_fac
2316 REAL,
DIMENSION(nc,nc) :: idq2, idq2dt, idq2dt2, idq4, idq4dt, idq4dt2
2317 REAL,
DIMENSION(nc,nc,nc) :: idq3, idq3dt, idq3dt2
2318 REAL :: fdq2, fdq2dt, fdq2dt2, fdq3, fdq3dt, fdq3dt2
2325 CALL perturbation_parameter
2333 sumseg = z0t / (pi/6.0)
2342 zges = (pges * 1.e-30)/(kbol*t*rho)
2344 dfdr = pges/(eta*rho*(kbol*t)/1.e-30)
2345 ddfdrdr = pgesdz/(eta*rho*(kbol*t)/1.e-30) - dfdr*2.0/eta - 1.0/eta**2
2358 dhsdt(i)=parame(i,2) *(-3.0*parame(i,3)/t/t)*0.12*exp(-3.0*parame(i,3)/t)
2359 dhsdt2(i) = dhsdt(i)*3.0*parame(i,3)/t/t &
2360 + 6.0*parame(i,2)*parame(i,3)/t**3 *0.12*exp(-3.0*parame(i,3)/t)
2367 z1tdt = z1tdt + x(i) * mseg(i) * dhsdt(i)
2368 z2tdt = z2tdt + x(i) * mseg(i) * 2.0*dhs(i)*dhsdt(i)
2369 z3tdt = z3tdt + x(i) * mseg(i) * 3.0*dhs(i)*dhs(i)*dhsdt(i)
2371 z1dt = pi / 6.0*z1tdt *rho
2372 z2dt = pi / 6.0*z2tdt *rho
2373 z3dt = pi / 6.0*z3tdt *rho
2380 z1tdt2 = z1tdt2 + x(i)*mseg(i)*dhsdt2(i)
2381 z2tdt2 = z2tdt2 + x(i)*mseg(i)*2.0 *( dhsdt(i)*dhsdt(i) +dhs(i)*dhsdt2(i) )
2382 z3tdt2 = z3tdt2 + x(i)*mseg(i)*3.0 *( 2.0*dhs(i)*dhsdt(i)* &
2383 dhsdt(i) +dhs(i)*dhs(i)*dhsdt2(i) )
2385 z1dt2 = pi / 6.0*z1tdt2 *rho
2386 z2dt2 = pi / 6.0*z2tdt2 *rho
2387 z3dt2 = pi / 6.0*z3tdt2 *rho
2393 fhsdt = 6.0/pi/rho*( 3.0*(z1dt*z2+z1*z2dt)/zms + 3.0*z1*z2*z3dt/zms/zms &
2394 + 3.0*z2*z2*z2dt/z3/zms/zms &
2395 + z2**3 *(2.0*z3*z3dt-z3dt*zms)/(z3*z3*zms**3 ) &
2396 + (3.0*z2*z2*z2dt*z3-2.0*z2**3 *z3dt)/z3**3 *log(zms) &
2397 + (z0-z2**3 /z3/z3)*z3dt/zms )
2399 fhsdt2= 6.0/pi/rho*( 3.0*(z1dt2*z2+2.0*z1dt*z2dt+z1*z2dt2)/zms &
2400 + 6.0*(z1dt*z2+z1*z2dt)*z3dt/zms/zms &
2401 + 3.0*z1*z2*z3dt2/zms/zms + 6.0*z1*z2*z3dt*z3dt/zms**3 &
2402 + 3.0*z2*(2.0*z2dt*z2dt+z2*z2dt2)/z3/zms/zms &
2403 - z2*z2*(6.0*z2dt*z3dt+z2*z3dt2)/(z3*z3*zms*zms) &
2404 + 2.0*z2**3 *z3dt*z3dt/(z3**3 *zms*zms) &
2405 - 4.0*z2**3 *z3dt*z3dt/(z3*z3 *zms**3 ) &
2406 + (12.0*z2*z2*z2dt*z3dt+2.0*z2**3 *z3dt2)/(z3*zms**3 ) &
2407 + 6.0*z2**3 *z3dt*z3dt/(z3*zms**4 ) &
2408 - 2.0*(3.0*z2*z2*z2dt/z3/z3-2.0*z2**3 *z3dt/z3**3 ) *z3dt/zms &
2409 -(z2**3 /z3/z3-z0)*(z3dt2/zms+z3dt*z3dt/zms/zms) &
2410 + ( (6.0*z2*z2dt*z2dt+3.0*z2*z2*z2dt2)/z3/z3 &
2411 - (12.0*z2*z2*z2dt*z3dt+2.0*z2**3 *z3dt2)/z3**3 &
2412 + 6.0*z2**3 *z3dt*z3dt/z3**4 )* log(zms) )
2422 dij=dhs(i)*dhs(j)/(dhs(i)+dhs(j))
2423 dijdt =(dhsdt(i)*dhs(j) + dhs(i)*dhsdt(j)) / (dhs(i)+dhs(j)) &
2424 - dhs(i)*dhs(j)/(dhs(i)+dhs(j))**2 *(dhsdt(i)+dhsdt(j))
2425 dijdt2=(dhsdt2(i)*dhs(j) + 2.0*dhsdt(i)*dhsdt(j) &
2426 + dhs(i)*dhsdt2(j)) / (dhs(i)+dhs(j)) &
2427 - 2.0*(dhsdt(i)*dhs(j) + dhs(i)*dhsdt(j)) &
2428 / (dhs(i)+dhs(j))**2 *(dhsdt(i)+dhsdt(j)) &
2429 + 2.0* dhs(i)*dhs(j) / (dhs(i)+dhs(j))**3 &
2430 * (dhsdt(i)+dhsdt(j))**2 &
2431 - dhs(i)*dhs(j)/(dhs(i)+dhs(j))**2 *(dhsdt2(i)+dhsdt2(j))
2432 gij1dt = z3dt/zms/zms
2433 gij2dt = 3.0*( z2dt*dij+z2*dijdt )/zms/zms +6.0*z2*dij*z3dt/zms**3
2434 gij3dt = 4.0*(dij*z2)* (dijdt*z2 + dij*z2dt) /zms**3 &
2435 + 6.0*(dij*z2)**2 * z3dt /zms**4
2436 gij1dt2 = z3dt2/zms/zms +2.0*z3dt*z3dt/zms**3
2437 gij2dt2 = 3.0*( z2dt2*dij+2.0*z2dt*dijdt+z2*dijdt2 )/zms/zms &
2438 + 6.0*( z2dt*dij+z2*dijdt )/zms**3 * z3dt &
2439 + 6.0*(z2dt*dij*z3dt+z2*dijdt*z3dt+z2*dij*z3dt2) /zms**3 &
2440 + 18.0*z2*dij*z3dt*z3dt/zms**4
2441 gij3dt2 = 4.0*(dijdt*z2+dij*z2dt)**2 /zms**3 &
2442 + 4.0*(dij*z2)* (dijdt2*z2+2.0*dijdt*z2dt+dij*z2dt2) /zms**3 &
2443 + 24.0*(dij*z2) *(dijdt*z2+dij*z2dt)/zms**4 *z3dt &
2444 + 6.0*(dij*z2)**2 * z3dt2 /zms**4 &
2445 + 24.0*(dij*z2)**2 * z3dt*z3dt /zms**5
2446 gijdt(i,j) = gij1dt + gij2dt + gij3dt
2447 gijdt2(i,j) = gij1dt2 + gij2dt2 + gij3dt2
2452 gii = 1.0/zms + 3.0*dhs(i)/2.0*z2/zms/zms + 2.0*dhs(i)*dhs(i)/4.0*z2*z2/zms**3
2453 fchdt = fchdt + x(i) * (1.0-mseg(i)) * gijdt(i,i) / gii
2454 fchdt2= fchdt2+ x(i) * (1.0-mseg(i)) &
2455 * (gijdt2(i,i) / gii - (gijdt(i,i)/gii)**2 )
2469 i1 = i1 + apar(m)*z3**
REAL(m)
2470 i2 = i2 + bpar(m)*z3**
REAL(m)
2471 i1dt = i1dt + apar(m)*z3dt*
REAL(m)*z3**
REAL(m-1)
2472 i2dt = i2dt + bpar(m)*z3dt*
REAL(m)*z3**
REAL(m-1)
2473 i1dt2= i1dt2+ apar(m)*z3dt2*
REAL(m)*z3**
REAL(m-1) &
2474 + apar(m)*z3dt*z3dt *
REAL(m)*
REAL(m-1)*z3**
REAL(m-2)
2475 i2dt2= i2dt2+ bpar(m)*z3dt2*
REAL(m)*z3**
REAL(m-1) &
2476 + bpar(m)*z3dt*z3dt *
REAL(m)*
REAL(m-1)*z3**
REAL(m-2)
2479 c1_con= 1.0/ ( 1.0 + sumseg*(8.0*z3-2.0*z3**2 )/zms**4 &
2480 + (1.0 - sumseg)*(20.0*z3-27.0*z3**2 +12.0*z3**3 -2.0*z3**4 ) &
2481 /(zms*(2.0-z3))**2 )
2482 c2_con= - c1_con*c1_con *(sumseg*(-4.0*z3**2 +20.0*z3+8.0)/zms**5 &
2483 + (1.0 - sumseg) *(2.0*z3**3 +12.0*z3**2 -48.0*z3+40.0) &
2484 /(zms*(2.0-z3))**3 )
2485 c3_con= 2.0 * c2_con*c2_con/c1_con - c1_con*c1_con &
2486 *( sumseg*(-12.0*z3**2 +72.0*z3+60.0)/zms**6 + (1.0 - sumseg) &
2487 *(-6.0*z3**4 -48.0*z3**3 +288.0*z3**2 -480.0*z3+264.0) &
2488 /(zms*(2.0-z3))**4 )
2490 c1_dt2 = c3_con*z3dt*z3dt + c2_con*z3dt2
2492 fdspdt = - 2.0*pi*rho*(i1dt-i1/t)*order1 &
2493 - pi*rho*sumseg*(c1_dt*i2+c1_con*i2dt-2.0*c1_con*i2/t)*order2
2495 fdspdt2 = - 2.0*pi*rho*(i1dt2-2.0*i1dt/t+2.0*i1/t/t)*order1 &
2496 - pi*rho*sumseg*order2*( c1_dt2*i2 +2.0*c1_dt*i2dt -4.0*c1_dt*i2/t &
2497 + 6.0*c1_con*i2/t/t -4.0*c1_con*i2dt/t +c1_con*i2dt2)
2507 IF ( nhb_typ(i) /= 0 ) assoc = .true.
2512 IF ( nhb_typ(i) /= 0 )
THEN 2513 kap_hb(i,i) = parame(i,13)
2517 eps_hb(i,i,j,k) = parame(i,(14+no))
2528 eps_hb(i,i,k,l) = 0.0
2536 IF ( i /= j .AND. (nhb_typ(i) /= 0 .AND. nhb_typ(j) /= 0) )
THEN 2537 kap_hb(i,j)= (kap_hb(i,i)*kap_hb(j,j))**0.5 &
2538 *((parame(i,2)*parame(j,2))**3 )**0.5 &
2539 /(0.5*(parame(i,2)+parame(j,2)))**3
2544 eps_hb(i,j,k,l)=(eps_hb(i,i,k,l)+eps_hb(j,j,l,k))/2.0
2552 IF (nhb_typ(1) == 3)
THEN 2553 eps_hb(1,2,3,1)=0.5*(eps_hb(1,1,3,1)+eps_hb(2,2,1,2))
2554 eps_hb(2,1,1,3) = eps_hb(1,2,3,1)
2556 IF (nhb_typ(2) == 3)
THEN 2557 eps_hb(2,1,3,1)=0.5*(eps_hb(2,2,3,1)+eps_hb(1,1,1,2))
2558 eps_hb(1,2,1,3) = eps_hb(2,1,3,1)
2562 DO k = 1, nhb_typ(i)
2564 DO l = 1, nhb_typ(j)
2566 ass_d_dt(i,j,k,l)= - kap_hb(i,j) *sig_ij(i,j)**3 * eps_hb(i,j,k,l)/t/t*exp(eps_hb(i,j,k,l)/t)
2567 ass_d_dt2(i,j,k,l)= - kap_hb(i,j) *sig_ij(i,j)**3 &
2568 * eps_hb(i,j,k,l)/t/t*exp(eps_hb(i,j,k,l)/t) &
2569 * (-2.0/t - eps_hb(i,j,k,l)/t/t)
2576 DO k = 1, nhb_typ(i)
2578 DO l = 1, nhb_typ(j)
2579 delta(i,j,k,l)=gij(i,j)*ass_d(i,j,k,l)
2580 deltadt(i,j,k,l) = gijdt(i,j)*ass_d(i,j,k,l) + gij(i,j)*ass_d_dt(i,j,k,l)
2581 deltadt2(i,j,k,l)= gijdt2(i,j)*ass_d(i,j,k,l) &
2582 + 2.0*gijdt(i,j)*ass_d_dt(i,j,k,l) +gij(i,j)*ass_d_dt2(i,j,k,l)
2592 IF (eta < 0.2) tol = 1.e-11
2593 IF (eta < 0.01) tol = 1.e-12
2598 DO k = 1, nhb_typ(i)
2606 DO ass_cnt = 1, max_eval
2609 DO k = 1, nhb_typ(i)
2614 DO l = 1, nhb_typ(j)
2615 suma = suma + x(j)*nhb_no(j,l)* mx(j,l) *delta(i,j,k,l)
2616 sumdt = sumdt + x(j)*nhb_no(j,l)*( mx(j,l) *deltadt(i,j,k,l) &
2617 + mxdt(j,l)*delta(i,j,k,l) )
2618 sumdt2 = sumdt2 + x(j)*nhb_no(j,l)*( mx(j,l)*deltadt2(i,j,k,l) &
2619 + 2.0*mxdt(j,l)*deltadt(i,j,k,l) + mxdt2(j,l)*delta(i,j,k,l) )
2622 mx_itr(i,k) = 1.0 / (1.0 + suma * rho)
2623 mx_itrdt(i,k)= - mx_itr(i,k)**2 * sumdt*rho
2624 mx_itrdt2(i,k)= +2.0*mx_itr(i,k)**3 * (sumdt*rho)**2 - mx_itr(i,k)**2 *sumdt2*rho
2630 DO k = 1, nhb_typ(i)
2631 err_sum = err_sum + abs(mx_itr(i,k) - mx(i,k)) &
2632 + abs(mx_itrdt(i,k) - mxdt(i,k)) + abs(mx_itrdt2(i,k) - mxdt2(i,k))
2633 mx(i,k) = mx_itr(i,k) * attenu + mx(i,k) * (1.0 - attenu)
2634 mxdt(i,k)=mx_itrdt(i,k)*attenu +mxdt(i,k)* (1.0 - attenu)
2635 mxdt2(i,k)=mx_itrdt2(i,k)*attenu +mxdt2(i,k)* (1.0 - attenu)
2638 IF(err_sum <= tol)
GO TO 10
2641 WRITE(6,*)
'CAL_PCSAFT: max_eval violated err_sum = ',err_sum,tol
2646 DO k = 1, nhb_typ(i)
2648 fhbdt = fhbdt + x(i)*nhb_no(i,k) *mxdt(i,k)*(1.0/mx(i,k)-0.5)
2649 fhbdt2= fhbdt2 + x(i)*nhb_no(i,k) *(mxdt2(i,k)*(1.0/mx(i,k)-0.5) &
2650 -(mxdt(i,k)/mx(i,k))**2 )
2665 IF ( parame(i,6) /= 0.0 .AND. uij(i,i) /= 0.0 )
THEN 2667 my2dd(i) = (parame(i,6))**2 *1.e-49 / (uij(i,i)*kbol*mseg(i)*sig_ij(i,i)**3 *1.e-30)
2672 IF (dipole == 1)
THEN 2682 idd2(i,j) = idd2(i,j) +ddp2(i,j,m)*z3**
REAL(m)
2683 idd4(i,j) = idd4(i,j) +ddp4(i,j,m)*z3**
REAL(m)
2684 idd2dt(i,j)= idd2dt(i,j) +ddp2(i,j,m)*z3dt*
REAL(m)*z3**
REAL(m-1)
2685 idd4dt(i,j)= idd4dt(i,j) +ddp4(i,j,m)*z3dt*
REAL(m)*z3**
REAL(m-1)
2686 idd2dt2(i,j)=idd2dt2(i,j)+ddp2(i,j,m)*z3dt2*
REAL(m)*z3**
REAL(m-1) &
2687 + ddp2(i,j,m)*z3dt**2 *
REAL(m)*
REAL(m-1)*z3**
REAL(m-2)
2688 idd4dt2(i,j)=idd4dt2(i,j)+ddp4(i,j,m)*z3dt2*
REAL(m)*z3**
REAL(m-1) &
2689 + ddp4(i,j,m)*z3dt**2 *
REAL(m)*
REAL(m-1)*z3**
REAL(m-2)
2696 idd3(i,j,k) = idd3(i,j,k) +ddp3(i,j,k,m)*z3**
REAL(m)
2697 idd3dt(i,j,k) = idd3dt(i,j,k)+ddp3(i,j,k,m)*z3dt*
REAL(m) *z3**
REAL(m-1)
2698 idd3dt2(i,j,k)= idd3dt2(i,j,k)+ddp3(i,j,k,m)*z3dt2*
REAL(m) &
2699 *z3**
REAL(m-1) +ddp3(i,j,k,m)*z3dt**2 *
REAL(m)*
REAL(m-1) *z3**
REAL(m-2)
2707 factor3= -4.0/3.0*pi**2 * rho**2
2717 xijmt = x(i)*parame(i,3)*parame(i,2)**3 *x(j)*parame(j,3)*parame(j,2)**3 &
2718 / ((parame(i,2)+parame(j,2))/2.0)**3 *my2dd(i)*my2dd(j)
2719 eij = (parame(i,3)*parame(j,3))**0.5
2720 fdd2 = fdd2 +factor2* xijmt/t/t*(idd2(i,j)+eij/t*idd4(i,j))
2721 fdd2dt= fdd2dt+ factor2* xijmt/t/t*(idd2dt(i,j)-2.0*idd2(i,j)/t &
2722 +eij/t*idd4dt(i,j)-3.0*eij/t/t*idd4(i,j))
2723 fdd2dt2=fdd2dt2+factor2*xijmt/t/t*(idd2dt2(i,j)-4.0*idd2dt(i,j)/t &
2724 +6.0*idd2(i,j)/t/t+eij/t*idd4dt2(i,j) &
2725 -6.0*eij/t/t*idd4dt(i,j)+12.0*eij/t**3 *idd4(i,j))
2727 xijkmt=x(i)*parame(i,3)*parame(i,2)**3 &
2728 *x(j)*parame(j,3)*parame(j,2)**3 &
2729 *x(k)*parame(k,3)*parame(k,2)**3 &
2730 /((parame(i,2)+parame(j,2))/2.0) /((parame(i,2)+parame(k,2))/2.0) &
2731 /((parame(j,2)+parame(k,2))/2.0) *my2dd(i)*my2dd(j)*my2dd(k)
2732 fdd3 =fdd3 +factor3*xijkmt/t**3 *idd3(i,j,k)
2733 fdd3dt =fdd3dt +factor3*xijkmt/t**3 * (idd3dt(i,j,k)-3.0*idd3(i,j,k)/t)
2734 fdd3dt2=fdd3dt2+factor3*xijkmt/t**3 &
2735 *( idd3dt2(i,j,k)-6.0*idd3dt(i,j,k)/t+12.0*idd3(i,j,k)/t/t )
2740 IF ( fdd2 < -1.e-100 .AND. fdd3 /= 0.0 )
THEN 2741 fdddt = fdd2* (fdd2*fdd2dt - 2.0*fdd3*fdd2dt+fdd2*fdd3dt) / (fdd2-fdd3)**2
2742 fdddt2 = ( 2.0*fdd2*fdd2dt*fdd2dt +fdd2*fdd2*fdd2dt2 &
2743 -2.0*fdd2dt**2 *fdd3 -2.0*fdd2*fdd2dt2*fdd3 +fdd2*fdd2*fdd3dt2 ) &
2744 /(fdd2-fdd3)**2 + fdddt * 2.0*(fdd3dt-fdd2dt)/(fdd2-fdd3)
2756 qq2(i) = (parame(i,7))**2 *1.e-69 / (uij(i,i)*kbol*mseg(i)*sig_ij(i,i)**5 *1.e-50)
2757 IF (qq2(i) /= 0.0) qudpole = 1
2760 IF (qudpole == 1)
THEN 2771 iqq2(i,j) = iqq2(i,j) + qqp2(i,j,m)*z3**
REAL(m)
2772 iqq4(i,j) = iqq4(i,j) + qqp4(i,j,m)*z3**
REAL(m)
2773 iqq2dt(i,j) = iqq2dt(i,j)+ qqp2(i,j,m)*z3dt*
REAL(m)*z3**
REAL(m-1)
2774 iqq4dt(i,j) = iqq4dt(i,j)+ qqp4(i,j,m)*z3dt*
REAL(m)*z3**
REAL(m-1)
2775 iqq2dt2(i,j)= iqq2dt2(i,j)+qqp2(i,j,m)*z3dt2*
REAL(m)*z3**
REAL(m-1) &
2776 + qqp2(i,j,m)*z3dt**2 *
REAL(m)*
REAL(m-1)*z3**
REAL(m-2)
2777 iqq4dt2(i,j)= iqq4dt2(i,j)+qqp4(i,j,m)*z3dt2*
REAL(m)*z3**
REAL(m-1) &
2778 + qqp4(i,j,m)*z3dt**2 *
REAL(m)*
REAL(m-1)*z3**
REAL(m-2)
2785 iqq3(i,j,k) = iqq3(i,j,k) + qqp3(i,j,k,m)*z3**
REAL(m)
2786 iqq3dt(i,j,k) = iqq3dt(i,j,k)+ qqp3(i,j,k,m)*z3dt*
REAL(m) * z3**
REAL(m-1)
2787 iqq3dt2(i,j,k)= iqq3dt2(i,j,k)+qqp3(i,j,k,m)*z3dt2*
REAL(m) * z3**
REAL(m-1) &
2788 + qqp3(i,j,k,m)*z3dt**2 *
REAL(m)*
REAL(m-1)*z3**
REAL(m-2)
2794 factor2 = -9.0/16.0 * pi *rho
2795 factor3 = 9.0/16.0 * pi**2 * rho**2
2805 xijmt = x(i)*uij(i,i)*qq2(i)*sig_ij(i,i)**5 &
2806 * x(j)*uij(j,j)*qq2(j)*sig_ij(j,j)**5 /sig_ij(i,j)**7.0
2807 eij = (parame(i,3)*parame(j,3))**0.5
2808 fqq2 = fqq2 +factor2* xijmt/t/t*(iqq2(i,j)+eij/t*iqq4(i,j))
2809 fqq2dt= fqq2dt +factor2* xijmt/t/t*(iqq2dt(i,j)-2.0*iqq2(i,j)/t &
2810 + eij/t*iqq4dt(i,j)-3.0*eij/t/t*iqq4(i,j))
2811 fqq2dt2=fqq2dt2+factor2*xijmt/t/t*(iqq2dt2(i,j)-4.0*iqq2dt(i,j)/t &
2812 + 6.0*iqq2(i,j)/t/t+eij/t*iqq4dt2(i,j) &
2813 - 6.0*eij/t/t*iqq4dt(i,j)+12.0*eij/t**3 *iqq4(i,j))
2815 xijkmt = x(i)*uij(i,i)*qq2(i)*sig_ij(i,i)**5 /sig_ij(i,j)**3 &
2816 * x(j)*uij(j,j)*qq2(j)*sig_ij(j,j)**5 /sig_ij(i,k)**3 &
2817 * x(k)*uij(k,k)*qq2(k)*sig_ij(k,k)**5 /sig_ij(j,k)**3
2818 fqq3 = fqq3 +factor3*xijkmt/t**3 *iqq3(i,j,k)
2819 fqq3dt = fqq3dt +factor3*xijkmt/t**3 *(iqq3dt(i,j,k)-3.0*iqq3(i,j,k)/t)
2820 fqq3dt2= fqq3dt2+factor3*xijkmt/t**3 &
2821 * ( iqq3dt2(i,j,k)-6.0*iqq3dt(i,j,k)/t+12.0*iqq3(i,j,k)/t/t )
2826 IF ( fqq2 /= 0.0 .AND. fqq3 /= 0.0 )
THEN 2827 fqqdt = fqq2* (fqq2*fqq2dt - 2.0*fqq3*fqq2dt+fqq2*fqq3dt) / (fqq2-fqq3)**2
2828 fqqdt2 = ( 2.0*fqq2*fqq2dt*fqq2dt +fqq2*fqq2*fqq2dt2 &
2829 - 2.0*fqq2dt**2 *fqq3 -2.0*fqq2*fqq2dt2*fqq3 +fqq2*fqq2*fqq3dt2 ) &
2830 / (fqq2-fqq3)**2 + fqqdt * 2.0*(fqq3dt-fqq2dt)/(fqq2-fqq3)
2844 IF (parame(i,6) /= 0.0 .AND. parame(j,7) /= 0.0) dip_quad = 1
2846 myfac(i) = parame(i,3)*parame(i,2)**4 *my0(i)
2847 q_fac(i) = parame(i,3)*parame(i,2)**4 *qq2(i)
2850 IF (dip_quad == 1)
THEN 2860 IF ( myfac(i) /= 0.0 .AND. q_fac(j) /= 0.0 )
THEN 2862 idq2(i,j) = idq2(i,j) + dqp2(i,j,m)*z3**
REAL(m)
2863 idq4(i,j) = idq4(i,j) + dqp4(i,j,m)*z3**
REAL(m)
2864 idq2dt(i,j) = idq2dt(i,j)+ dqp2(i,j,m)*z3dt*
REAL(m)*z3**
REAL(m-1)
2865 idq4dt(i,j) = idq4dt(i,j)+ dqp4(i,j,m)*z3dt*
REAL(m)*z3**
REAL(m-1)
2866 idq2dt2(i,j)= idq2dt2(i,j)+dqp2(i,j,m)*z3dt2*
REAL(m)*z3**
REAL(m-1) &
2867 + dqp2(i,j,m)*z3dt**2 *
REAL(m)*
REAL(m-1)*z3**
REAL(m-2)
2868 idq4dt2(i,j)= idq4dt2(i,j)+dqp4(i,j,m)*z3dt2*
REAL(m)*z3**
REAL(m-1) &
2869 + dqp4(i,j,m)*z3dt**2 *
REAL(m)*
REAL(m-1)*z3**
REAL(m-2)
2876 IF ( myfac(k) /= 0.0 .OR. q_fac(k) /= 0.0 )
THEN 2878 idq3(i,j,k) = idq3(i,j,k) + dqp3(i,j,k,m)*z3**
REAL(m)
2879 idq3dt(i,j,k)= idq3dt(i,j,k)+ dqp3(i,j,k,m)*z3dt*
REAL(m) *z3**
REAL(m-1)
2880 idq3dt2(i,j,k)= idq3dt2(i,j,k)+dqp3(i,j,k,m)*z3dt2*
REAL(m) *z3**
REAL(m-1) &
2881 + dqp3(i,j,k,m)*z3dt**2 *
REAL(m)*
REAL(m-1)*z3**
REAL(m-2)
2889 factor2= -9.0/4.0 * pi * rho
2890 factor3= pi**2 * rho**2
2900 IF ( myfac(i) /= 0.0 .AND. q_fac(j) /= 0.0 )
THEN 2901 xijmt = x(i)*myfac(i) * x(j)*q_fac(j) /sig_ij(i,j)**5
2902 eij = (parame(i,3)*parame(j,3))**0.5
2903 fdq2 = fdq2 + factor2* xijmt/t/t*(idq2(i,j)+eij/t*idq4(i,j))
2904 fdq2dt= fdq2dt+ factor2* xijmt/t/t*(idq2dt(i,j)-2.0*idq2(i,j)/t &
2905 + eij/t*idq4dt(i,j)-3.0*eij/t/t*idq4(i,j))
2906 fdq2dt2 = fdq2dt2+factor2*xijmt/t/t*(idq2dt2(i,j)-4.0*idq2dt(i,j)/t &
2907 + 6.0*idq2(i,j)/t/t+eij/t*idq4dt2(i,j) &
2908 - 6.0*eij/t/t*idq4dt(i,j)+12.0*eij/t**3 *idq4(i,j))
2910 IF ( myfac(k) /= 0.0 .OR. q_fac(k) /= 0.0 )
THEN 2911 xijkmt= x(i)*x(j)*x(k)/(sig_ij(i,j)*sig_ij(i,k)*sig_ij(j,k))**2 &
2912 * ( myfac(i)*q_fac(j)*myfac(k) &
2913 + myfac(i)*q_fac(j)*q_fac(k)*1.193735 )
2915 fdq3 =fdq3 + factor3*xijkmt/t**3 *idq3(i,j,k)
2916 fdq3dt=fdq3dt+ factor3*xijkmt/t**3 * (idq3dt(i,j,k)-3.0*idq3(i,j,k)/t)
2917 fdq3dt2=fdq3dt2+factor3*xijkmt/t**3 &
2918 *( idq3dt2(i,j,k)-6.0*idq3dt(i,j,k)/t+12.0*idq3(i,j,k)/t/t )
2925 IF (fdq2 /= 0.0 .AND. fdq3 /= 0.0)
THEN 2926 fdqdt = fdq2* (fdq2*fdq2dt - 2.0*fdq3*fdq2dt+fdq2*fdq3dt) / (fdq2-fdq3)**2
2927 fdqdt2 = ( 2.0*fdq2*fdq2dt*fdq2dt +fdq2*fdq2*fdq2dt2 &
2928 - 2.0*fdq2dt**2 *fdq3 -2.0*fdq2*fdq2dt2*fdq3 +fdq2*fdq2*fdq3dt2 ) &
2929 / (fdq2-fdq3)**2 + fdqdt * 2.0*(fdq3dt-fdq2dt)/(fdq2-fdq3)
2942 df_dt = fhsdt + fchdt + fdspdt + fhbdt + fdddt + fqqdt + fdqdt
2950 df_dt2 = fhsdt2 +fchdt2 +fdspdt2 +fhbdt2 +fdddt2 +fqqdt2 +fdqdt2
2963 dist = t * 100.e-5 * fact
2968 t = t_tmp - 2.0*dist
2969 CALL perturbation_parameter
2972 fdr1 = pges / (eta*rho_0*(kbol*t)/1.e-30)
2974 CALL perturbation_parameter
2977 fdr2 = pges / (eta*rho_0*(kbol*t)/1.e-30)
2980 CALL perturbation_parameter
2983 fdr3 = pges / (eta*rho_0*(kbol*t)/1.e-30)
2985 t = t_tmp + 2.0*dist
2986 CALL perturbation_parameter
2989 fdr4 = pges / (eta*rho_0*(kbol*t)/1.e-30)
2992 CALL perturbation_parameter
2997 df_drdt = (-fdr4+8.0*fdr3-8.0*fdr2+fdr1)/(12.0*dist)
3007 s_res = ( - df_dt *t - fres )*rgas + rgas * log(zges)
3008 h_res = ( - t*df_dt + zges-1.0 ) * rgas *t
3009 cv_res = - ( t*df_dt2 + 2.0*df_dt ) * rgas*t
3010 cp_res = cv_res - rgas + rgas*(zges +eta*t*df_drdt)**2 &
3011 / (1.0 + 2.0*eta*dfdr +eta**2 *ddfdrdr)
3018 END SUBROUTINE h_eos
3029 SUBROUTINE h_eos_num
3036 REAL :: dist, fact, rho_0
3037 REAL :: fres1, fres2, fres3, fres4, fres5
3038 REAL :: f_1, f_2, f_3, f_4
3039 REAL :: cv_res, t_tmp, zges
3040 REAL :: df_dt, df_dtdt, df_drdt, dfdr, ddfdrdr
3045 CALL perturbation_parameter
3050 dist = t * 100.e-5 * fact
3054 t = t_tmp - 2.0*dist
3055 CALL perturbation_parameter
3060 CALL perturbation_parameter
3065 CALL perturbation_parameter
3069 t = t_tmp + 2.0*dist
3070 CALL perturbation_parameter
3075 CALL perturbation_parameter
3081 zges = (p * 1.e-30)/(kbol*t*rho_0)
3084 df_dt = (-fres4+8.0*fres3-8.0*fres2+fres1)/(12.0*dist)
3085 df_dtdt = (-fres4+16.0*fres3-3.d1*fres5+16.0*fres2-fres1) &
3089 s_res = (- df_dt -fres/t)*rgas*t + rgas * log(zges)
3090 h_res = ( - t*df_dt + zges-1.0 ) * rgas*t
3091 cv_res = -( t*df_dtdt + 2.0*df_dt ) * rgas*t
3095 t = t_tmp - 2.0*dist
3096 CALL perturbation_parameter
3099 f_1 = pges/(eta*rho_0*(kbol*t)/1.e-30)
3102 CALL perturbation_parameter
3105 f_2 = pges/(eta*rho_0*(kbol*t)/1.e-30)
3108 CALL perturbation_parameter
3111 f_3 = pges/(eta*rho_0*(kbol*t)/1.e-30)
3113 t = t_tmp + 2.0*dist
3114 CALL perturbation_parameter
3117 f_4 = pges/(eta*rho_0*(kbol*t)/1.e-30)
3120 CALL perturbation_parameter
3124 dfdr = pges / (eta*rho_0*(kbol*t)/1.e-30)
3125 ddfdrdr = pgesdz/(eta*rho_0*(kbol*t)/1.e-30) - dfdr*2.0/eta - 1.0/eta**2
3127 df_drdt = ( -f_4 +8.0*f_3 -8.0*f_2 +f_1) / (12.0*dist)
3129 cp_res = cv_res - rgas +rgas*(zges+eta*t*df_drdt)**2 &
3130 * 1.0/(1.0 + 2.0*eta*dfdr + eta**2 *ddfdrdr)
3139 END SUBROUTINE h_eos_num
3150 SUBROUTINE density_iteration
3157 INTEGER :: i, start, max_i
3158 REAL :: eta_iteration
3159 REAL :: error, dydx, acc_i, delta_eta
3163 IF ( densav(phas) /= 0.0 .AND. eta_start == denold(phas) )
THEN 3164 denold(phas) = eta_start
3165 eta_start = densav(phas)
3167 denold(phas) = eta_start
3168 densav(phas) = eta_start
3174 density_error(:) = 0.0
3177 eta_iteration = eta_start
3187 IF ( num == 0 )
THEN 3189 ELSE IF ( num == 1 )
THEN 3191 ELSE IF ( num == 2 )
THEN 3192 WRITE(*,*)
'CRITICAL RENORM NOT INCLUDED YET' 3195 write (*,*)
'define calculation option (num)' 3198 error = (pges / p ) - 1.0
3201 IF ( pgesdz < 0.0 .AND. i < max_i )
THEN 3202 IF ( error > 0.0 .AND. pgesd2 > 0.0 )
THEN 3203 CALL pressure_spinodal
3205 error = (pges / p ) - 1.0
3206 IF ( ((pges/p ) -1.0) > 0.0 ) eta_iteration = 0.001
3207 IF ( ((pges/p ) -1.0) <=0.0 ) eta_iteration = eta_iteration * 1.1
3208 ELSE IF ( error < 0.0 .AND. pgesd2 < 0.0 )
THEN 3209 CALL pressure_spinodal
3211 error = (pges / p ) - 1.0
3212 IF ( ((pges/p ) -1.0) < 0.0 ) eta_iteration = 0.5
3213 IF ( ((pges/p ) -1.0) >=0.0 ) eta_iteration = eta_iteration * 0.9
3215 eta_iteration = (eta_iteration + eta_start) / 2.0
3216 IF (eta_iteration == eta_start) eta_iteration = eta_iteration + 0.2
3218 cycle iterate_density
3223 delta_eta = error/ dydx
3224 IF ( delta_eta > 0.05 ) delta_eta = 0.05
3225 IF ( delta_eta < -0.05 ) delta_eta = -0.05
3227 eta_iteration = eta_iteration - delta_eta
3229 IF (eta_iteration > 0.9) eta_iteration = 0.6
3230 IF (eta_iteration <= 0.0) eta_iteration = 1.e-16
3233 IF ( abs(error*p/pgesdz) < 1.d-12 ) start = 0
3234 IF ( abs(error) < acc_i ) start = 0
3235 IF ( i > max_i )
THEN 3237 density_error(phas) = abs( error )
3241 IF (start /= 1)
EXIT iterate_density
3243 END DO iterate_density
3247 IF ((eta > 0.3 .AND. densav(phas) > 0.3) .OR. &
3248 (eta < 0.1 .AND. densav(phas) < 0.1)) densav(phas) = eta
3250 END SUBROUTINE density_iteration
3270 INTEGER :: i, j, k, l, m, n
3271 REAL :: z0, z1, z2, z3
3273 REAL :: i1,i2, c1_con
3274 REAL :: fhs, fdsp, fhc
3277 INTEGER :: ass_cnt,max_eval
3278 REAL :: delta(nc,nc,nsite,nsite)
3279 REAL :: mx_itr(nc,nsite), err_sum, sum, attenu, tol, fhb
3280 REAL :: ass_s1, ass_s2
3282 REAL :: fdd, fqq, fdq
3295 m_mean = z0t / ( pi / 6.0 )
3313 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
3321 fhs= m_mean*( 3.0*z1*z2/zms + z2**3 /z3/zms/zms + (z2**3 /z3/z3-z0)*log(zms) )/z0
3329 fhc = fhc + x(i) *(1.0- mseg(i)) *log(gij(i,i))
3341 i1 = i1 + apar(m)*eta**
REAL(m)
3342 i2 = i2 + bpar(m)*eta**
REAL(m)
3345 c1_con= 1.0/ ( 1.0 + m_mean*(8.0*eta-2.0*eta**2 )/zms**4 &
3346 + (1.0 - m_mean)*(20.0*eta-27.0*eta**2 &
3347 + 12.0*eta**3 -2.0*eta**4 ) /(zms*(2.0-eta))**2 )
3349 fdsp = -2.0*pi*rho*i1*order1 - pi*rho*c1_con*m_mean*i2*order2
3359 fdsp = fdsp + dnm(n,m) * (um/t)**
REAL(n) *(eta/tau)**
REAL(m)
3362 fdsp = m_mean * fdsp
3373 IF (nhb_typ(i) /= 0) assoc = .true.
3378 DO k = 1, nhb_typ(i)
3379 IF (mx(i,k) == 0.0) mx(i,k) = 1.0
3381 DO l = 1, nhb_typ(j)
3382 delta(i,j,k,l) = gij(i,j) * ass_d(i,j,k,l)
3392 IF (eta < 0.2) tol = 1.d-12
3393 IF (eta < 0.01) tol = 1.d-13
3400 ass_cnt = ass_cnt + 1
3403 DO k = 1, nhb_typ(i)
3406 DO l = 1, nhb_typ(j)
3407 sum = sum + x(j)* mx(j,l)*nhb_no(j,l) *delta(i,j,k,l)
3411 mx_itr(i,k) = 1.0 / (1.0 + sum * rho)
3418 DO k = 1, nhb_typ(i)
3419 err_sum = err_sum + abs(mx_itr(i,k) - mx(i,k))
3420 mx(i,k) = mx_itr(i,k) * attenu + mx(i,k) * (1.0 - attenu)
3424 IF ( err_sum <= tol .OR. ass_cnt >= max_eval )
THEN 3425 IF (ass_cnt >= max_eval)
WRITE(*,
'(a,2G15.7)')
'F_EOS: Max_eval violated (mx) Err_Sum = ',err_sum,tol
3434 DO k = 1, nhb_typ(i)
3435 ass_s1 = ass_s1 + nhb_no(i,k) * ( 1.0 - mx(i,k) )
3436 ass_s2 = ass_s2 + nhb_no(i,k) * log( mx(i,k) )
3438 fhb = fhb + x(i) * ( ass_s2 + ass_s1 / 2.0 )
3448 CALL f_polar ( fdd, fqq, fdq )
3454 fres = fhs + fhc + fdsp + fhb + fdd + fqq + fdq
3458 END SUBROUTINE f_eos
3465 SUBROUTINE f_numerical
3470 disp_term, hb_term, lc_term, branch_term, &
3471 ii_term, id_term, subtract1, subtract2
3479 REAL :: fid, fhs, fdsp, fhc
3480 REAL :: fhb, fdd, fqq, fdq
3485 REAL :: eps_kij, k_kij
3504 CALL perturbation_parameter
3515 order1 = order1 + x(i)*x(j)* mseg(i)*mseg(j) * sig_ij(i,j)**3 * uij(i,j)/t
3516 order2 = order2 + x(i)*x(j)* mseg(i)*mseg(j) * sig_ij(i,j)**3 * (uij(i,j)/t)**2
3521 order1 = order1 + x(i)*mseg(i)/t*( x(j)*mseg(j) &
3522 *sig_ij(i,j)*(uij(i,i)*uij(j,j))**(1.0/6.0) )**3 *lij(i,j)
3535 m_mean2 = m_mean2 + x(i)*x(j)*(mseg(i)+mseg(j))/2.0
3545 IF ( ideal_gas ==
'yes' )
CALL f_ideal_gas ( fid )
3549 IF ( hard_sphere ==
'CSBM' )
CALL f_hard_sphere ( m_mean2, fhs )
3553 IF ( chain_term ==
'TPT1' )
CALL f_chain_tpt1 ( fhc )
3554 IF ( chain_term ==
'TPT2' )
CALL f_chain_tpt_d ( fhc )
3555 IF ( chain_term ==
'HuLiu' )
CALL f_chain_hu_liu ( fhc )
3556 IF ( chain_term ==
'HuLiu_rc' )
CALL f_chain_hu_liu_rc ( fhs, fhc )
3558 IF ( chain_term ==
'SPT' )
WRITE(*,*)
'SPT NOT INCLUDED YET' 3562 IF ( disp_term ==
'PC-SAFT')
CALL f_disp_pcsaft ( fdsp )
3563 IF ( disp_term ==
'CK')
CALL f_disp_cksaft ( fdsp )
3564 IF ( disp_term(1:2) ==
'PT')
CALL f_pert_theory ( fdsp )
3568 IF ( hb_term ==
'TPT1_Chap')
CALL f_association( eps_kij, k_kij, fhb )
3572 CALL f_polar ( fdd, fqq, fdq )
3576 IF ( id_term ==
'TBH')
CALL f_ion_dipole_tbh ( fhend )
3580 IF ( ii_term ==
'primMSA')
CALL f_ion_ion_primmsa ( fcc )
3581 IF ( ii_term ==
'nprMSA')
CALL f_ion_ion_nonprimmsa ( fdd, fqq, fdq, fcc )
3585 IF ( lc_term ==
'MSaupe')
CALL f_lc_mayersaupe ( flc )
3588 IF ( lc_term ==
'OVL')
WRITE(*,*)
'OVL NOT INCLUDED YET' 3591 IF ( lc_term ==
'SPT')
WRITE(*,*)
'SPT NOT INCLUDED YET' 3607 fres = fid + fhs + fhc + fdsp + fhb + fdd + fqq + fdq + fcc + flc
3611 END SUBROUTINE f_numerical
3630 INTEGER :: i, j, k, l, m, n
3631 INTEGER :: ass_cnt,max_eval
3633 REAL :: z0, z1, z2, z3
3635 REAL :: zges, zgesdz, zgesd2, zgesd3
3636 REAL :: zhs, zhsdz, zhsd2, zhsd3
3637 REAL :: zhc, zhcdz, zhcd2, zhcd3
3638 REAL,
DIMENSION(nc,nc) :: dgijdz, dgijd2, dgijd3, dgijd4
3639 REAL :: zdsp, zdspdz, zdspd2, zdspd3
3640 REAL :: c1_con, c2_con, c3_con, c4_con, c5_con
3641 REAL :: i2, edi1dz, edi2dz, edi1d2, edi2d2
3642 REAL :: edi1d3, edi2d3, edi1d4, edi2d4
3643 REAL :: fdspdz,fdspd2
3644 REAL :: zhb, zhbdz, zhbd2, zhbd3
3645 REAL,
DIMENSION(nc,nc,nsite,nsite) :: delta, dq_dz, dq_d2, dq_d3, dq_d4
3646 REAL,
DIMENSION(nc,nsite) :: mx_itr, dmx_dz, ndmxdz, dmx_d2, ndmxd2
3647 REAL,
DIMENSION(nc,nsite) :: dmx_d3, ndmxd3, dmx_d4, ndmxd4
3648 REAL :: err_sum, sum0, sum1, sum2, sum3, sum4, attenu, tol
3649 REAL :: sum_d1, sum_d2, sum_d3, sum_d4
3650 REAL :: zdd, zddz, zddz2, zddz3
3651 REAL :: zqq, zqqz, zqqz2, zqqz3
3652 REAL :: zdq, zdqz, zdqz2, zdqz3
3665 m_mean = z0t/(pi/6.0)
3685 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
3686 dgijdz(i,j)= 1.0/zms/zms + 3.0*dij_ab(i,j)*z2*(1.0+z3)/z3/zms**3 &
3687 + (dij_ab(i,j)*z2/zms/zms)**2 *(4.0+2.0*z3)/z3
3688 dgijd2(i,j) = 2.0/zms**3 &
3689 + 6.0*dij_ab(i,j)*z2/z3/zms**4 *(2.0+z3) &
3690 + (2.0*dij_ab(i,j)*z2/z3)**2 /zms**5 *(1.0+4.0*z3+z3*z3)
3691 dgijd3(i,j) = 6.0/zms**4 &
3692 + 18.0*dij_ab(i,j)*z2/z3/zms**5 *(3.0+z3) &
3693 + 12.0*(dij_ab(i,j)*z2/z3/zms**3 )**2 *(3.0+6.0*z3+z3*z3)
3694 dgijd4(i,j) = 24.0/zms**5 &
3695 + 72.0*dij_ab(i,j)*z2/z3/zms**6 *(4.0+z3) &
3696 + 48.0*(dij_ab(i,j)*z2/z3)**2 /zms**7 *(6.0+8.0*z3+z3*z3)
3704 zhs = m_mean* ( z3/zms + 3.0*z1*z2/z0/zms/zms + z2**3 /z0*(3.0-z3)/zms**3 )
3705 zhsdz = m_mean*( 1.0/zms/zms + 3.0*z1*z2/z0/z3*(1.0+z3)/zms**3 &
3706 + 6.0*z2**3 /z0/z3/zms**4 )
3707 zhsd2 = m_mean*( 2.0/zms**3 + 6.0*z1*z2/z0/z3*(2.0+z3)/zms**4 &
3708 + 6.0*z2**3 /z0/z3/z3*(1.0+3.0*z3)/zms**5 )
3709 zhsd3 = m_mean*( 6.0/zms**4 + 18.0*z1*z2/z0/z3*(3.0+z3)/zms**5 &
3710 + 24.0*z2**3 /z0/z3/z3*(2.0+3.0*z3)/zms**6 )
3721 zhc = zhc + x(i)*(1.0-mseg(i))*eta/gij(i,i)* dgijdz(i,i)
3722 zhcdz = zhcdz + x(i)*(1.0-mseg(i)) *(-eta*(dgijdz(i,i)/gij(i,i))**2 &
3723 + dgijdz(i,i)/gij(i,i) + eta/gij(i,i)*dgijd2(i,i))
3724 zhcd2 = zhcd2 + x(i)*(1.0-mseg(i)) &
3725 *( 2.0*eta*(dgijdz(i,i)/gij(i,i))**3 &
3726 -2.0*(dgijdz(i,i)/gij(i,i))**2 &
3727 -3.0*eta/gij(i,i)**2 *dgijdz(i,i)*dgijd2(i,i) &
3728 +2.0/gij(i,i)*dgijd2(i,i) +eta/gij(i,i)*dgijd3(i,i) )
3729 zhcd3 = zhcd3 + x(i)*(1.0-mseg(i)) *( 6.0*(dgijdz(i,i)/gij(i,i))**3 &
3730 -6.0*eta*(dgijdz(i,i)/gij(i,i))**4 &
3731 +12.0*eta/gij(i,i)**3 *dgijdz(i,i)**2 *dgijd2(i,i) &
3732 -9.0/gij(i,i)**2 *dgijdz(i,i)*dgijd2(i,i) +3.0/gij(i,i)*dgijd3(i,i) &
3733 -3.0*eta*(dgijd2(i,i)/gij(i,i))**2 &
3734 -4.0*eta/gij(i,i)**2 *dgijdz(i,i)*dgijd3(i,i) &
3735 +eta/gij(i,i)*dgijd4(i,i) )
3754 i2 = i2 + bpar(m)*z3**
REAL(m)
3755 edi1dz= edi1dz+apar(m)*
REAL(m+1)*z3**
REAL(m)
3756 edi2dz= edi2dz+bpar(m)*
REAL(m+1)*z3**
REAL(m)
3757 edi1d2= edi1d2+apar(m)*
REAL((m+1)*m)*z3**
REAL(m-1)
3758 edi2d2= edi2d2+bpar(m)*
REAL((m+1)*m)*z3**
REAL(m-1)
3759 edi1d3= edi1d3+apar(m)*
REAL((m+1)*m*(m-1))*z3**
REAL(m-2)
3760 edi2d3= edi2d3+bpar(m)*
REAL((m+1)*m*(m-1))*z3**
REAL(m-2)
3761 edi1d4= edi1d4+apar(m)*
REAL((m+1)*m*(m-1)*(m-2))*z3**
REAL(m-3)
3762 edi2d4= edi2d4+bpar(m)*
REAL((m+1)*m*(m-1)*(m-2))*z3**
REAL(m-3)
3765 c1_con= 1.0/ ( 1.0 + m_mean*(8.0*eta-2.0*eta**2 )/zms**4 &
3766 + (1.0 - m_mean)*(20.0*eta-27.0*eta**2 &
3767 + 12.0*eta**3 -2.0*eta**4 ) /(zms*(2.0-eta))**2 )
3768 c2_con= - c1_con*c1_con &
3769 *(m_mean*(-4.0*eta**2 +20.0*eta+8.0)/zms**5 + (1.0 - m_mean) &
3770 *(2.0*eta**3 +12.0*eta**2 -48.0*eta+40.0) &
3771 /(zms*(2.0-eta))**3 )
3772 c3_con= 2.0 * c2_con*c2_con/c1_con - c1_con*c1_con &
3773 *( m_mean*(-12.0*eta**2 +72.0*eta+60.0)/zms**6 &
3775 *(-6.0*eta**4 -48.0*eta**3 +288.0*eta**2 &
3776 -480.0*eta+264.0) /(zms*(2.0-eta))**4 )
3777 c4_con= 6.0*c2_con*c3_con/c1_con -6.0*c2_con**3 /c1_con**2 &
3779 *( m_mean*(-48.0*eta**2 +336.0*eta+432.0)/zms**7 &
3781 *(24.0*eta**5 +240.0*eta**4 -1920.0*eta**3 &
3782 +4800.0*eta**2 -5280.0*eta+2208.0) /(zms*(2.0-eta))**5 )
3783 c5_con= 6.0*c3_con**2 /c1_con - 36.0*c2_con**2 /c1_con**2 *c3_con &
3784 + 8.0*c2_con/c1_con*c4_con+24.0*c2_con**4 /c1_con**3 &
3786 *( m_mean*(-240.0*eta**2 +1920.0*eta+3360.0)/zms**8 &
3788 *(-120.0*eta**6 -1440.0*eta**5 +14400.0*eta**4 &
3789 -48000.0*eta**3 +79200.0*eta**2 -66240.0*eta+22560.0) &
3790 /(zms*(2.0-eta))**6 )
3792 zdsp = - 2.0*pi*rho*edi1dz*order1 &
3793 - pi*rho*order2*m_mean*(c2_con*i2*eta + c1_con*edi2dz)
3794 zdspdz= zdsp/eta - 2.0*pi*rho*edi1d2*order1 &
3795 - pi*rho*order2*m_mean*(c3_con*i2*eta &
3796 + 2.0*c2_con*edi2dz + c1_con*edi2d2)
3797 zdspd2= -2.0*zdsp/eta/eta +2.0*zdspdz/eta &
3798 - 2.0*pi*rho*edi1d3*order1 - pi*rho*order2*m_mean*(c4_con*i2*eta &
3799 + 3.0*c3_con*edi2dz +3.0*c2_con*edi2d2 +c1_con*edi2d3)
3800 zdspd3= 6.0*zdsp/eta**3 -6.0*zdspdz/eta/eta &
3801 + 3.0*zdspd2/eta - 2.0*pi*rho*edi1d4*order1 &
3802 - pi*rho*order2*m_mean*(c5_con*i2*eta &
3803 + 4.0*c4_con*edi2dz +6.0*c3_con*edi2d2 &
3804 + 4.0*c2_con*edi2d3 + c1_con*edi2d4)
3816 fdspdz = fdspdz + m_mean/tau * dnm(n,m) * (um/t)**
REAL(n) *
REAL(m)*(eta/tau)**
REAL(m-1)
3819 fdspd2= fdspd2 + m_mean/tau * dnm(n,m)*(um/t)**
REAL(n) *
REAL(m)*
REAL(m-1) &
3820 * (eta/tau)**
REAL(m-2) * 1.0/tau
3824 zdspdz = (2.0*fdspdz + eta*fdspd2) - zdsp/z3
3839 IF (nhb_typ(i) /= 0) assoc = .true.
3844 DO i = 1, nhb_typ(j)
3846 DO l = 1, nhb_typ(k)
3847 delta(j,k,i,l) = gij(j,k) * ass_d(j,k,i,l)
3848 dq_dz(j,k,i,l) = dgijdz(j,k) * ass_d(j,k,i,l)
3849 dq_d2(j,k,i,l) = dgijd2(j,k) * ass_d(j,k,i,l)
3850 dq_d3(j,k,i,l) = dgijd3(j,k) * ass_d(j,k,i,l)
3851 dq_d4(j,k,i,l) = dgijd4(j,k) * ass_d(j,k,i,l)
3860 IF ( eta < 0.2 ) tol = 1.d-12
3861 IF ( eta < 0.01 ) tol = 1.d-13
3862 IF ( eta < 1.e-6) tol = 1.d-15
3867 DO j = 1, nhb_typ(i)
3879 DO WHILE ( err_sum > tol .AND. ass_cnt <= max_eval)
3880 ass_cnt = ass_cnt + 1
3882 DO j = 1, nhb_typ(i)
3889 DO l = 1, nhb_typ(k)
3890 sum0 =sum0 +x(k)*nhb_no(k,l)* mx(k,l)* delta(i,k,j,l)
3891 sum1 =sum1 +x(k)*nhb_no(k,l)*( mx(k,l)* dq_dz(i,k,j,l) &
3892 + dmx_dz(k,l)* delta(i,k,j,l))
3893 sum2 =sum2 +x(k)*nhb_no(k,l)*( mx(k,l)* dq_d2(i,k,j,l) &
3894 + 2.0*dmx_dz(k,l)* dq_dz(i,k,j,l) &
3895 + dmx_d2(k,l)* delta(i,k,j,l))
3896 sum3 =sum3 +x(k)*nhb_no(k,l)*( mx(k,l)* dq_d3(i,k,j,l) &
3897 + 3.0*dmx_dz(k,l)* dq_d2(i,k,j,l) &
3898 + 3.0*dmx_d2(k,l)* dq_dz(i,k,j,l) &
3899 + dmx_d3(k,l)* delta(i,k,j,l))
3900 sum4 =sum4 + x(k)*nhb_no(k,l)*( mx(k,l)* dq_d4(i,k,j,l) &
3901 + 4.0*dmx_dz(k,l)* dq_d3(i,k,j,l) &
3902 + 6.0*dmx_d2(k,l)* dq_d2(i,k,j,l) &
3903 + 4.0*dmx_d3(k,l)* dq_dz(i,k,j,l) &
3904 + dmx_d4(k,l)* delta(i,k,j,l))
3907 mx_itr(i,j)= 1.0 / (1.0 + sum0 * rho)
3908 ndmxdz(i,j)= -(mx_itr(i,j)*mx_itr(i,j))* (sum0/z3t +sum1*rho)
3909 ndmxd2(i,j)= + 2.0/mx_itr(i,j)*ndmxdz(i,j)*ndmxdz(i,j) &
3910 - (mx_itr(i,j)*mx_itr(i,j)) * (2.0/z3t*sum1 + rho*sum2)
3911 ndmxd3(i,j)= - 6.0/mx_itr(i,j)**2 *ndmxdz(i,j)**3 &
3912 + 6.0/mx_itr(i,j)*ndmxdz(i,j)*ndmxd2(i,j) - mx_itr(i,j)*mx_itr(i,j) &
3913 * (3.0/z3t*sum2 + rho*sum3)
3914 ndmxd4(i,j)= 24.0/mx_itr(i,j)**3 *ndmxdz(i,j)**4 &
3915 -36.0/mx_itr(i,j)**2 *ndmxdz(i,j)**2 *ndmxd2(i,j) &
3916 +6.0/mx_itr(i,j)*ndmxd2(i,j)**2 &
3917 +8.0/mx_itr(i,j)*ndmxdz(i,j)*ndmxd3(i,j) - mx_itr(i,j)**2 &
3918 *(4.0/z3t*sum3 + rho*sum4)
3924 DO j = 1, nhb_typ(i)
3925 err_sum = err_sum + abs(mx_itr(i,j) - mx(i,j)) &
3926 + abs(ndmxdz(i,j) - dmx_dz(i,j)) + abs(ndmxd2(i,j) - dmx_d2(i,j))
3927 mx(i,j) = mx_itr(i,j)*attenu + mx(i,j) * (1.0-attenu)
3928 dmx_dz(i,j) = ndmxdz(i,j)*attenu + dmx_dz(i,j) * (1.0-attenu)
3929 dmx_d2(i,j) = ndmxd2(i,j)*attenu + dmx_d2(i,j) * (1.0-attenu)
3930 dmx_d3(i,j) = ndmxd3(i,j)*attenu + dmx_d3(i,j) * (1.0-attenu)
3931 dmx_d4(i,j) = ndmxd4(i,j)*attenu + dmx_d4(i,j) * (1.0-attenu)
3936 IF ( ass_cnt >= max_eval .AND. err_sum > sqrt(tol) )
THEN 3937 WRITE (*,
'(a,2G15.7)')
'P_EOS: Max_eval violated (mx) Err_Sum= ',err_sum,tol
3948 DO j = 1, nhb_typ(i)
3949 sum_d1= sum_d1 +nhb_no(i,j)* dmx_dz(i,j)*(1.0/mx(i,j)-0.5)
3950 sum_d2= sum_d2 +nhb_no(i,j)*(dmx_d2(i,j)*(1.0/mx(i,j)-0.5) &
3951 -(dmx_dz(i,j)/mx(i,j))**2 )
3952 sum_d3= sum_d3 +nhb_no(i,j)*(dmx_d3(i,j)*(1.0/mx(i,j)-0.5) &
3953 -3.0/mx(i,j)**2 *dmx_dz(i,j)*dmx_d2(i,j) + 2.0*(dmx_dz(i,j)/mx(i,j))**3 )
3954 sum_d4= sum_d4 +nhb_no(i,j)*(dmx_d4(i,j)*(1.0/mx(i,j)-0.5) &
3955 -4.0/mx(i,j)**2 *dmx_dz(i,j)*dmx_d3(i,j) &
3956 + 12.0/mx(i,j)**3 *dmx_dz(i,j)**2 *dmx_d2(i,j) &
3957 - 3.0/mx(i,j)**2 *dmx_d2(i,j)**2 - 6.0*(dmx_dz(i,j)/mx(i,j))**4 )
3959 zhb = zhb + x(i) * eta * sum_d1
3960 zhbdz = zhbdz + x(i) * eta * sum_d2
3961 zhbd2 = zhbd2 + x(i) * eta * sum_d3
3962 zhbd3 = zhbd3 + x(i) * eta * sum_d4
3964 zhbdz = zhbdz + zhb/eta
3965 zhbd2 = zhbd2 + 2.0/eta*zhbdz-2.0/eta**2 *zhb
3966 zhbd3 = zhbd3 - 6.0/eta**2 *zhbdz+3.0/eta*zhbd2 + 6.0/eta**3 *zhb
3974 CALL p_polar ( zdd, zddz, zddz2, zddz3, zqq, zqqz, zqqz2, zqqz3, zdq, zdqz, zdqz2, zdqz3 )
3981 zges = 1.0 + zhs + zhc + zdsp + zhb + zdd + zqq + zdq
3982 zgesdz = zhsdz + zhcdz + zdspdz + zhbdz + zddz + zqqz + zdqz
3983 zgesd2 = zhsd2 + zhcd2 + zdspd2 + zhbd2 + zddz2 +zqqz2+zdqz2
3984 zgesd3 = zhsd3 + zhcd3 + zdspd3 + zhbd3 + zddz3 +zqqz3+zdqz3
3986 pges = zges *rho *(kbol*t)/1.d-30
3987 pgesdz = ( zgesdz*rho + zges*rho/z3 )*(kbol*t)/1.d-30
3988 pgesd2 = ( zgesd2*rho + 2.0*zgesdz*rho/z3 )*(kbol*t)/1.d-30
3989 pgesd3 = ( zgesd3*rho + 3.0*zgesd2*rho/z3 )*(kbol*t)/1.d-30
3991 END SUBROUTINE p_eos
4000 SUBROUTINE phi_polar ( k, z3_rk, fdd_rk, fqq_rk, fdq_rk )
4002 USE eos_variables, ONLY: ncomp, parame, dd_term, qq_term, dq_term
4006 INTEGER,
INTENT(IN) :: k
4007 REAL,
INTENT(IN) :: z3_rk
4008 REAL,
INTENT(OUT) :: fdd_rk, fqq_rk, fdq_rk
4012 INTEGER :: quadrupole
4013 INTEGER :: dipole_quad
4023 IF ( sum( parame(1:ncomp,6) ) /= 0.0 ) dipole = 1
4024 IF ( sum( parame(1:ncomp,7) ) /= 0.0 ) quadrupole = 1
4025 IF ( dipole == 1 .AND. quadrupole == 1 ) dipole_quad = 1
4030 IF (dipole == 1)
THEN 4032 IF (dd_term ==
'GV')
CALL phi_dd_gross_vrabec( k, z3_rk, fdd_rk )
4035 IF (dd_term /=
'GV' .AND. dd_term /=
'SF')
write (*,*)
'specify dipole term !' 4042 IF (quadrupole == 1)
THEN 4045 IF (qq_term ==
'JG')
CALL phi_qq_gross( k, z3_rk, fqq_rk )
4047 IF (qq_term /=
'JG' .AND. qq_term /=
'SF')
write (*,*)
'specify quadrupole term !' 4054 IF (dipole_quad == 1)
THEN 4056 IF (dq_term ==
'VG')
CALL phi_dq_vrabec_gross( k, z3_rk, fdq_rk )
4058 IF (dq_term /=
'VG' )
write (*,*)
'specify DQ-cross term !' 4062 END SUBROUTINE phi_polar
4069 SUBROUTINE phi_dd_gross_vrabec( k, z3_rk, fdd_rk )
4072 USE eos_variables, ONLY: nc, ncomp, uij, parame, mseg, sig_ij, rho, eta, x, t
4077 INTEGER,
INTENT(IN) :: k
4078 REAL,
INTENT(IN) :: z3_rk
4079 REAL,
INTENT(IN OUT) :: fdd_rk
4082 INTEGER :: i, j, l, m
4084 REAL :: factor2, factor3, z3
4085 REAL :: xijfa, xijkfa, xijfa_x, xijkf_x, eij
4086 REAL :: fdd2, fdd3, fdd2x, fdd3x
4087 REAL,
DIMENSION(nc) :: my2dd
4088 REAL,
DIMENSION(nc,nc) :: idd2, idd4, idd2x, idd4x
4089 REAL,
DIMENSION(nc,nc,nc) :: idd3, idd3x
4096 IF ( uij(i,i) == 0.0 )
write (*,*)
'PHI_DD_GROSS_VRABEC: do not use dimensionless units' 4097 IF ( uij(i,i) == 0.0 ) stop
4098 my2dd(i) = (parame(i,6))**2 *1.e-49 / (uij(i,i)*kbol*mseg(i)*sig_ij(i,i)**3 *1.e-30)
4107 IF (parame(i,6) /= 0.0 .AND. parame(j,6) /= 0.0)
THEN 4109 idd2(i,j) =idd2(i,j) + ddp2(i,j,m)*z3**m
4110 idd4(i,j) =idd4(i,j) + ddp4(i,j,m)*z3**m
4111 idd2x(i,j) =idd2x(i,j)+ ddp2(i,j,m)*
REAL(m)*z3**(m-1)*z3_rk
4112 idd4x(i,j) =idd4x(i,j)+ ddp4(i,j,m)*
REAL(m)*z3**(m-1)*z3_rk
4117 IF (parame(l,6) /= 0.0)
THEN 4119 idd3(i,j,l) =idd3(i,j,l) +ddp3(i,j,l,m)*z3**m
4120 idd3x(i,j,l)=idd3x(i,j,l)+ddp3(i,j,l,m)*
REAL(m)*z3**(m-1)*z3_rk
4129 factor3= -4.0/3.0*pi**2
4136 xijfa_x = 2.0*x(i)*rho*uij(i,i)*my2dd(i)*sig_ij(i,i)**3 /t &
4137 *uij(k,k)*my2dd(k)*sig_ij(k,k)**3 /t/sig_ij(i,k)**3
4138 eij = (parame(i,3)*parame(k,3))**0.5
4139 fdd2x = fdd2x + factor2*xijfa_x*( idd2(i,k) + eij/t*idd4(i,k) )
4141 IF (parame(i,6) /= 0.0 .AND. parame(j,6) /= 0.0)
THEN 4142 xijfa =x(i)*rho*uij(i,i)*my2dd(i)*sig_ij(i,i)**3 /t &
4143 *x(j)*rho*uij(j,j)*my2dd(j)*sig_ij(j,j)**3 /t/sig_ij(i,j)**3
4144 eij = (parame(i,3)*parame(j,3))**0.5
4145 fdd2= fdd2 +factor2*xijfa*(idd2(i,j) +eij/t*idd4(i,j) )
4146 fdd2x =fdd2x +factor2*xijfa*(idd2x(i,j)+eij/t*idd4x(i,j))
4148 xijkf_x=x(i)*rho*uij(i,i)*my2dd(i)*sig_ij(i,i)**3 /t/sig_ij(i,j) &
4149 *x(j)*rho*uij(j,j)*my2dd(j)*sig_ij(j,j)**3 /t/sig_ij(i,k) &
4150 *3.0* uij(k,k)*my2dd(k)*sig_ij(k,k)**3 /t/sig_ij(j,k)
4151 fdd3x=fdd3x+factor3*xijkf_x*idd3(i,j,k)
4153 IF (parame(l,6) /= 0.0)
THEN 4154 xijkfa= x(i)*rho*uij(i,i)/t*my2dd(i)*sig_ij(i,i)**3 &
4155 *x(j)*rho*uij(j,j)/t*my2dd(j)*sig_ij(j,j)**3 &
4156 *x(l)*rho*uij(l,l)/t*my2dd(l)*sig_ij(l,l)**3 &
4157 /sig_ij(i,j)/sig_ij(i,l)/sig_ij(j,l)
4158 fdd3 =fdd3 + factor3 * xijkfa *idd3(i,j,l)
4159 fdd3x =fdd3x + factor3 * xijkfa *idd3x(i,j,l)
4166 IF (fdd2 < -1.e-50 .AND. fdd3 /= 0.0 .AND. fdd2x /= 0.0 .AND. fdd3x /= 0.0)
THEN 4168 fdd_rk = fdd2* (fdd2*fdd2x - 2.0*fdd3*fdd2x+fdd2*fdd3x) / (fdd2-fdd3)**2
4172 END SUBROUTINE phi_dd_gross_vrabec
4180 SUBROUTINE phi_qq_gross( k, z3_rk, fqq_rk )
4183 USE eos_variables, ONLY: nc, ncomp, uij, parame, mseg, sig_ij, rho, eta, x, t
4188 INTEGER,
INTENT(IN) :: k
4189 REAL,
INTENT(IN) :: z3_rk
4190 REAL,
INTENT(IN OUT) :: fqq_rk
4193 INTEGER :: i, j, l, m
4195 REAL :: factor2, factor3, z3
4196 REAL :: xijfa, xijkfa, xijfa_x, xijkf_x, eij
4197 REAL :: fqq2, fqq3, fqq2x, fqq3x
4198 REAL,
DIMENSION(nc) :: qq2
4199 REAL,
DIMENSION(nc,nc) :: iqq2, iqq4, iqq2x, iqq4x
4200 REAL,
DIMENSION(nc,nc,nc) :: iqq3, iqq3x
4207 IF ( uij(i,i) == 0.0 )
write (*,*)
'PHI_QQ_GROSS: do not use dimensionless units' 4208 qq2(i) = (parame(i,7))**2 *1.e-69 / (uij(i,i)*kbol*mseg(i)*sig_ij(i,i)**5 *1.e-50)
4217 IF (parame(i,7) /= 0.0 .AND. parame(j,7) /= 0.0)
THEN 4219 iqq2(i,j) = iqq2(i,j) + qqp2(i,j,m) * z3**m
4220 iqq4(i,j) = iqq4(i,j) + qqp4(i,j,m) * z3**m
4221 iqq2x(i,j) = iqq2x(i,j) + qqp2(i,j,m) *
REAL(m)*z3**(m-1)*z3_rk
4222 iqq4x(i,j) = iqq4x(i,j) + qqp4(i,j,m) *
REAL(m)*z3**(m-1)*z3_rk
4227 IF (parame(l,7) /= 0.0)
THEN 4229 iqq3(i,j,l) = iqq3(i,j,l) + qqp3(i,j,l,m)*z3**m
4230 iqq3x(i,j,l) = iqq3x(i,j,l) + qqp3(i,j,l,m)*
REAL(m) *z3**(m-1)*z3_rk
4238 factor2= -9.0/16.0*pi
4239 factor3= 9.0/16.0*pi**2
4246 xijfa_x = 2.0*x(i)*rho*uij(i,i)*qq2(i)*sig_ij(i,i)**5 /t &
4247 *uij(k,k)*qq2(k)*sig_ij(k,k)**5 /t/sig_ij(i,k)**7.0
4248 eij = (parame(i,3)*parame(k,3))**0.5
4249 fqq2x =fqq2x +factor2*xijfa_x*(iqq2(i,k)+eij/t*iqq4(i,k))
4251 IF (parame(i,7) /= 0.0 .AND. parame(j,7) /= 0.0)
THEN 4252 xijfa =x(i)*rho*uij(i,i)*qq2(i)*sig_ij(i,i)**5 /t &
4253 *x(j)*rho*uij(j,j)*qq2(j)*sig_ij(j,j)**5 /t/sig_ij(i,j)**7.0
4254 eij = (parame(i,3)*parame(j,3))**0.5
4255 fqq2= fqq2 +factor2*xijfa*(iqq2(i,j) +eij/t*iqq4(i,j) )
4256 fqq2x =fqq2x +factor2*xijfa*(iqq2x(i,j)+eij/t*iqq4x(i,j))
4258 xijkf_x=x(i)*rho*uij(i,i)*qq2(i)*sig_ij(i,i)**5 /t/sig_ij(i,j)**3 &
4259 *x(j)*rho*uij(j,j)*qq2(j)*sig_ij(j,j)**5 /t/sig_ij(i,k)**3 &
4260 *3.0* uij(k,k)*qq2(k)*sig_ij(k,k)**5 /t/sig_ij(j,k)**3
4261 fqq3x = fqq3x + factor3*xijkf_x*iqq3(i,j,k)
4263 IF (parame(l,7) /= 0.0)
THEN 4264 xijkfa=x(i)*rho*uij(i,i)*qq2(i)*sig_ij(i,i)**5 /t/sig_ij(i,j)**3 &
4265 *x(j)*rho*uij(j,j)*qq2(j)*sig_ij(j,j)**5 /t/sig_ij(i,l)**3 &
4266 *x(l)*rho*uij(l,l)*qq2(l)*sig_ij(l,l)**5 /t/sig_ij(j,l)**3
4267 fqq3 =fqq3 + factor3 * xijkfa *iqq3(i,j,l)
4268 fqq3x =fqq3x + factor3 * xijkfa *iqq3x(i,j,l)
4275 IF (fqq2 < -1.e-50 .AND. fqq3 /= 0.0 .AND. fqq2x /= 0.0 .AND. fqq3x /= 0.0)
THEN 4276 fqq_rk = fqq2* (fqq2*fqq2x - 2.0*fqq3*fqq2x+fqq2*fqq3x) / (fqq2-fqq3)**2
4279 END SUBROUTINE phi_qq_gross
4285 SUBROUTINE phi_dq_vrabec_gross( k, z3_rk, fdq_rk )
4288 USE eos_variables, ONLY: nc, ncomp, uij, parame, mseg, sig_ij, rho, eta, x, t
4293 INTEGER,
INTENT(IN) :: k
4294 REAL,
INTENT(IN) :: z3_rk
4295 REAL,
INTENT(IN OUT) :: fdq_rk
4298 INTEGER :: i, j, l, m
4300 REAL :: factor2, factor3, z3
4301 REAL :: xijfa, xijkfa, xijfa_x, xijkf_x, eij
4302 REAL :: fdq2, fdq3, fdq2x, fdq3x
4303 REAL,
DIMENSION(nc) :: my2dd, myfac, qq2, q_fac
4304 REAL,
DIMENSION(nc,nc) :: idq2, idq4, idq2x, idq4x
4305 REAL,
DIMENSION(nc,nc,nc) :: idq3, idq3x
4311 my2dd(i) = (parame(i,6))**2 *1.e-49 / (uij(i,i)*kbol*mseg(i)*sig_ij(i,i)**3 *1.e-30)
4312 myfac(i) = parame(i,3)/t*parame(i,2)**4 *my2dd(i)
4313 qq2(i) = (parame(i,7))**2 *1.e-69 / (uij(i,i)*kbol*mseg(i)*sig_ij(i,i)**5 *1.e-50)
4314 q_fac(i) = parame(i,3)/t*parame(i,2)**4 *qq2(i)
4324 idq2(i,j) = idq2(i,j) + dqp2(i,j,m)*z3**m
4325 idq4(i,j) = idq4(i,j) + dqp4(i,j,m)*z3**m
4326 idq2x(i,j) = idq2x(i,j) + dqp2(i,j,m)*
REAL(m)*z3**(m-1) *z3_rk
4327 idq4x(i,j) = idq4x(i,j) + dqp4(i,j,m)*
REAL(m)*z3**(m-1) *z3_rk
4333 idq3(i,j,l) =idq3(i,j,l) +dqp3(i,j,l,m)*z3**m
4334 idq3x(i,j,l)=idq3x(i,j,l)+dqp3(i,j,l,m)*
REAL(m)*z3**(m-1)*z3_rk
4340 factor2= -9.0/4.0*pi
4348 xijfa_x = x(i)*rho*( myfac(i)*q_fac(k) + myfac(k)*q_fac(i) ) / sig_ij(i,k)**5
4349 eij = (parame(i,3)*parame(k,3))**0.5
4350 fdq2x =fdq2x +factor2*xijfa_x*(idq2(i,k)+eij/t*idq4(i,k))
4352 xijfa =x(i)*rho*myfac(i) * x(j)*rho*q_fac(j) /sig_ij(i,j)**5
4353 eij = (parame(i,3)*parame(j,3))**0.5
4354 fdq2= fdq2 +factor2*xijfa*(idq2(i,j) +eij/t*idq4(i,j) )
4355 fdq2x =fdq2x +factor2*xijfa*(idq2x(i,j) +eij/t*idq4x(i,j))
4357 xijkf_x=x(i)*rho*x(j)*rho/(sig_ij(i,j)*sig_ij(i,k)*sig_ij(j,k))**2 &
4358 *( myfac(i)*q_fac(j)*myfac(k) + myfac(i)*q_fac(k)*myfac(j) &
4359 + myfac(k)*q_fac(i)*myfac(j) +myfac(i)*q_fac(j)*q_fac(k)*1.1937350 &
4360 +myfac(i)*q_fac(k)*q_fac(j)*1.193735 &
4361 +myfac(k)*q_fac(i)*q_fac(j)*1.193735 )
4362 fdq3x = fdq3x + factor3*xijkf_x*idq3(i,j,k)
4364 xijkfa=x(i)*rho*x(j)*rho*x(l)*rho/(sig_ij(i,j)*sig_ij(i,l)*sig_ij(j,l))**2 &
4365 *( myfac(i)*q_fac(j)*myfac(l) &
4366 +myfac(i)*q_fac(j)*q_fac(l)*1.193735 )
4367 fdq3 =fdq3 + factor3 * xijkfa *idq3(i,j,l)
4368 fdq3x =fdq3x + factor3 * xijkfa *idq3x(i,j,l)
4373 IF (fdq2 < -1.e-50 .AND. fdq3 /= 0.0 .AND. fdq2x /= 0.0 .AND. fdq3x /= 0.0)
THEN 4375 fdq_rk = fdq2* (fdq2*fdq2x - 2.0*fdq3*fdq2x+fdq2*fdq3x) / (fdq2-fdq3)**2
4379 END SUBROUTINE phi_dq_vrabec_gross
4387 SUBROUTINE p_numerical
4393 REAL :: dzetdv, eta_0, dist, fact
4394 REAL :: fres1, fres2, fres3, fres4, fres5
4395 REAL :: df_dr, df_drdr, pideal, dpiddz
4396 REAL :: tfr_1, tfr_2, tfr_3, tfr_4, tfr_5
4402 ELSE IF (eta <= 0.1 .AND. eta > 0.01)
THEN 4407 dist = eta*3.d-3 *fact
4414 eta = eta_0 - 2.0*dist
4426 eta = eta_0 + 2.0*dist
4452 df_dr = (-fres4+8.0*fres3-8.0*fres2+fres1) / (12.0*dist)
4453 df_drdr = (-fres4+16.0*fres3-3.d1*fres5+16.0*fres2-fres1) &
4459 pges = (-fres4+8.0*fres3-8.0*fres2+fres1) &
4460 /(12.0*dist) *dzetdv*(kbol*t)/1.e-30
4462 dpiddz = 1.0/z3t*(kbol*t)/1.e-30
4463 pgesdz = (-fres4+16.0*fres3-3.d1*fres5+16.0*fres2-fres1) &
4464 /(12.0*(dist**2 ))* dzetdv*(kbol*t)/1.e-30 &
4465 + (-fres4+8.0*fres3-8.0*fres2+fres1) /(12.0*dist) * 2.0 *rho &
4466 *(kbol*t)/1.e-30 + dpiddz
4468 pgesd2 = (fres4-2.0*fres3+2.0*fres2-fres1) /(2.0*dist**3 ) &
4469 * dzetdv*(kbol*t)/1.e-30 &
4470 + (-fres4+16.0*fres3-3.d1*fres5+16.0*fres2-fres1) /(12.0*(dist**2 )) &
4471 * 4.0 *rho *(kbol*t)/1.e-30 + (-fres4+8.0*fres3-8.0*fres2+fres1) &
4472 /(12.0*dist) * 2.0 /z3t *(kbol*t)/1.e-30
4473 pgesd3 = (fres4-4.0*fres3+6.0*fres5-4.0*fres2+fres1) /(dist**4 ) &
4474 * dzetdv*(kbol*t)/1.e-30 + (fres4-2.0*fres3+2.0*fres2-fres1) &
4475 /(2.0*dist**3 ) * 6.0 *rho *(kbol*t)/1.e-30 &
4476 + (-fres4+16.0*fres3-3.d1*fres5+16.0*fres2-fres1) &
4477 /(12.0*dist**2 )* 6.0 /z3t *(kbol*t)/1.e-30
4480 pideal = rho * (kbol*t) / 1.e-30
4483 pges = pideal + pges
4485 END SUBROUTINE p_numerical
4492 SUBROUTINE h_numerical
4499 REAL :: dist, fact, rho_0
4500 REAL :: fres1,fres2,fres3,fres4,fres5
4501 REAL :: f_1, f_2, f_3, f_4
4502 REAL :: cv_res, t_tmp, zges
4503 REAL :: f_dt, f_dtdt, f_dr, f_drdr, f_drdt
4507 CALL perturbation_parameter
4512 dist = t * 100.e-5 * fact
4516 t = t_tmp - 2.0*dist
4517 CALL perturbation_parameter
4522 CALL perturbation_parameter
4527 CALL perturbation_parameter
4531 t = t_tmp + 2.0*dist
4532 CALL perturbation_parameter
4537 CALL perturbation_parameter
4543 zges = (p * 1.e-30)/(kbol*t*rho_0)
4546 f_dt = (-fres4+8.0*fres3-8.0*fres2+fres1)/(12.0*dist)
4547 f_dtdt = (-fres4+16.0*fres3-3.d1*fres5+16.0*fres2-fres1) /(12.0*(dist**2 ))
4549 s_res = (- f_dt -fres/t)*rgas*t + rgas * log(zges)
4550 h_res = ( - t*f_dt + zges-1.0 ) * rgas*t
4551 cv_res = -( t*f_dtdt + 2.0*f_dt ) * rgas*t
4555 t = t_tmp - 2.0*dist
4556 CALL perturbation_parameter
4559 f_1 = pges/(eta*rho_0*(kbol*t)/1.e-30)
4562 CALL perturbation_parameter
4565 f_2 = pges/(eta*rho_0*(kbol*t)/1.e-30)
4568 CALL perturbation_parameter
4571 f_3 = pges/(eta*rho_0*(kbol*t)/1.e-30)
4573 t = t_tmp + 2.0*dist
4574 CALL perturbation_parameter
4577 f_4 = pges/(eta*rho_0*(kbol*t)/1.e-30)
4580 CALL perturbation_parameter
4584 f_dr = pges / (eta*rho_0*(kbol*t)/1.e-30)
4585 f_drdr = pgesdz/ (eta*rho_0*(kbol*t)/1.e-30) - f_dr*2.0/eta - 1.0/eta**2
4587 f_drdt = ( - f_4 + 8.0*f_3 - 8.0*f_2 + f_1 ) / ( 12.0*dist )
4589 cp_res = cv_res - rgas + rgas*( zges + eta*t*f_drdt)**2 / (1.0 + 2.0*eta*f_dr + eta**2 *f_drdr)
4593 END SUBROUTINE h_numerical
4609 SUBROUTINE f_polar ( fdd, fqq, fdq )
4611 USE eos_variables, ONLY: ncomp, parame, dd_term, qq_term, dq_term
4615 REAL,
INTENT(OUT) :: fdd, fqq, fdq
4619 INTEGER :: quadrupole
4620 INTEGER :: dipole_quad
4630 IF ( sum( parame(1:ncomp,6) ) /= 0.0 ) dipole = 1
4631 IF ( sum( parame(1:ncomp,7) ) /= 0.0 ) quadrupole = 1
4632 IF ( dipole == 1 .AND. quadrupole == 1 ) dipole_quad = 1
4637 IF (dipole == 1)
THEN 4639 IF (dd_term ==
'GV')
CALL f_dd_gross_vrabec( fdd )
4641 IF (dd_term /=
'GV' .AND. dd_term /=
'SF')
write (*,*)
'specify dipole term !' 4648 IF (quadrupole == 1)
THEN 4651 IF (qq_term ==
'JG')
CALL f_qq_gross( fqq )
4652 IF (qq_term /=
'JG' .AND. qq_term /=
'SF')
write (*,*)
'specify quadrupole term !' 4659 IF (dipole_quad == 1)
THEN 4661 IF (dq_term ==
'VG')
CALL f_dq_vrabec_gross( fdq )
4662 IF (dq_term /=
'VG' )
write (*,*)
'specify DQ-cross term !' 4666 END SUBROUTINE f_polar
4678 SUBROUTINE pressure_spinodal
4686 REAL :: eta_iteration
4687 REAL :: error, acc_i, delta_eta
4694 eta_iteration = eta_start
4701 DO WHILE ( abs(error) > acc_i .AND. i < max_i )
4706 IF ( num == 0 )
THEN 4708 ELSE IF ( num == 1 )
THEN 4710 ELSE IF ( num == 2 )
THEN 4711 WRITE(*,*)
'CRITICAL RENORM NOT INCLUDED YET' 4714 write (*,*)
'define calculation option (num)' 4719 delta_eta = error/ pgesd2
4720 IF ( delta_eta > 0.02 ) delta_eta = 0.02
4721 IF ( delta_eta < -0.02 ) delta_eta = -0.02
4723 eta_iteration = eta_iteration - delta_eta
4726 IF (eta_iteration > 0.9) eta_iteration = 0.5
4727 IF (eta_iteration <= 0.0) eta_iteration = 1.e-16
4733 END SUBROUTINE pressure_spinodal
4740 SUBROUTINE f_ideal_gas ( fid )
4742 USE eos_variables, ONLY: nc, ncomp, t, x, rho, pi, kbol, nav
4746 REAL,
INTENT(IN OUT) :: fid
4749 REAL,
DIMENSION(nc) :: rhoi
4754 rhoi(i) = x(i) * rho
4759 fid = fid + x(i) * ( log(rhoi(i)) - 1.0 )
4762 END SUBROUTINE f_ideal_gas
4768 SUBROUTINE f_hard_sphere ( m_mean2, fhs )
4774 REAL,
INTENT(IN) :: m_mean2
4775 REAL,
INTENT(IN OUT) :: fhs
4777 REAL :: z0, z1, z2, z3, zms
4787 fhs= m_mean2*( 3.0*z1*z2/zms + z2**3 /z3/zms/zms + (z2**3 /z3/z3-z0)*log(zms) )/z0
4790 END SUBROUTINE f_hard_sphere
4796 SUBROUTINE f_chain_tpt1 ( fhc )
4798 USE eos_variables, ONLY: nc, ncomp, mseg, x, z0t, z1t, z2t, z3t, &
4799 rho, eta, dij_ab, gij
4803 REAL,
INTENT(IN OUT) :: fhc
4806 REAL :: z0, z1, z2, z3, zms
4818 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
4824 fhc = fhc + x(i) *(1.0- mseg(i)) *log(gij(i,i))
4827 END SUBROUTINE f_chain_tpt1
4834 SUBROUTINE f_chain_tpt_d ( fhc )
4836 USE eos_variables, ONLY: nc, ncomp, mseg, x, z0t, z1t, z2t, z3t, rho, eta, &
4837 dhs, mseg, dij_ab, gij
4841 REAL,
INTENT(OUT) :: fhc
4844 REAL,
DIMENSION(nc) :: gij_hd
4845 REAL :: z0, z1, z2, z3, zms
4857 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
4862 gij_hd(i) = 1.0/(2.0*zms) + 3.0*dij_ab(i,i)*z2 / zms**2
4867 IF ( mseg(i) >= 2.0 )
THEN 4868 fhc = fhc - x(i) * ( mseg(i)/2.0 * log( gij(i,i) ) + ( mseg(i)/2.0 - 1.0 ) * log( gij_hd(i)) )
4870 fhc = fhc + x(i) * ( 1.0 - mseg(i) ) * log( gij(i,i) )
4874 END SUBROUTINE f_chain_tpt_d
4881 SUBROUTINE f_chain_hu_liu ( fhc )
4888 REAL,
INTENT(OUT) :: fhc
4890 REAL :: a2, b2, c2, a3, b3, c3
4891 REAL :: a20, b20, c20, a30, b30, c30
4892 REAL :: sum1, sum2, am, bm, cm
4898 sum1 = sum( x(1:ncomp)*(mseg(1:ncomp)-1.0) )
4899 sum2 = sum( x(1:ncomp)/mseg(1:ncomp)*(mseg(1:ncomp)-1.0)*(mseg(1:ncomp)-2.0) )
4907 a20 = - a2 + b2 - 3.0*c2
4908 b20 = - a2 - b2 + c2
4910 a30 = - a3 + b3 - 3.0*c3
4911 b30 = - a3 - b3 + c3
4913 am = (3.0 + a20) * sum1 + a30 * sum2
4914 bm = (1.0 + b20) * sum1 + b30 * sum2
4915 cm = (1.0 + c20) * sum1 + c30 * sum2
4917 fhc = - ( (am*eta - bm) / (2.0*zms) + bm/2.0/zms**2 - cm *log(zms) )
4920 END SUBROUTINE f_chain_hu_liu
4926 SUBROUTINE f_chain_hu_liu_rc ( fhs, fhc )
4933 REAL,
INTENT(IN) :: fhs
4934 REAL,
INTENT(OUT) :: fhc
4936 REAL :: a2, b2, c2, a3, b3, c3
4937 REAL :: para1,para2,para3,para4
4948 para2 = 0.299154629727814
4949 para3 = 1.087271036653154
4950 para4 = -0.708979110326831
4951 a3 = para1 + para2*chir(1) + para3*chir(1)**2 + para4*chir(1)**3
4952 b3 = 3.49695 - (3.49695 + 0.317719074806190)*chir(1)
4953 c3 = 4.83207 - (4.83207 - 3.480163780334421)*chir(1)
4955 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 )
4956 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 )
4957 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 )
4959 fhc = ((3.0 + alh - blh + 3.0*clh)*eta - (1.0 + alh + blh - clh)) / (2.0*(1.0-eta)) + &
4960 (1.0 + alh + blh - clh) / ( 2.0*(1.0-eta)**2 ) + (clh - 1.0)*log(1.0-eta)
4964 END SUBROUTINE f_chain_hu_liu_rc
4973 SUBROUTINE f_disp_pcsaft ( fdsp )
4975 USE eos_variables, ONLY: pi, rho, eta, z0t, apar, bpar, order1, order2
4979 REAL,
INTENT(IN OUT) :: fdsp
4982 REAL :: i1, i2, c1_con, z3, zms, m_mean
4987 m_mean = z0t / ( pi / 6.0 )
4992 i1 = i1 + apar(m) * z3**m
4993 i2 = i2 + bpar(m) * z3**m
4996 c1_con= 1.0/ ( 1.0 + m_mean*(8.0*z3-2.0*z3**2 )/zms**4 &
4997 + (1.0-m_mean)*( 20.0*z3-27.0*z3**2 +12.0*z3**3 -2.0*z3**4 )/(zms*(2.0-z3))**2 )
4999 fdsp = -2.0*pi*rho*i1*order1 - pi*rho*c1_con*m_mean*i2*order2
5001 END SUBROUTINE f_disp_pcsaft
5007 SUBROUTINE f_disp_cksaft ( fdsp )
5009 USE eos_variables, ONLY: nc, ncomp, pi, tau, t, rho, eta, x, z0t, mseg, vij, uij, parame, um
5014 REAL,
INTENT(IN OUT) :: fdsp
5016 INTEGER :: i, j, n, m
5017 REAL :: zmr, nmr, m_mean
5020 m_mean = z0t / ( pi / 6.0 )
5024 vij(i,j)=(0.5*((parame(i,2)*(1.0-0.12 *exp(-3.0*parame(i,3)/t))**3 )**(1.0/3.0) &
5025 +(parame(j,2)*(1.0-0.12 *exp(-3.0*parame(j,3)/t))**3 )**(1.0/3.0)))**3
5032 zmr = zmr + x(i)*x(j)*mseg(i)*mseg(j)*vij(i,j)*uij(i,j)
5033 nmr = nmr + x(i)*x(j)*mseg(i)*mseg(j)*vij(i,j)
5040 fdsp = fdsp + dnm(n,m) * (um/t)**n *(eta/tau)**m
5043 fdsp = m_mean * fdsp
5046 END SUBROUTINE f_disp_cksaft
5052 SUBROUTINE f_association ( eps_kij, k_kij, fhb )
5054 USE eos_variables, ONLY: nc, nsite, ncomp, t, z0t, z1t, z2t, z3t, rho, eta, x, &
5055 parame, sig_ij, dij_ab, gij, nhb_typ, mx, nhb_no
5059 REAL,
INTENT(IN) :: eps_kij, k_kij
5060 REAL,
INTENT(IN OUT) :: fhb
5063 INTEGER :: i, j, k, l, no, ass_cnt, max_eval
5064 REAL,
DIMENSION(nc,nc) :: kap_hb
5065 REAL,
DIMENSION(nc,nc,nsite,nsite) :: eps_hb
5066 REAL,
DIMENSION(nc,nsite,nc,nsite) :: delta
5067 REAL,
DIMENSION(nc,nsite) :: mx_itr
5068 REAL :: err_sum, sum0, amix, tol, ass_s1, ass_s2
5069 REAL :: z0, z1, z2, z3, zms
5074 IF (nint(parame(i,12)) /= 0) assoc = .true.
5087 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
5093 IF ( nint(parame(i,12)) /= 0 )
THEN 5094 nhb_typ(i) = nint( parame(i,12) )
5095 kap_hb(i,i) = parame(i,13)
5099 eps_hb(i,i,k,l) = parame(i,(14+no))
5104 nhb_no(i,k) = parame(i,(14+no))
5112 eps_hb(i,i,k,l) = 0.0
5120 IF ( i /= j .AND. (nhb_typ(i) /= 0.AND.nhb_typ(j) /= 0) )
THEN 5123 kap_hb(i,j) = (kap_hb(i,i)*kap_hb(j,j))**0.5 &
5124 *((parame(i,2)*parame(j,2))**3 )**0.5 &
5125 / (0.5*(parame(i,2)+parame(j,2)))**3
5126 kap_hb(i,j)= kap_hb(i,j)*(1.0-k_kij)
5129 IF ( k /= l .AND. nhb_typ(i) >= 2 .AND. nhb_typ(j) >= 2 )
THEN 5130 eps_hb(i,j,k,l) = (eps_hb(i,i,k,l)+eps_hb(j,j,l,k))/2.0
5132 eps_hb(i,j,k,l) = eps_hb(i,j,k,l)*(1.0-eps_kij)
5133 ELSE IF ( nhb_typ(i) == 1 .AND. l > k )
THEN 5134 eps_hb(i,j,k,l) = (eps_hb(i,i,k,k)+eps_hb(j,j,l,k))/2.0
5135 eps_hb(j,i,l,k) = (eps_hb(i,i,k,k)+eps_hb(j,j,l,k))/2.0
5136 eps_hb(i,j,k,l) = eps_hb(i,j,k,l)*(1.0-eps_kij)
5137 eps_hb(j,i,l,k) = eps_hb(j,i,l,k)*(1.0-eps_kij)
5147 IF ( parame(i,10) /= 0) kap_hb(i,i)=0.0
5149 IF ( parame(i,10) /= 0 .AND. parame(j,10) /= 0 ) kap_hb(i,j) = 0.0
5163 DO k = 1, nhb_typ(i)
5165 DO l = 1, nhb_typ(j)
5166 delta(i,k,j,l)=gij(i,j)*kap_hb(i,j)*(exp(eps_hb(i,j,k,l)/t)-1.0) *sig_ij(i,j)**3
5172 IF ( mx(i,k) == 0.0 ) mx(i,k) = 1.0
5179 IF (eta < 0.2) tol = 1.e-11
5180 IF (eta < 0.01) tol = 1.e-14
5187 ass_cnt = ass_cnt + 1
5190 DO k = 1, nhb_typ(i)
5193 DO l = 1, nhb_typ(j)
5194 sum0 = sum0 + x(j)* mx(j,l)*nhb_no(j,l) *delta(i,k,j,l)
5197 mx_itr(i,k) = 1.0 / (1.0 + sum0*rho)
5203 DO k = 1, nhb_typ(i)
5204 err_sum = err_sum + abs( mx_itr(i,k) - mx(i,k) )
5205 mx(i,k) = mx_itr(i,k) * amix + mx(i,k) * (1.0 - amix)
5206 IF ( mx(i,k) <= 0.0 ) mx(i,k)=1.e-50
5207 IF ( mx(i,k) > 1.0 ) mx(i,k)=1.0
5211 IF ( err_sum <= tol .OR. ass_cnt >= max_eval )
THEN 5212 IF ( ass_cnt >= max_eval )
WRITE(*,*)
'F_NUMERICAL: Max_eval violated = ',err_sum,tol
5221 DO k = 1, nhb_typ(i)
5222 ass_s1 = ass_s1 + nhb_no(i,k) * ( 1.0 - mx(i,k) )
5223 ass_s2 = ass_s2 + nhb_no(i,k) * log(mx(i,k))
5225 fhb = fhb + x(i) * ( ass_s2 + ass_s1/2.0 )
5230 END SUBROUTINE f_association
5237 SUBROUTINE f_ion_dipole_tbh ( fhend )
5239 USE eos_variables, ONLY: nc, pi, kbol, nav, ncomp, t, rho, eta, x, z0t, parame, uij, sig_ij
5243 REAL,
INTENT(IN OUT) :: fhend
5245 INTEGER :: i, dipole, ions
5247 REAL :: fh32, fh2, fh52, fh3
5248 REAL :: e_elem, eps_cc0, rho_sol, dielec
5249 REAL :: polabil, ydd, kappa, x_dipol, x_ions
5250 REAL,
DIMENSION(nc) :: my2dd, z_ii, e_cd, x_dd, x_ii
5251 REAL :: sig_c, sig_d, sig_cd, r_s
5252 REAL :: i0cc, i1cc, i2cc, icd, idd
5253 REAL :: iccc, iccd, icdd, iddd
5256 m_mean = z0t / ( pi / 6.0 )
5259 e_elem = 1.602189246e-19
5260 eps_cc0 = 8.854187818e-22
5265 rho_sol = rho * 18.015 * 1.e27/ nav
5266 rho_sol = rho_sol/1000.0
5267 dielec = 1.0+(7.62571/(t/293.15))*rho_sol +(2.44e2/(t/293.15)-1.41e2 &
5268 + 2.78e1*(t/293.15))*rho_sol**2 &
5269 + (-9.63e1/(t/293.15)+4.18e1*(t/293.15) &
5270 - 1.02e1*(t/293.15)**2 )*rho_sol**3 +(-4.52e1/(t/293.15)**2 &
5271 + 8.46e1/(t/293.15)-3.59e1)*rho_sol**4
5286 IF ( parame(i,6) /= 0.0 .AND. uij(i,i) /= 0.0 .AND. x(i) > 0.0 )
THEN 5287 my2dd(i) = (parame(i,6))**2 *1.e-49 / (uij(i,i)*kbol*sig_ij(i,i)**3 *1.e-30)
5293 z_ii(i) = parame(i,10)
5294 IF ( z_ii(i) /= 0.0 .AND. uij(i,i) /= 0.0 .AND. x(i) > 0.0 )
THEN 5295 e_cd(i) = ( parame(i,10)*e_elem* 1.e5 / sqrt(1.11265005) )**2 &
5296 / ( uij(i,i)*kbol*sig_ij(i,i)*1.e-10 )
5304 IF ( dipole == 1 .AND. ions == 1 )
THEN 5312 ydd = ydd + x(i)*(parame(i,6))**2 *1.e-49/ (kbol*t*1.e-30)
5313 kappa = kappa + x(i) &
5314 *(parame(i,10)*e_elem* 1.e5/sqrt(1.11265005))**2 /(kbol*t*1.e-10)
5315 IF (parame(i,10) /= 0.0)
THEN 5316 x_ions = x_ions + x(i)
5318 polabil = polabil + 4.0*pi*x(i)*rho*1.4573 *1.e-30 &
5319 / (sig_ij(3,3)**3 *1.e-30)
5321 IF (parame(i,6) /= 0.0) x_dipol= x_dipol+ x(i)
5323 ydd = ydd * 4.0/9.0 * pi * rho
5324 kappa = sqrt( 4.0 * pi * rho * kappa )
5332 IF(parame(i,10) /= 0.0 .AND. x_ions /= 0.0) x_ii(i) = x(i)/x_ions
5333 IF(parame(i,6) /= 0.0 .AND. x_dipol /= 0.0) x_dd(i) = x(i)/x_dipol
5334 sig_c = sig_c + x_ii(i)*parame(i,2)
5335 sig_d = sig_d + x_dd(i)*parame(i,2)
5337 sig_cd = 0.5 * (sig_c + sig_d)
5343 r_s = eta*6.0 / pi / m_mean
5345 i0cc = - (1.0 + 0.97743 * r_s + 0.05257*r_s*r_s) &
5346 /(1.0 + 1.43613 * r_s + 0.41580*r_s*r_s)
5348 i1cc = - (10.0 - 2.0*r_s*pi/6.0 + r_s*r_s*pi/6.0*pi/6.0) &
5349 /20.0/(1.0 + 2.0*r_s*pi/6.0)
5353 i2cc = -0.33331+0.7418*r_s - 1.2047*r_s*r_s &
5354 + 1.6139*r_s**3 - 1.5487*r_s**4 + 0.6626*r_s**5
5356 icd = (1.0 + 0.79576 *r_s + 0.104556 *r_s*r_s) &
5357 /(1.0 + 0.486704*r_s - 0.0222903*r_s*r_s)
5358 idd = (1.0 + 0.18158*r_s - 0.11467*r_s*r_s) &
5359 /3.0/(1.0 - 0.49303*r_s + 0.06293*r_s*r_s)
5360 iccc= 3.0*(1.0 - 1.05560*r_s + 0.26591*r_s*r_s) &
5361 /2.0/(1.0 + 0.53892*r_s - 0.94236*r_s*r_s)
5362 iccd= 11.0*(1.0 + 2.25642 *r_s + 0.05679 *r_s*r_s) &
5363 /6.0/(1.0 + 2.64178 *r_s + 0.79783 *r_s*r_s)
5364 icdd= 0.94685*(1.0 + 2.97323 *r_s + 3.11931 *r_s*r_s) &
5365 /(1.0 + 2.70186 *r_s + 1.22989 *r_s*r_s)
5366 iddd= 5.0*(1.0 + 1.12754*r_s + 0.56192*r_s*r_s) &
5367 /24.0/(1.0 - 0.05495*r_s + 0.13332*r_s*r_s)
5369 IF ( sig_c <= 0.0 )
WRITE (*,*)
'error in Henderson ion term' 5371 fh32= - kappa**3 /(12.0*pi*rho)
5372 fh2 = - 3.0* kappa**2 * ydd*icd /(8.0*pi*rho) / sig_cd &
5373 - kappa**4 *sig_c/(16.0*pi*rho)*i0cc
5374 IF (sig_d /= 0.0) fh2 = fh2 - 27.0* ydd * ydd*idd &
5375 /(8.0*pi*rho) / sig_d**3
5376 fh52= (3.0*kappa**3 * ydd + kappa**5 *sig_c**2 *i1cc) &
5378 fh3 = - kappa**6 * sig_c**3 /(8.0*pi*rho) *(i2cc-iccc/6.0) &
5379 + kappa**4 * ydd *sig_c/(16.0*pi*rho) &
5380 *( (6.0+5.0/3.0*sig_d/sig_c)*i0cc + 3.0*sig_d/sig_c*iccd ) &
5381 + 3.0*kappa**2 * ydd*ydd /(8.0*pi*rho) / sig_cd &
5382 *( (2.0-3.21555*sig_d/sig_cd)*icd +3.0*sig_d/sig_cd*icdd )
5383 IF (sig_d /= 0.0) fh3 = fh3 + 27.0*ydd**3 &
5384 /(16.0*pi*rho)/sig_d**3 *iddd
5386 fhend = ( fh32 + (fh32*fh32*fh3-2.0*fh32*fh2*fh52+fh2**3 ) &
5387 /(fh2*fh2-fh32*fh52) ) &
5388 / ( 1.0 + (fh32*fh3-fh2*fh52) /(fh2*fh2-fh32*fh52) &
5389 - (fh2*fh3-fh52*fh52) /(fh2*fh2-fh32*fh52) )
5406 END SUBROUTINE f_ion_dipole_tbh
5413 SUBROUTINE f_ion_ion_primmsa ( fcc )
5415 USE eos_variables, ONLY: nc, pi, kbol, nav, ncomp, t, rho, x, parame, mx
5419 REAL,
INTENT(IN OUT) :: fcc
5421 INTEGER :: i, j, cc_it, ions
5422 REAL :: e_elem, eps_cc0, rho_sol, dielec
5424 REAL :: cc_sig1, cc_sig2, cc_sig3
5425 REAL,
DIMENSION(nc) :: z_ii, x_ii, sigm_i, my2dd
5426 REAL :: alpha_2, kappa, ii_par
5427 REAL :: cc_omeg, p_n, q2_i, cc_q2, cc_gam
5428 REAL :: cc_error(2), cc_delt
5429 REAL :: rhs,
lambda, lam_s
5433 e_elem = 1.602189246e-19
5434 eps_cc0 = 8.854187818e-22
5439 rho_sol = rho * 18.015 * 1.e27/ nav
5440 rho_sol = rho_sol/1000.0
5441 dielec = 1.0+(7.62571/(t/293.15))*rho_sol +(2.44e2/(t/293.15)-1.41e2 &
5442 +2.78e1*(t/293.15))*rho_sol**2 &
5443 +(-9.63e1/(t/293.15)+4.18e1*(t/293.15) &
5444 -1.02e1*(t/293.15)**2 )*rho_sol**3 +(-4.52e1/(t/293.15)**2 &
5445 +8.46e1/(t/293.15)-3.59e1)*rho_sol**4
5475 z_ii(i) = parame(i,10)
5476 IF (z_ii(i) /= 0.0)
THEN 5477 sigm_i(i) = parame(i,2)
5481 IF (z_ii(i) /= 0.0) ions = 1
5482 IF (z_ii(i) /= 0.0) x_ions = x_ions + x(i)
5485 IF (ions == 1 .AND. x_ions > 0.0)
THEN 5491 IF (z_ii(i) /= 0.0)
THEN 5492 x_ii(i) = x(i)/x_ions
5496 cc_sig1 = cc_sig1 +x_ii(i)*sigm_i(i)
5497 cc_sig2 = cc_sig2 +x_ii(i)*sigm_i(i)**2
5498 cc_sig3 = cc_sig3 +x_ii(i)*sigm_i(i)**3
5503 alpha_2 = e_elem**2 /eps_cc0 / dielec / kbol/t
5506 kappa = kappa + x(i)*z_ii(i)*z_ii(i)*mx(i,1)
5508 kappa = sqrt( rho * alpha_2 * kappa )
5509 ii_par= kappa * cc_sig1
5522 cc_delt = cc_delt + x(i)*mx(i,1)*rho*sigm_i(i)**3
5524 cc_delt= 1.0 - pi/6.0*cc_delt
5534 cc_omeg = cc_omeg +x(i)*mx(i,1)*sigm_i(i)**3 /(1.0+cc_gam*sigm_i(i))
5536 cc_omeg = 1.0 + pi/2.0 / cc_delt * rho * cc_omeg
5539 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))
5544 q2_i = q2_i + rho*x(i)*mx(i,1)*( (z_ii(i)-pi/2.0/cc_delt*sigm_i(i)**2 *p_n) &
5545 /(1.0+cc_gam*sigm_i(i)) )**2
5546 cc_q2 = cc_q2 + x(i)*mx(i,1)*rho*z_ii(i)**2 / (1.0+cc_gam*sigm_i(i))
5548 q2_i = q2_i*alpha_2 / 4.0
5550 cc_error(j) = cc_gam - sqrt(q2_i)
5551 IF (j == 1) cc_gam = cc_gam*1.000001
5552 IF (j == 2) cc_gam = cc_gam - cc_error(2)* (cc_gam-cc_gam/1.000001)/(cc_error(2)-cc_error(1))
5554 IF ( j == 1 .AND. abs(cc_error(1)) > 1.e-15 )
GO TO 131
5555 IF ( cc_it >= 10 )
THEN 5556 WRITE (*,*)
' cc error' 5559 IF ( j /= 1 )
GO TO 13
5561 fcc= - alpha_2 / pi/4.0 /rho* (cc_gam*cc_q2 &
5562 + pi/2.0/cc_delt *cc_omeg*p_n**2 ) + cc_gam**3 /pi/3.0/rho
5570 my2dd(3) = (parame(3,6))**2 *1.e-19 /(kbol*t)
5571 my2dd(3) = (1.84)**2 *1.e-19 /(kbol*t)
5573 rhs = 12.0 * pi * rho * x(3) * my2dd(3)
5576 lambda = (rhs/((lam_s+2.0)**2 ) + 16.0/((1.0+lam_s)**4 ) )**0.5
5577 IF ( abs(lam_s-
lambda) > 1.e-10 )
THEN 5578 lam_s = (
lambda + lam_s ) / 2.0
5591 END SUBROUTINE f_ion_ion_primmsa
5598 SUBROUTINE f_ion_ion_nonprimmsa ( fdd, fqq, fdq, fcc )
5604 REAL,
INTENT(IN OUT) :: fdd, fqq, fdq, fcc
5608 REAL,
DIMENSION(nc) :: x_export, msegm
5612 IF ( sum( parame(1:ncomp,6) ) > 1.e-5 ) dipole = 1
5614 IF ( dipole /= 0 )
THEN 5625 write (*,*)
'why are individual contrib. A_CC,A_CD,A_DD not used' 5629 END SUBROUTINE f_ion_ion_nonprimmsa
5636 SUBROUTINE f_lc_mayersaupe ( flc )
5638 USE eos_variables, ONLY: nc, pi, kbol, nav, ncomp, phas, t, rho, eta, &
5639 x, mseg, parame, e_lc, s_lc, dhs
5643 REAL,
INTENT(IN OUT) :: flc
5646 INTEGER :: liq_crystal, count_lc, steps_lc
5647 REAL :: alpha_lc, tolerance, deltay
5648 REAL :: integrand1, integrand2, accel_lc
5649 REAL :: error_lc, u_term, sphase
5650 REAL,
DIMENSION(nc) :: z_lc, s_lc1, s_lc2, sumu
5651 REAL,
DIMENSION(nc,nc) :: u_lc, klc
5654 COMMON /stabil / stabil
5663 IF ( eta < 0.35 ) accel_lc = 1.3
5664 IF ( eta < 0.15 ) accel_lc = 1.0
5669 e_lc(i,j) = (e_lc(i,i)*e_lc(j,j))**0.5 *(1.0-klc(i,j))
5672 IF (e_lc(i,j) /= 0.0) liq_crystal = 1
5680 IF ( liq_crystal == 1 .AND. phas == 1 .AND. stabil == 0 )
THEN 5686 deltay = 1.0 /
REAL(steps_lc)
5692 u_lc(i,j) = 2.0/3.0*pi*mseg(i)*mseg(j) *(0.5*(dhs(i)+dhs(j)))**3 &
5693 *(e_lc(i,j)/t+s_lc(i,j))*rho
5707 IF (s_lc2(i) <= 0.3) s_lc1(i) = s_lc2(i)
5708 IF (s_lc2(i) > 0.3) s_lc1(i) = s_lc1(i) + (s_lc2(i)-s_lc1(i))*accel_lc
5711 count_lc = count_lc + 1
5718 sumu(i) = sumu(i) + x(j)*u_lc(i,j)*s_lc1(j)
5724 integrand1 = exp(-0.5*sumu(i))
5726 integrand2 = exp(0.5*sumu(i)*(3.0*(deltay*
REAL(k)) **2 -1.0))
5727 z_lc(i) = z_lc(i) + (integrand1 + integrand2)/2.0*deltay
5728 integrand1 = integrand2
5737 integrand1 = -1.0/z_lc(i)*0.5*exp(-0.5*sumu(i))
5739 integrand2 = 1.0/z_lc(i)*0.5*(3.0*(deltay*
REAL(k)) &
5740 **2 -1.0)*exp(0.5*sumu(i)*(3.0 *(deltay*
REAL(k))**2 -1.0))
5741 s_lc2(i) = s_lc2(i) + (integrand1 + integrand2)/2.0*deltay
5742 integrand1 = integrand2
5744 error_lc = error_lc + abs(s_lc2(i)-s_lc1(i))
5749 sphase = sphase + s_lc2(i)
5751 IF (sphase < 1.e-4)
THEN 5761 IF (error_lc > tolerance .AND. count_lc < 400)
GO TO 1
5764 IF (count_lc == 400)
WRITE (*,*)
'LC iteration not converg.' 5771 u_term = u_term + 0.5*x(i)*x(j)*s_lc2(i) *s_lc2(j)*u_lc(i,j)
5777 IF (z_lc(i) /= 0.0) flc = flc - x(i) * log(z_lc(i))
5786 END SUBROUTINE f_lc_mayersaupe
5794 SUBROUTINE p_polar ( zdd, zddz, zddz2, zddz3, zqq, zqqz, zqqz2, zqqz3, zdq, zdqz, zdqz2, zdqz3 )
5796 USE eos_variables, ONLY: ncomp, parame, dd_term, qq_term, dq_term
5800 REAL,
INTENT(OUT) :: zdd, zddz, zddz2, zddz3
5801 REAL,
INTENT(OUT) :: zqq, zqqz, zqqz2, zqqz3
5802 REAL,
INTENT(OUT) :: zdq, zdqz, zdqz2, zdqz3
5806 INTEGER :: quadrupole
5807 INTEGER :: dipole_quad
5826 IF ( sum( parame(1:ncomp,6) ) /= 0.0 ) dipole = 1
5827 IF ( sum( parame(1:ncomp,7) ) /= 0.0 ) quadrupole = 1
5828 IF ( dipole == 1 .AND. quadrupole == 1 ) dipole_quad = 1
5833 IF (dipole == 1)
THEN 5835 IF (dd_term ==
'GV')
CALL p_dd_gross_vrabec( zdd, zddz, zddz2, zddz3 )
5837 IF (dd_term /=
'GV' .AND. dd_term /=
'SF')
write (*,*)
'specify dipole term !' 5844 IF (quadrupole == 1)
THEN 5847 IF (qq_term ==
'JG')
CALL p_qq_gross( zqq, zqqz, zqqz2, zqqz3 )
5848 IF (qq_term /=
'JG' .AND. qq_term /=
'SF')
write (*,*)
'specify quadrupole term !' 5855 IF (dipole_quad == 1)
THEN 5857 IF (dq_term ==
'VG')
CALL p_dq_vrabec_gross( zdq, zdqz, zdqz2, zdqz3 )
5858 IF (dq_term /=
'VG' )
write (*,*)
'specify DQ-cross term !' 5862 END SUBROUTINE p_polar
5869 SUBROUTINE p_dd_gross_vrabec( zdd, zddz, zddz2, zddz3 )
5872 USE eos_variables, ONLY: nc, ncomp, uij, parame, mseg, sig_ij, rho, eta, x, t
5877 REAL,
INTENT(IN OUT) :: zdd, zddz, zddz2, zddz3
5879 INTEGER :: i, j, k, m
5880 REAL :: factor2, factor3, z3
5881 REAL :: xijfa, xijkfa, eij
5882 REAL :: fdddr, fddd2, fddd3, fddd4
5883 REAL :: fdd2, fdd2z, fdd2z2, fdd2z3, fdd2z4
5884 REAL :: fdd3, fdd3z, fdd3z2, fdd3z3, fdd3z4
5885 REAL,
DIMENSION(nc) :: my2dd
5886 REAL,
DIMENSION(nc,nc) :: idd2, idd2z, idd2z2, idd2z3, idd2z4
5887 REAL,
DIMENSION(nc,nc) :: idd4, idd4z, idd4z2, idd4z3, idd4z4
5888 REAL,
DIMENSION(nc,nc,nc) :: idd3, idd3z, idd3z2, idd3z3, idd3z4
5898 my2dd(i) = (parame(i,6))**2 *1.e-49 / (uij(i,i)*kbol*mseg(i)*sig_ij(i,i)**3 *1.e-30)
5915 idd2(i,j) =idd2(i,j) + ddp2(i,j,m) *z3**(m+1)
5916 idd4(i,j) =idd4(i,j) + ddp4(i,j,m) *z3**(m+1)
5917 idd2z(i,j) =idd2z(i,j) +ddp2(i,j,m)*
REAL(m+1) *z3**m
5918 idd4z(i,j) =idd4z(i,j) +ddp4(i,j,m)*
REAL(m+1) *z3**m
5919 idd2z2(i,j)=idd2z2(i,j)+ddp2(i,j,m)*
REAL((m+1)*m) *z3**(m-1)
5920 idd4z2(i,j)=idd4z2(i,j)+ddp4(i,j,m)*
REAL((m+1)*m) *z3**(m-1)
5921 idd2z3(i,j)=idd2z3(i,j)+ddp2(i,j,m)*
REAL((m+1)*m*(m-1)) *z3**(m-2)
5922 idd4z3(i,j)=idd4z3(i,j)+ddp4(i,j,m)*
REAL((m+1)*m*(m-1)) *z3**(m-2)
5923 idd2z4(i,j)=idd2z4(i,j)+ddp2(i,j,m)*
REAL((m+1)*m*(m-1)*(m-2)) *z3**(m-3)
5924 idd4z4(i,j)=idd4z4(i,j)+ddp4(i,j,m)*
REAL((m+1)*m*(m-1)*(m-2)) *z3**(m-3)
5934 idd3(i,j,k) =idd3(i,j,k) +ddp3(i,j,k,m)*z3**(m+2)
5935 idd3z(i,j,k) =idd3z(i,j,k) +ddp3(i,j,k,m)*
REAL(m+2)*z3**(m+1)
5936 idd3z2(i,j,k)=idd3z2(i,j,k)+ddp3(i,j,k,m)*
REAL((m+2)*(m+1))*z3**m
5937 idd3z3(i,j,k)=idd3z3(i,j,k)+ddp3(i,j,k,m)*
REAL((m+2)*(m+1)*m)*z3**(m-1)
5938 idd3z4(i,j,k)=idd3z4(i,j,k)+ddp3(i,j,k,m)*
REAL((m+2)*(m+1)*m*(m-1)) *z3**(m-2)
5946 factor2= -pi *rho/z3
5947 factor3= -4.0/3.0*pi**2 * (rho/z3)**2
5962 xijfa = x(i)*parame(i,3)/t*parame(i,2)**3 *x(j)*parame(j,3)/t*parame(j,2)**3 &
5963 / ((parame(i,2)+parame(j,2))/2.0)**3 *my2dd(i)*my2dd(j)
5964 eij = (parame(i,3)*parame(j,3))**0.5
5965 fdd2 = fdd2 +factor2*xijfa*(idd2(i,j) +eij/t*idd4(i,j))
5966 fdd2z = fdd2z +factor2*xijfa*(idd2z(i,j) +eij/t*idd4z(i,j))
5967 fdd2z2 = fdd2z2+factor2*xijfa*(idd2z2(i,j)+eij/t*idd4z2(i,j))
5968 fdd2z3 = fdd2z3+factor2*xijfa*(idd2z3(i,j)+eij/t*idd4z3(i,j))
5969 fdd2z4 = fdd2z4+factor2*xijfa*(idd2z4(i,j)+eij/t*idd4z4(i,j))
5972 xijkfa= x(i)*parame(i,3)/t*parame(i,2)**3 *x(j)*parame(j,3)/t*parame(j,2)**3 &
5973 *x(k)*parame(k,3)/t*parame(k,2)**3 /((parame(i,2)+parame(j,2))/2.0) &
5974 /((parame(i,2)+parame(k,2))/2.0) /((parame(j,2)+parame(k,2))/2.0) &
5975 *my2dd(i)*my2dd(j)*my2dd(k)
5976 fdd3 = fdd3 + factor3 * xijkfa*idd3(i,j,k)
5977 fdd3z = fdd3z + factor3 * xijkfa*idd3z(i,j,k)
5978 fdd3z2 = fdd3z2 + factor3 * xijkfa*idd3z2(i,j,k)
5979 fdd3z3 = fdd3z3 + factor3 * xijkfa*idd3z3(i,j,k)
5980 fdd3z4 = fdd3z4 + factor3 * xijkfa*idd3z4(i,j,k)
5986 IF (fdd2 < -1.e-50 .AND. fdd3 /= 0.0 .AND. fdd2z /= 0.0 .AND. fdd3z /= 0.0)
THEN 5988 fdddr= fdd2* (fdd2*fdd2z - 2.0*fdd3*fdd2z+fdd2*fdd3z) / (fdd2-fdd3)**2
5989 fddd2=(2.0*fdd2*fdd2z*fdd2z +fdd2*fdd2*fdd2z2 &
5990 -2.0*fdd2z**2 *fdd3-2.0*fdd2*fdd2z2*fdd3+fdd2*fdd2*fdd3z2) &
5991 /(fdd2-fdd3)**2 + fdddr * 2.0*(fdd3z-fdd2z)/(fdd2-fdd3)
5992 fddd3=(2.0*fdd2z**3 +6.0*fdd2*fdd2z*fdd2z2+fdd2*fdd2*fdd2z3 &
5993 -6.0*fdd2z*fdd2z2*fdd3-2.0*fdd2z**2 *fdd3z &
5994 -2.0*fdd2*fdd2z3*fdd3 -2.0*fdd2*fdd2z2*fdd3z &
5995 +2.0*fdd2*fdd2z*fdd3z2+fdd2*fdd2*fdd3z3) /(fdd2-fdd3)**2 &
5996 + 2.0/(fdd2-fdd3)* ( 2.0*fddd2*(fdd3z-fdd2z) &
5997 + fdddr*(fdd3z2-fdd2z2) &
5998 - fdddr/(fdd2-fdd3)*(fdd3z-fdd2z)**2 )
5999 fddd4=( 12.0*fdd2z**2 *fdd2z2+6.0*fdd2*fdd2z2**2 &
6000 +8.0*fdd2*fdd2z*fdd2z3+fdd2*fdd2*fdd2z4-6.0*fdd2z2**2 *fdd3 &
6001 -12.0*fdd2z*fdd2z2*fdd3z -8.0*fdd2z*fdd2z3*fdd3 &
6002 -2.0*fdd2*fdd2z4*fdd3-4.0*fdd2*fdd2z3*fdd3z &
6003 +4.0*fdd2*fdd2z*fdd3z3+fdd2**2 *fdd3z4 ) /(fdd2-fdd3)**2 &
6004 + 6.0/(fdd2-fdd3)* ( fddd3*(fdd3z-fdd2z) &
6005 -fddd2/(fdd2-fdd3)*(fdd3z-fdd2z)**2 &
6006 - fdddr/(fdd2-fdd3)*(fdd3z-fdd2z)*(fdd3z2-fdd2z2) &
6007 + fddd2*(fdd3z2-fdd2z2) +1.0/3.0*fdddr*(fdd3z3-fdd2z3) )
6009 zddz = fddd2*eta + fdddr
6010 zddz2 = fddd3*eta + 2.0* fddd2
6011 zddz3 = fddd4*eta + 3.0* fddd3
6016 END SUBROUTINE p_dd_gross_vrabec
6024 SUBROUTINE p_qq_gross( zqq, zqqz, zqqz2, zqqz3 )
6027 USE eos_variables, ONLY: nc, ncomp, uij, parame, mseg, sig_ij, rho, eta, x, t
6032 REAL,
INTENT(IN OUT) :: zqq, zqqz, zqqz2, zqqz3
6034 INTEGER :: i, j, k, m
6035 REAL :: factor2, factor3, z3
6036 REAL :: xijfa, xijkfa, eij
6037 REAL :: fqqdr, fqqd2, fqqd3, fqqd4
6038 REAL :: fqq2, fqq2z, fqq2z2, fqq2z3, fqq2z4
6039 REAL :: fqq3, fqq3z, fqq3z2, fqq3z3, fqq3z4
6040 REAL,
DIMENSION(nc) :: qq2
6041 REAL,
DIMENSION(nc,nc) :: iqq2, iqq2z, iqq2z2, iqq2z3, iqq2z4
6042 REAL,
DIMENSION(nc,nc) :: iqq4, iqq4z, iqq4z2, iqq4z3, iqq4z4
6043 REAL,
DIMENSION(nc,nc,nc) :: iqq3, iqq3z, iqq3z2, iqq3z3, iqq3z4
6052 qq2(i) = (parame(i,7))**2 *1.e-69 / (uij(i,i)*kbol*mseg(i)*sig_ij(i,i)**5 *1.e-50)
6067 IF (parame(i,7) /= 0.0 .AND. parame(j,7) /= 0.0)
THEN 6069 iqq2(i,j) =iqq2(i,j) + qqp2(i,j,m)*z3**(m+1)
6070 iqq4(i,j) =iqq4(i,j) + qqp4(i,j,m)*z3**(m+1)
6071 iqq2z(i,j) =iqq2z(i,j) +qqp2(i,j,m)*
REAL(m+1)*z3**m
6072 iqq4z(i,j) =iqq4z(i,j) +qqp4(i,j,m)*
REAL(m+1)*z3**m
6073 iqq2z2(i,j)=iqq2z2(i,j)+qqp2(i,j,m)*
REAL((m+1)*m)*z3**(m-1)
6074 iqq4z2(i,j)=iqq4z2(i,j)+qqp4(i,j,m)*
REAL((m+1)*m)*z3**(m-1)
6075 iqq2z3(i,j)=iqq2z3(i,j)+qqp2(i,j,m)*
REAL((m+1)*m*(m-1)) *z3**(m-2)
6076 iqq4z3(i,j)=iqq4z3(i,j)+qqp4(i,j,m)*
REAL((m+1)*m*(m-1)) *z3**(m-2)
6077 iqq2z4(i,j)=iqq2z4(i,j)+qqp2(i,j,m)*
REAL((m+1)*m*(m-1)*(m-2)) *z3**(m-3)
6078 iqq4z4(i,j)=iqq4z4(i,j)+qqp4(i,j,m)*
REAL((m+1)*m*(m-1)*(m-2)) *z3**(m-3)
6086 IF (parame(k,7) /= 0.0)
THEN 6088 iqq3(i,j,k) =iqq3(i,j,k) + qqp3(i,j,k,m)*z3**(m+2)
6089 iqq3z(i,j,k)=iqq3z(i,j,k)+qqp3(i,j,k,m)*
REAL(m+2)*z3**(m+1)
6090 iqq3z2(i,j,k)=iqq3z2(i,j,k)+qqp3(i,j,k,m)*
REAL((m+2)*(m+1)) *z3**m
6091 iqq3z3(i,j,k)=iqq3z3(i,j,k)+qqp3(i,j,k,m)*
REAL((m+2)*(m+1)*m) *z3**(m-1)
6092 iqq3z4(i,j,k)=iqq3z4(i,j,k)+qqp3(i,j,k,m) *
REAL((m+2)*(m+1)*m*(m-1)) *z3**(m-2)
6101 factor2= -9.0/16.0*pi *rho/z3
6102 factor3= 9.0/16.0*pi**2 * (rho/z3)**2
6116 IF (parame(i,7) /= 0.0 .AND. parame(j,7) /= 0.0)
THEN 6117 xijfa =x(i)*uij(i,i)*qq2(i)*sig_ij(i,i)**5 /t &
6118 *x(j)*uij(j,j)*qq2(j)*sig_ij(j,j)**5 /t/sig_ij(i,j)**7.0
6119 eij = (parame(i,3)*parame(j,3))**0.5
6120 fqq2= fqq2 +factor2*xijfa*(iqq2(i,j) +eij/t*iqq4(i,j) )
6121 fqq2z =fqq2z +factor2*xijfa*(iqq2z(i,j) +eij/t*iqq4z(i,j) )
6122 fqq2z2=fqq2z2+factor2*xijfa*(iqq2z2(i,j)+eij/t*iqq4z2(i,j))
6123 fqq2z3=fqq2z3+factor2*xijfa*(iqq2z3(i,j)+eij/t*iqq4z3(i,j))
6124 fqq2z4=fqq2z4+factor2*xijfa*(iqq2z4(i,j)+eij/t*iqq4z4(i,j))
6126 IF (parame(k,7) /= 0.0)
THEN 6127 xijkfa=x(i)*uij(i,i)*qq2(i)*sig_ij(i,i)**5 /t/sig_ij(i,j)**3 &
6128 *x(j)*uij(j,j)*qq2(j)*sig_ij(j,j)**5 /t/sig_ij(i,k)**3 &
6129 *x(k)*uij(k,k)*qq2(k)*sig_ij(k,k)**5 /t/sig_ij(j,k)**3
6130 fqq3 = fqq3 + factor3 * xijkfa*iqq3(i,j,k)
6131 fqq3z = fqq3z + factor3 * xijkfa*iqq3z(i,j,k)
6132 fqq3z2 = fqq3z2 + factor3 * xijkfa*iqq3z2(i,j,k)
6133 fqq3z3 = fqq3z3 + factor3 * xijkfa*iqq3z3(i,j,k)
6134 fqq3z4 = fqq3z4 + factor3 * xijkfa*iqq3z4(i,j,k)
6140 IF (fqq2 < -1.e-50 .AND. fqq3 /= 0.0 .AND. fqq2z /= 0.0 .AND. fqq3z /= 0.0)
THEN 6141 fqqdr = fqq2* (fqq2*fqq2z - 2.0*fqq3*fqq2z+fqq2*fqq3z) /(fqq2-fqq3)**2
6142 fqqd2= (2.0*fqq2*fqq2z*fqq2z +fqq2*fqq2*fqq2z2 &
6143 -2.0*fqq2z**2 *fqq3-2.0*fqq2*fqq2z2*fqq3+fqq2*fqq2*fqq3z2) &
6144 /(fqq2-fqq3)**2 + fqqdr * 2.0*(fqq3z-fqq2z)/(fqq2-fqq3)
6145 fqqd3=(2.0*fqq2z**3 +6.0*fqq2*fqq2z*fqq2z2+fqq2*fqq2*fqq2z3 &
6146 -6.0*fqq2z*fqq2z2*fqq3-2.0*fqq2z**2 *fqq3z &
6147 -2.0*fqq2*fqq2z3*fqq3 -2.0*fqq2*fqq2z2*fqq3z &
6148 +2.0*fqq2*fqq2z*fqq3z2+fqq2*fqq2*fqq3z3) /(fqq2-fqq3)**2 &
6149 + 2.0/(fqq2-fqq3)* ( 2.0*fqqd2*(fqq3z-fqq2z) &
6150 + fqqdr*(fqq3z2-fqq2z2) - fqqdr/(fqq2-fqq3)*(fqq3z-fqq2z)**2 )
6151 fqqd4=( 12.0*fqq2z**2 *fqq2z2+6.0*fqq2*fqq2z2**2 &
6152 +8.0*fqq2*fqq2z*fqq2z3+fqq2*fqq2*fqq2z4-6.0*fqq2z2**2 *fqq3 &
6153 -12.0*fqq2z*fqq2z2*fqq3z -8.0*fqq2z*fqq2z3*fqq3 &
6154 -2.0*fqq2*fqq2z4*fqq3-4.0*fqq2*fqq2z3*fqq3z &
6155 +4.0*fqq2*fqq2z*fqq3z3+fqq2**2 *fqq3z4 ) /(fqq2-fqq3)**2 &
6156 + 6.0/(fqq2-fqq3)* ( fqqd3*(fqq3z-fqq2z) &
6157 -fqqd2/(fqq2-fqq3)*(fqq3z-fqq2z)**2 &
6158 - fqqdr/(fqq2-fqq3)*(fqq3z-fqq2z)*(fqq3z2-fqq2z2) &
6159 + fqqd2*(fqq3z2-fqq2z2) +1.0/3.0*fqqdr*(fqq3z3-fqq2z3) )
6161 zqqz = fqqd2*eta + fqqdr
6162 zqqz2 = fqqd3*eta + 2.0* fqqd2
6163 zqqz3 = fqqd4*eta + 3.0* fqqd3
6167 END SUBROUTINE p_qq_gross
6173 SUBROUTINE p_dq_vrabec_gross( zdq, zdqz, zdqz2, zdqz3 )
6176 USE eos_variables, ONLY: nc, ncomp, uij, parame, mseg, sig_ij, rho, eta, x, t
6181 REAL,
INTENT(IN OUT) :: zdq, zdqz, zdqz2, zdqz3
6183 INTEGER :: i, j, k, m
6184 REAL :: factor2, factor3, z3
6185 REAL :: xijfa, xijkfa, eij
6186 REAL :: fdqdr, fdqd2, fdqd3, fdqd4
6187 REAL :: fdq2, fdq2z, fdq2z2, fdq2z3, fdq2z4
6188 REAL :: fdq3, fdq3z, fdq3z2, fdq3z3, fdq3z4
6189 REAL,
DIMENSION(nc) :: my2dd, myfac, qq2, q_fac
6190 REAL,
DIMENSION(nc,nc) :: idq2, idq2z, idq2z2, idq2z3, idq2z4
6191 REAL,
DIMENSION(nc,nc) :: idq4, idq4z, idq4z2, idq4z3, idq4z4
6192 REAL,
DIMENSION(nc,nc,nc) :: idq3, idq3z, idq3z2, idq3z3, idq3z4
6201 my2dd(i) = (parame(i,6))**2 *1.e-49 / (uij(i,i)*kbol*mseg(i)*sig_ij(i,i)**3 *1.e-30)
6202 myfac(i) = parame(i,3)/t*parame(i,2)**4 *my2dd(i)
6203 qq2(i) = (parame(i,7))**2 *1.e-69 / (uij(i,i)*kbol*mseg(i)*sig_ij(i,i)**5 *1.e-50)
6204 q_fac(i) = parame(i,3)/t*parame(i,2)**4 *qq2(i)
6220 IF (myfac(i) /= 0.0 .AND. q_fac(j) /= 0.0)
THEN 6222 idq2(i,j) =idq2(i,j) + dqp2(i,j,m)*z3**(m+1)
6223 idq4(i,j) =idq4(i,j) + dqp4(i,j,m)*z3**(m+1)
6224 idq2z(i,j) =idq2z(i,j) +dqp2(i,j,m)*
REAL(m+1)*z3**m
6225 idq4z(i,j) =idq4z(i,j) +dqp4(i,j,m)*
REAL(m+1)*z3**m
6226 idq2z2(i,j)=idq2z2(i,j)+dqp2(i,j,m)*
REAL((m+1)*m)*z3**(m-1)
6227 idq4z2(i,j)=idq4z2(i,j)+dqp4(i,j,m)*
REAL((m+1)*m)*z3**(m-1)
6228 idq2z3(i,j)=idq2z3(i,j)+dqp2(i,j,m)*
REAL((m+1)*m*(m-1)) *z3**(m-2)
6229 idq4z3(i,j)=idq4z3(i,j)+dqp4(i,j,m)*
REAL((m+1)*m*(m-1)) *z3**(m-2)
6230 idq2z4(i,j)=idq2z4(i,j)+dqp2(i,j,m)*
REAL((m+1)*m*(m-1)*(m-2)) *z3**(m-3)
6231 idq4z4(i,j)=idq4z4(i,j)+dqp4(i,j,m)*
REAL((m+1)*m*(m-1)*(m-2)) *z3**(m-3)
6239 IF (myfac(k) /= 0.0.OR.q_fac(k) /= 0.0)
THEN 6241 idq3(i,j,k) =idq3(i,j,k) + dqp3(i,j,k,m)*z3**(m+2)
6242 idq3z(i,j,k)=idq3z(i,j,k)+dqp3(i,j,k,m)*
REAL(m+2)*z3**(m+1)
6243 idq3z2(i,j,k)=idq3z2(i,j,k)+dqp3(i,j,k,m)*
REAL((m+2)*(m+1)) *z3**m
6244 idq3z3(i,j,k)=idq3z3(i,j,k)+dqp3(i,j,k,m)*
REAL((m+2)*(m+1)*m) *z3**(m-1)
6245 idq3z4(i,j,k)=idq3z4(i,j,k)+dqp3(i,j,k,m) &
6246 *
REAL((m+2)*(m+1)*m*(m-1)) *z3**(m-2)
6255 factor2= -9.0/4.0*pi *rho/z3
6256 factor3= pi**2 * (rho/z3)**2
6270 IF (myfac(i) /= 0.0 .AND. q_fac(j) /= 0.0)
THEN 6271 xijfa =x(i)*myfac(i) * x(j)*q_fac(j) /sig_ij(i,j)**5
6272 eij = (parame(i,3)*parame(j,3))**0.5
6273 fdq2= fdq2 +factor2*xijfa*(idq2(i,j) +eij/t*idq4(i,j) )
6274 fdq2z =fdq2z +factor2*xijfa*(idq2z(i,j) +eij/t*idq4z(i,j) )
6275 fdq2z2=fdq2z2+factor2*xijfa*(idq2z2(i,j)+eij/t*idq4z2(i,j))
6276 fdq2z3=fdq2z3+factor2*xijfa*(idq2z3(i,j)+eij/t*idq4z3(i,j))
6277 fdq2z4=fdq2z4+factor2*xijfa*(idq2z4(i,j)+eij/t*idq4z4(i,j))
6279 IF (myfac(k) /= 0.0.OR.q_fac(k) /= 0.0)
THEN 6280 xijkfa=x(i)*x(j)*x(k)/(sig_ij(i,j)*sig_ij(i,k)*sig_ij(j,k))**2 &
6281 *( myfac(i)*q_fac(j)*myfac(k) + myfac(i)*q_fac(j)*q_fac(k)*1.193735 )
6282 fdq3 =fdq3 + factor3 * xijkfa*idq3(i,j,k)
6283 fdq3z =fdq3z + factor3 * xijkfa*idq3z(i,j,k)
6284 fdq3z2=fdq3z2 + factor3 * xijkfa*idq3z2(i,j,k)
6285 fdq3z3=fdq3z3 + factor3 * xijkfa*idq3z3(i,j,k)
6286 fdq3z4=fdq3z4 + factor3 * xijkfa*idq3z4(i,j,k)
6292 IF (fdq2 < -1.e-50 .AND. fdq3 /= 0.0 .AND. fdq2z /= 0.0 .AND. fdq3z /= 0.0)
THEN 6293 fdqdr = fdq2* (fdq2*fdq2z - 2.0*fdq3*fdq2z+fdq2*fdq3z) /(fdq2-fdq3)**2
6294 fdqd2= (2.0*fdq2*fdq2z*fdq2z +fdq2*fdq2*fdq2z2 &
6295 -2.0*fdq2z**2 *fdq3-2.0*fdq2*fdq2z2*fdq3+fdq2*fdq2*fdq3z2) &
6296 /(fdq2-fdq3)**2 + fdqdr * 2.0*(fdq3z-fdq2z)/(fdq2-fdq3)
6297 fdqd3=(2.0*fdq2z**3 +6.0*fdq2*fdq2z*fdq2z2+fdq2*fdq2*fdq2z3 &
6298 -6.0*fdq2z*fdq2z2*fdq3-2.0*fdq2z**2 *fdq3z &
6299 -2.0*fdq2*fdq2z3*fdq3 -2.0*fdq2*fdq2z2*fdq3z &
6300 +2.0*fdq2*fdq2z*fdq3z2+fdq2*fdq2*fdq3z3) /(fdq2-fdq3)**2 &
6301 + 2.0/(fdq2-fdq3)* ( 2.0*fdqd2*(fdq3z-fdq2z) &
6302 + fdqdr*(fdq3z2-fdq2z2) - fdqdr/(fdq2-fdq3)*(fdq3z-fdq2z)**2 )
6303 fdqd4=( 12.0*fdq2z**2 *fdq2z2+6.0*fdq2*fdq2z2**2 &
6304 +8.0*fdq2*fdq2z*fdq2z3+fdq2*fdq2*fdq2z4-6.0*fdq2z2**2 *fdq3 &
6305 -12.0*fdq2z*fdq2z2*fdq3z -8.0*fdq2z*fdq2z3*fdq3 &
6306 -2.0*fdq2*fdq2z4*fdq3-4.0*fdq2*fdq2z3*fdq3z &
6307 +4.0*fdq2*fdq2z*fdq3z3+fdq2**2 *fdq3z4 ) /(fdq2-fdq3)**2 &
6308 + 6.0/(fdq2-fdq3)* ( fdqd3*(fdq3z-fdq2z) &
6309 -fdqd2/(fdq2-fdq3)*(fdq3z-fdq2z)**2 &
6310 - fdqdr/(fdq2-fdq3)*(fdq3z-fdq2z)*(fdq3z2-fdq2z2) &
6311 + fdqd2*(fdq3z2-fdq2z2) +1.0/3.0*fdqdr*(fdq3z3-fdq2z3) )
6313 zdqz = fdqd2*eta + fdqdr
6314 zdqz2 = fdqd3*eta + 2.0* fdqd2
6315 zdqz3 = fdqd4*eta + 3.0* fdqd3
6319 END SUBROUTINE p_dq_vrabec_gross
6327 SUBROUTINE f_pert_theory ( fdsp )
6330 x, z0t, mseg, parame, order1, order2
6336 REAL,
INTENT(IN OUT) :: fdsp
6339 REAL :: z3, zms, c1_con, m_mean
6345 IF (disp_term ==
'PT1')
THEN 6347 CALL f_dft ( i1, i2)
6350 fdsp = + ( - 2.0*pi*rho*i1*order1 )
6352 ELSEIF (disp_term ==
'PT2')
THEN 6354 CALL f_dft ( i1, i2)
6357 m_mean = z0t / ( pi / 6.0 )
6358 c1_con = 1.0/ ( 1.0 + m_mean*(8.0*z3-2.0*z3**2 )/zms**4 &
6359 + (1.0 - m_mean)*( 20.0*z3 -27.0*z3**2 +12.0*z3**3 -2.0*z3**4 ) &
6360 /(zms*(2.0-z3))**2 )
6361 fdsp = + ( - 2.0*pi*rho*i1*order1 - pi*rho*c1_con*m_mean*i2*order2 )
6363 ELSEIF (disp_term ==
'PT_MIX')
THEN 6365 CALL f_pert_theory_mix ( fdsp )
6367 ELSEIF (disp_term ==
'PT_MF')
THEN 6370 i1 = - ( - 8.0/9.0 - 4.0/9.0*(rc**(-9) -3.0*rc**(-3) ) - tau_cut/3.0*(rc**3 -1.0) )
6371 fdsp = + ( - 2.0*pi*rho*i1*order1 )
6372 write (*,*)
'caution: not thoroughly checked and tested' 6375 write (*,*)
'define the type of perturbation theory' 6382 END SUBROUTINE f_pert_theory
6391 SUBROUTINE f_pert_theory_mix ( fdsp )
6393 USE eos_variables, ONLY: nc, pi, ncomp, t, rho, eta, x, parame, mseg, dhs, sig_ij, uij
6398 REAL,
INTENT(IN OUT) :: fdsp
6404 REAL :: ua, ua_c, rm
6405 REAL,
DIMENSION(nc,nc) :: i1
6406 REAL :: int10, int11
6407 REAL :: d_ij, dzr_local
6408 REAL :: rad, xg, rdf
6409 REAL :: dg_dz3, dg_dr
6414 ua_c = 4.0 * ( rc**(-12) - rc**(-6) )
6424 int10 = rc * rc * ua_c
6430 DO WHILE ( rad /= 1.0 )
6433 IF ( rad - dzr_local <= 1.0 ) dzr_local = rad - 1.0
6435 rad = rad - dzr_local
6439 d_ij = 0.5*(dhs(l)+dhs(m)) / sig_ij(l,m)
6444 IF ( rad <= rg )
THEN 6445 IF ( l == 1 .AND. m == 1 )
CALL bi_cub_spline (z3,xg,ya_11,x1a_11,x2a_11,y1a_11,y2a_11,y12a_11, &
6446 c_bicub_11,rdf,dg_dz3,dg_dr,den_step,ih,k)
6447 IF ( l /= m )
CALL bi_cub_spline (z3,xg,ya_12,x1a_12,x2a_12,y1a_12,y2a_12,y12a_12, &
6448 c_bicub_12,rdf,dg_dz3,dg_dr,den_step,ih,k)
6449 IF ( l == 2 .AND. m == 2 )
CALL bi_cub_spline (z3,xg,ya_22,x1a_22,x2a_22,y1a_22,y2a_22,y12a_22, &
6450 c_bicub_22,rdf,dg_dz3,dg_dr,den_step,ih,k)
6453 ua = 4.0 * ( rad**(-12) - rad**(-6) )
6455 int11 = rdf * rad * rad * ua
6456 i1(l,m) = i1(l,m) + dzr_local * ( int11 + int10 ) / 2.0
6482 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)
6502 END SUBROUTINE f_pert_theory_mix
6509 SUBROUTINE mu_pert_theory_mix ( mu_dsp )
6511 USE eos_variables, ONLY: nc, pi, ncomp, t, rho, eta, x, parame, mseg, dhs, sig_ij, uij
6516 REAL,
INTENT(IN OUT) :: mu_dsp(nc)
6522 REAL :: ua, ua_c, rm
6523 REAL,
DIMENSION(nc,nc) :: i1, i2
6524 REAL :: int1_0, int1_1, int2_0, int2_1
6525 REAL :: d_ij, dzr_local
6526 REAL :: rad, xg, rdf
6527 REAL :: dg_dz3, dg_dr
6528 REAL :: term1(nc), term2
6533 ua_c = 4.0 * ( rc**(-12) - rc**(-6) )
6547 int1_0 = rc * rc * ua_c
6553 DO WHILE ( rad /= 1.0 )
6556 IF ( rad - dzr_local <= 1.0 ) dzr_local = rad - 1.0
6558 rad = rad - dzr_local
6561 d_ij = 0.5*(dhs(l)+dhs(m)) / sig_ij(l,m)
6566 IF ( rad <= rg )
THEN 6567 IF ( l == 1 .AND. m == 1 )
CALL bi_cub_spline (z3,xg,ya_11,x1a_11,x2a_11,y1a_11,y2a_11,y12a_11, &
6568 c_bicub_11,rdf,dg_dz3,dg_dr,den_step,ih,k)
6569 IF ( l /= m )
CALL bi_cub_spline (z3,xg,ya_12,x1a_12,x2a_12,y1a_12,y2a_12,y12a_12, &
6570 c_bicub_12,rdf,dg_dz3,dg_dr,den_step,ih,k)
6571 IF ( l == 2 .AND. m == 2 )
CALL bi_cub_spline (z3,xg,ya_22,x1a_22,x2a_22,y1a_22,y2a_22,y12a_22, &
6572 c_bicub_22,rdf,dg_dz3,dg_dr,den_step,ih,k)
6575 ua = 4.0 * ( rad**(-12) - rad**(-6) )
6577 int1_1 = rdf * rad * rad * ua
6578 int2_1 = dg_dz3 * rad * rad * ua
6579 i1(l,m) = i1(l,m) + dzr_local * ( int1_1 + int1_0 ) / 2.0
6580 i2(l,m) = i2(l,m) + dzr_local * ( int2_1 + int2_0 ) / 2.0
6585 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
6603 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)
6608 mu_dsp(l) = term1(l) + term2 * pi/ 6.0 * mseg(l)*dhs(l)**3
6611 END SUBROUTINE mu_pert_theory_mix
6618 SUBROUTINE f_dd_gross_vrabec( fdd )
6621 USE eos_variables, ONLY: nc, ncomp, uij, parame, mseg, sig_ij, rho, eta, x, t
6626 REAL,
INTENT(IN OUT) :: fdd
6628 INTEGER :: i, j, k, m
6629 INTEGER :: ddit, ddmax
6630 REAL :: factor2, factor3
6631 REAL :: xijfa, xijkfa, xijf_j, xijkf_j, eij
6633 REAL,
DIMENSION(nc) :: my2dd, my0, alph_tst, z1dd, z2dd, dderror
6634 REAL,
DIMENSION(nc) :: fdd2m, fdd3m, fdd2m2, fdd3m2, fddm, fddm2
6635 REAL,
DIMENSION(nc,nc) :: idd2, idd4
6636 REAL,
DIMENSION(nc,nc,nc) :: idd3
6644 IF ( uij(i,i) == 0.0 )
write (*,*)
'F_DD_GROSS_VRABEC: do not use dimensionless units' 6645 IF ( uij(i,i) == 0.0 ) stop
6646 my2dd(i) = (parame(i,6))**2 *1.e-49 /(uij(i,i)*kbol* mseg(i)*sig_ij(i,i)**3 *1.e-30)
6647 alph_tst(i) = parame(i,11) / (mseg(i)*sig_ij(i,i)**3 ) * t/parame(i,3)
6648 IF ( alph_tst(i) /= 0.0 ) ddmax = 25
6649 z1dd(i) = my2dd(i) + 3.0*alph_tst(i)
6650 z2dd(i) = 3.0*alph_tst(i)
6660 idd2(i,j) = idd2(i,j) + ddp2(i,j,m)*eta**m
6661 idd4(i,j) = idd4(i,j) + ddp4(i,j,m)*eta**m
6667 idd3(i,j,k) = idd3(i,j,k) + ddp3(i,j,k,m)*eta**m
6676 factor3 = -4.0/3.0*pi**2 * rho**2
6689 xijfa =x(i)*parame(i,3)/t*parame(i,2)**3 * x(j)*parame(j,3)/t*parame(j,2)**3 &
6690 /((parame(i,2)+parame(j,2))/2.0)**3 * (z1dd(i)*z1dd(j)-z2dd(i)*z2dd(j))
6691 eij = (parame(i,3)*parame(j,3))**0.5
6692 fdd2= fdd2 + factor2 * xijfa * ( idd2(i,j) + eij/t*idd4(i,j) )
6693 xijf_j = parame(i,3)/t*parame(i,2)**3 *x(j)*parame(j,3)/t*parame(j,2)**3 &
6694 /((parame(i,2)+parame(j,2))/2.0)**3
6695 fdd2m(i)=fdd2m(i)+4.0*sqrt(my2dd(i))*z1dd(j)*factor2* xijf_j *(idd2(i,j)+eij/t*idd4(i,j))
6696 fdd2m2(i)=fdd2m2(i) + 4.0*z1dd(j)*factor2* xijf_j *(idd2(i,j)+eij/t*idd4(i,j))
6697 IF (j == i) fdd2m2(i) =fdd2m2(i) +8.0*factor2* xijf_j*my2dd(i) *(idd2(i,j)+eij/t*idd4(i,j))
6700 xijkfa = x(i)*parame(i,3)/t*parame(i,2)**3 *x(j)*parame(j,3)/t*parame(j,2)**3 &
6701 *x(k)*parame(k,3)/t*parame(k,2)**3 / ((parame(i,2)+parame(j,2))/2.0) &
6702 /((parame(i,2)+parame(k,2))/2.0) / ((parame(j,2)+parame(k,2))/2.0) &
6703 *(z1dd(i)*z1dd(j)*z1dd(k)-z2dd(i)*z2dd(j)*z2dd(k))
6705 fdd3 = fdd3 + factor3 * xijkfa * idd3(i,j,k)
6706 xijkf_j = parame(i,3)/t*parame(i,2)**3 *x(j)*parame(j,3)/t*parame(j,2)**3 &
6707 *x(k)*parame(k,3)/t*parame(k,2)**3 /((parame(i,2)+parame(j,2))/2.0) &
6708 /((parame(i,2)+parame(k,2))/2.0) /((parame(j,2)+parame(k,2))/2.0)
6710 fdd3m(i)=fdd3m(i)+6.0*factor3*sqrt(my2dd(i))*z1dd(j)*z1dd(k) *xijkf_j*idd3(i,j,k)
6711 fdd3m2(i)=fdd3m2(i)+6.0*factor3*z1dd(j)*z1dd(k) *xijkf_j*idd3(i,j,k)
6712 IF(j == i) fdd3m2(i) =fdd3m2(i)+24.0*factor3*my2dd(i)*z1dd(k) *xijkf_j*idd3(i,j,k)
6719 IF (fdd2 < -1.e-50 .AND. fdd3 /= 0.0)
THEN 6720 fdd = fdd2 / ( 1.0 - fdd3/fdd2 )
6721 IF ( ddmax /= 0 )
THEN 6724 fddm(i) =fdd2*(fdd2*fdd2m(i) -2.0*fdd3*fdd2m(i)+fdd2*fdd3m(i)) /(fdd2-fdd3)**2
6725 fddm2(i) = fdd2m(i) * (fdd2*fdd2m(i)-2.0*fdd3*fdd2m(i) +fdd2*fdd3m(i)) / (fdd2-fdd3)**2 &
6726 + fdd2*(fdd2*fdd2m2(i) -2.0*fdd3*fdd2m2(i)+fdd2m(i)**2 &
6727 -fdd2m(i)*fdd3m(i) +fdd2*fdd3m2(i)) / (fdd2-fdd3)**2 &
6728 - 2.0*fdd2*(fdd2*fdd2m(i) -2.0*fdd3*fdd2m(i) +fdd2*fdd3m(i)) /(fdd2-fdd3)**3 &
6729 *(fdd2m(i)-fdd3m(i))
6730 dderror(i)= sqrt( my2dd(i) ) - sqrt( my0(i) ) + alph_tst(i)*fddm(i)
6731 my2dd(i) = ( sqrt( my2dd(i) ) - dderror(i) / (1.0+alph_tst(i)*fddm2(i)) )**2
6732 z1dd(i) = my2dd(i) + 3.0 * alph_tst(i)
6735 IF (abs(dderror(i)) > 1.e-11 .AND. ddit < ddmax)
GOTO 9
6737 fdd = fdd + sum( 0.5*x(1:ncomp)*alph_tst(1:ncomp)*fddm(1:ncomp)**2 )
6742 END SUBROUTINE f_dd_gross_vrabec
6750 SUBROUTINE f_qq_gross( fqq )
6753 USE eos_variables, ONLY: nc, ncomp, uij, parame, mseg, sig_ij, rho, eta, x, t
6758 REAL,
INTENT(IN OUT) :: fqq
6760 INTEGER :: i, j, k, m
6761 REAL :: factor2, factor3
6762 REAL :: xijfa, xijkfa, eij
6764 REAL,
DIMENSION(nc) :: qq2
6765 REAL,
DIMENSION(nc,nc) :: iqq2, iqq4
6766 REAL,
DIMENSION(nc,nc,nc) :: iqq3
6772 qq2(i) = (parame(i,7))**2 *1.e-69 / (uij(i,i)*kbol*mseg(i)*sig_ij(i,i)**5 *1.e-50)
6779 IF (parame(i,7) /= 0.0 .AND. parame(j,7) /= 0.0)
THEN 6781 iqq2(i,j) = iqq2(i,j) + qqp2(i,j,m)*eta**m
6782 iqq4(i,j) = iqq4(i,j) + qqp4(i,j,m)*eta**m
6786 IF (parame(k,7) /= 0.0)
THEN 6788 iqq3(i,j,k) = iqq3(i,j,k) + qqp3(i,j,k,m)*eta**m
6796 factor2 = -9.0/16.0*pi *rho
6797 factor3 = 9.0/16.0*pi**2 * rho**2
6803 IF (parame(i,7) /= 0.0 .AND. parame(j,7) /= 0.0)
THEN 6804 xijfa=x(i)*uij(i,i)*qq2(i)*sig_ij(i,i)**5 /t &
6805 *x(j)*uij(j,j)*qq2(j)*sig_ij(j,j)**5 /t/sig_ij(i,j)**7.0
6806 eij = (parame(i,3)*parame(j,3))**0.5
6807 fqq2= fqq2 +factor2* xijfa * (iqq2(i,j)+eij/t*iqq4(i,j))
6809 IF (parame(k,7) /= 0.0)
THEN 6810 xijkfa=x(i)*uij(i,i)*qq2(i)*sig_ij(i,i)**5 /t/sig_ij(i,j)**3 &
6811 *x(j)*uij(j,j)*qq2(j)*sig_ij(j,j)**5 /t/sig_ij(i,k)**3 &
6812 *x(k)*uij(k,k)*qq2(k)*sig_ij(k,k)**5 /t/sig_ij(j,k)**3
6813 fqq3 = fqq3 + factor3 * xijkfa * iqq3(i,j,k)
6820 IF ( fqq2 < -1.e-50 .AND. fqq3 /= 0.0 )
THEN 6821 fqq = fqq2 / ( 1.0 - fqq3/fqq2 )
6826 END SUBROUTINE f_qq_gross
6832 SUBROUTINE f_dq_vrabec_gross( fdq )
6835 USE eos_variables, ONLY: nc, ncomp, uij, parame, mseg, sig_ij, rho, eta, x, t
6840 REAL,
INTENT(IN OUT) :: fdq
6842 INTEGER :: i, j, k, m
6843 REAL :: factor2, factor3
6844 REAL :: xijfa, xijkfa, eij
6846 REAL,
DIMENSION(nc) :: my2dd, myfac, qq2, q_fac
6847 REAL,
DIMENSION(nc,nc) :: idq2, idq4
6848 REAL,
DIMENSION(nc,nc,nc) :: idq3
6854 my2dd(i) = (parame(i,6))**2 *1.e-49 /(uij(i,i)*kbol* mseg(i)*sig_ij(i,i)**3 *1.e-30)
6855 myfac(i) = parame(i,3)/t*parame(i,2)**4 *my2dd(i)
6857 qq2(i) = (parame(i,7))**2 *1.e-69 / (uij(i,i)*kbol*mseg(i)*sig_ij(i,i)**5 *1.e-50)
6858 q_fac(i) = parame(i,3)/t*parame(i,2)**4 *qq2(i)
6865 IF (myfac(i) /= 0.0 .AND. q_fac(j) /= 0.0)
THEN 6867 idq2(i,j) = idq2(i,j) + dqp2(i,j,m)*eta**m
6868 idq4(i,j) = idq4(i,j) + dqp4(i,j,m)*eta**m
6872 IF (myfac(k) /= 0.0 .OR. q_fac(k) /= 0.0)
THEN 6874 idq3(i,j,k) = idq3(i,j,k) + dqp3(i,j,k,m)*eta**m
6882 factor2 = -9.0/4.0 * pi *rho
6883 factor3 = pi**2 * rho**2
6889 IF (myfac(i) /= 0.0 .AND. q_fac(j) /= 0.0)
THEN 6890 xijfa = x(i)*myfac(i) * x(j)*q_fac(j) /sig_ij(i,j)**5
6891 eij = (parame(i,3)*parame(j,3))**0.5
6892 fdq2 = fdq2 +factor2* xijfa*(idq2(i,j)+eij/t*idq4(i,j))
6894 IF (myfac(k) /= 0.0 .OR. q_fac(k) /= 0.0)
THEN 6895 xijkfa=x(i)*x(j)*x(k)/(sig_ij(i,j)*sig_ij(i,k)*sig_ij(j,k))**2 &
6896 *( myfac(i)*q_fac(j)*myfac(k) + myfac(i)*q_fac(j)*q_fac(k)*1.1937350 )
6897 fdq3 = fdq3 + factor3*xijkfa*idq3(i,j,k)
6904 IF (fdq2 < -1.e-50 .AND. fdq3 /= 0.0)
THEN 6905 fdq = fdq2 / ( 1.0 - fdq3/fdq2 )
6908 END SUBROUTINE f_dq_vrabec_gross
6915 SUBROUTINE f_dft ( I1_dft, I2_dft )
6917 USE eos_variables, ONLY: nc, pi, ncomp, t, rho, eta, x, mseg, parame
6922 REAL,
INTENT(OUT) :: i1_dft
6923 REAL,
INTENT(OUT) :: i2_dft
6928 REAL :: ua, ua_c, ua_2, ua_c_2, rm
6929 REAL :: int10, int11, int20, int21
6931 REAL :: rad, xg, rdf, rho_st, msegm
6933 REAL :: dg_dr, dzr_org
6939 rho_st = rho * parame(1,2)**3
6941 ua_c = 4.0 * ( rc**(-12) - rc**(-6) )
6942 ua_c_2 = ua_c * ua_c
6946 int20 = rc*rc* ua_c_2
6951 sig_ij = parame(1,2)
6962 DO WHILE ( rad-dzr+1.e-9 >= 1.0 )
6969 ua = 4.0 * ( rad**(-12) - rad**(-6) )
6973 IF ( rad <= rg )
THEN 6974 CALL bi_cub_spline (rho_st,xg,ya,x1a,x2a,y1a,y2a,y12a, &
6975 c_bicub,rdf,dg_drho,dg_dr,den_step,ih,k)
6978 int11 = rdf*rad*rad* ua
6979 int21 = rdf*rad*rad* ua_2
6980 i1_dft= i1_dft + dzr*(int11+int10)/2.0
6981 i2_dft= i2_dft + dzr*(int21+int20)/2.0
6994 i1_dft= - i1_dft - ( 4.0/9.0 * rc**(-9) - 4.0/3.0 * rc**(-3) )
6999 i2_dft = i2_dft + 16.0/21.0 * rc**(-21) - 32.0/15.0 * rc**(-15) + 16.0/9.0 * rc**(-9)
7002 END SUBROUTINE f_dft
7011 REAL FUNCTION tangent_value2 ( optpara, n )
7019 INTEGER,
INTENT(IN) :: n
7020 REAL,
INTENT(IN) :: optpara(n)
7026 REAL :: lnphi(np,nc),ph_frac, gibbs_full(np),xlnx1,xlnx2
7027 REAL,
DIMENSION(nc) :: ni_1, ni_2
7033 IF ( optpara(i) < -300.0 )
THEN 7036 ni_2(i) = exp( optpara(i) )
7041 ni_1(i) = xif(i) - ni_2(i)
7042 IF ( ni_2(i) > xif(i) )
THEN 7044 ni_1(i) = xif(i) * 1.e-20
7048 xi(2,1:ncomp) = ni_2(1:ncomp) / sum( ni_2(1:ncomp) )
7049 lnx(2,1:ncomp) = optpara(1:ncomp) - log( sum( ni_2(1:ncomp) ) )
7051 ph_frac = sum( ni_1(1:ncomp) )
7052 xi(1,1:ncomp) = ni_1(1:ncomp) / ph_frac
7053 lnx(1,1:ncomp) = log( ni_1(1:ncomp) ) - log( ph_frac )
7058 CALL fugacity (lnphi)
7061 gibbs(1) = sum( xi(1,1:ncomp) * lnphi(1,1:ncomp) )
7062 gibbs(2) = sum( xi(2,1:ncomp) * lnphi(2,1:ncomp) )
7064 xlnx1 = sum( xi(1,1:ncomp)*lnx(1,1:ncomp) )
7065 xlnx2 = sum( xi(2,1:ncomp)*lnx(2,1:ncomp) )
7067 gibbs_full(1) = gibbs(1) + xlnx1
7068 gibbs_full(2) = gibbs(2) + xlnx2
7070 tangent_value2 = gibbs_full(1)*ph_frac + gibbs_full(2)*(1.0-ph_frac)
7077 END FUNCTION tangent_value2
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains constant...
double lambda
Latent heat of blowing agent, J/kg.
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW Module DFT_MODULE This module...
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW Module STARTING_VALUES This m...
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...