MoDeNa  1.0
Software framework facilitating sequential multi-scale modelling
eos_polar_second_deriv.f90
1 
2 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
3 !
4 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
5 
6 subroutine a_polar_drhoi_drhoj( n_comp, A_polar_rr )
7 
8  USE eos_variables, ONLY: ncomp, parame, dd_term, qq_term, dq_term
9  implicit none
10 
11  !-----------------------------------------------------------------------------
12  integer, intent(in) :: n_comp
13  real, dimension(n_comp,n_comp), intent(out) :: a_polar_rr
14 
15  !-----------------------------------------------------------------------------
16  integer :: dipole
17  integer :: quadrupole
18  integer :: dipole_quad
19  real, allocatable, dimension(:,:) :: add_rr, aqq_rr, adq_rr
20  !-----------------------------------------------------------------------------
21 
22  a_polar_rr(:,:) = 0.0
23 
24  dipole = 0
25  quadrupole = 0
26  dipole_quad = 0
27  if ( sum( parame(1:ncomp,6) ) /= 0.0 ) dipole = 1
28  if ( sum( parame(1:ncomp,7) ) /= 0.0 ) quadrupole = 1
29  if ( dipole == 1 .AND. quadrupole == 1 ) dipole_quad = 1
30 
31  !-----------------------------------------------------------------------------
32  ! dipole-dipole term
33  !-----------------------------------------------------------------------------
34  if (dipole == 1) then
35 
36  allocate( add_rr(ncomp, ncomp) )
37  if (dd_term == 'GV') CALL a_rr_dd_gross_vrabec( n_comp, add_rr )
38  if (dd_term /= 'GV' .AND. dd_term /= 'SF') write (*,*) 'specify dipole term !'
39 
40  a_polar_rr = a_polar_rr + add_rr
41  deallocate( add_rr )
42 
43  end if
44 
45  !-----------------------------------------------------------------------------
46  ! quadrupole-quadrupole term
47  !-----------------------------------------------------------------------------
48  if (quadrupole == 1) then
49 
50  allocate( aqq_rr(ncomp, ncomp) )
51  if (qq_term == 'JG') CALL a_rr_qq_gross( n_comp, aqq_rr )
52  if (qq_term /= 'JG' .AND. qq_term /= 'SF') write (*,*) 'specify quadrupole term !'
53 
54  a_polar_rr = a_polar_rr + aqq_rr
55  deallocate( aqq_rr )
56 
57  end if
58 
59  !-----------------------------------------------------------------------------
60  ! dipole-quadrupole cross term
61  !-----------------------------------------------------------------------------
62  if (dipole_quad == 1) then
63 
64  allocate( adq_rr(ncomp, ncomp) )
65  if (dq_term == 'VG') CALL a_rr_dq_vrabec_gross( n_comp, adq_rr )
66  if (dq_term /= 'VG' ) write (*,*) 'specify DQ-cross term !'
67 
68  a_polar_rr = a_polar_rr + adq_rr
69  deallocate( adq_rr )
70 
71  end if
72 
73 end subroutine a_polar_drhoi_drhoj
74 
75 
76 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
77 !
78 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
79 
80 subroutine a_rr_dd_gross_vrabec( n_comp, Add_rr )
81 
82  USE parameters, ONLY: pi, kbol
83  USE eos_variables, ONLY: nc, ncomp, uij, parame, mseg, sig_ij, rho, eta, x, t, dhs
84  USE eos_constants, ONLY: ddp2, ddp3, ddp4
85  IMPLICIT NONE
86 
87  !-----------------------------------------------------------------------------
88  integer, intent(in) :: n_comp
89  real, dimension(n_comp,n_comp), INTENT(IN OUT) :: add_rr
90 
91  !-----------------------------------------------------------------------------
92  integer :: i, j, k, l, m, o
93 
94  real :: factor2, factor3, eij
95  real :: z3, z3_m
96  real :: fdd2, fdd3
97  real, allocatable, dimension(:,:) :: xijfa
98  real, allocatable, dimension(:,:,:) :: xijkfa
99 
100  real, dimension(nc) :: my2dd, z3_r
101  real, dimension(nc) :: rhoi
102  real, allocatable, dimension(:,:) :: idd2
103  real, allocatable, dimension(:,:,:) :: idd3
104  real, allocatable, dimension(:,:,:) :: idd2_r
105  real, allocatable, dimension(:,:,:,:) :: idd3_r
106  real, allocatable, dimension(:) :: fdd2_r, fdd3_r, fdd_r
107  real, allocatable, dimension(:,:) :: idd2_rkrl
108  real, allocatable, dimension(:,:,:) :: idd3_rkrl
109  real :: fdd2_rkrl, fdd3_rkrl
110  !-----------------------------------------------------------------------------
111 
112  z3 = eta
113  do i = 1, ncomp
114  if ( uij(i,i) == 0.0 ) write (*,*) 'A_rr_DD_GROSS_VRABEC: do not use dimensionless units'
115  if ( uij(i,i) == 0.0 ) stop
116  my2dd(i) = (parame(i,6))**2 *1.e-49 / (uij(i,i)*kbol*mseg(i)*sig_ij(i,i)**3 *1.e-30)
117  z3_r(i) = pi/6.0 * mseg(i) * dhs(i)**3
118  end do
119 
120  !-----------------------------------------------------------------------------
121  ! some basic quantities
122  !-----------------------------------------------------------------------------
123 
124  allocate( xijfa(ncomp,ncomp), xijkfa(ncomp,ncomp,ncomp) )
125  allocate( idd2(ncomp,ncomp), idd2_r(ncomp,ncomp,ncomp) )
126  allocate( idd3(ncomp,ncomp,ncomp), idd3_r(ncomp,ncomp,ncomp,ncomp) )
127  allocate( fdd2_r(ncomp), fdd3_r(ncomp), fdd_r(ncomp) )
128  allocate( idd2_rkrl(ncomp,ncomp), idd3_rkrl(ncomp,ncomp,ncomp) )
129 
130  rhoi( 1:ncomp ) = x(1:ncomp ) * rho
131 
132  factor2= -pi
133  factor3= -4.0/3.0*pi*pi
134 
135  do i = 1, ncomp
136  do j = 1, ncomp
137  xijfa(i,j) = factor2* uij(i,i)*my2dd(i)*sig_ij(i,i)**3 /t &
138  *uij(j,j)*my2dd(j)*sig_ij(j,j)**3 /t / sig_ij(i,j)**3
139  do l = 1, ncomp
140  xijkfa(i,j,l)= factor3*uij(i,i)/t*my2dd(i)*sig_ij(i,i)**3 &
141  *uij(j,j)/t*my2dd(j)*sig_ij(j,j)**3 &
142  *uij(l,l)/t*my2dd(l)*sig_ij(l,l)**3 &
143  /sig_ij(i,j)/sig_ij(i,l)/sig_ij(j,l)
144  end do
145  end do
146  end do
147 
148  do i = 1, ncomp
149  do j = 1, ncomp
150  idd2(i,j) = 0.0
151  if (parame(i,6) /= 0.0 .AND. parame(j,6) /= 0.0) then
152  do m = 0, 4
153  eij = (parame(i,3)*parame(j,3))**0.5
154  idd2(i,j) =idd2(i,j) + ( ddp2(i,j,m) + eij/t*ddp4(i,j,m))*z3**m
155  end do
156  do l = 1, ncomp
157  idd3(i,j,l) = 0.0
158  if (parame(l,6) /= 0.0) then
159  do m = 0, 4
160  idd3(i,j,l) = idd3(i,j,l) + ddp3(i,j,l,m)*z3**m
161  end do
162  end if
163  end do
164  end if
165  end do
166  end do
167 
168  fdd2 = 0.0
169  fdd3 = 0.0
170  do i = 1, ncomp
171  do j = 1, ncomp
172  if (parame(i,6) /= 0.0 .AND. parame(j,6) /= 0.0) then
173  fdd2 = fdd2 + rhoi(i)*rhoi(j)*xijfa(i,j) * idd2(i,j)
174  do l=1,ncomp
175  if (parame(l,6) /= 0.0) then
176  fdd3 = fdd3 + rhoi(i)*rhoi(j)*rhoi(l)*xijkfa(i,j,l) * idd3(i,j,l)
177  end if
178  end do
179  end if
180  end do
181  end do
182 
183 
184  !-----------------------------------------------------------------------------
185  ! some first derivatives
186  !-----------------------------------------------------------------------------
187 
188  do k = 1, ncomp
189 
190  do i = 1, ncomp
191  do j = 1, ncomp
192  idd2_r(k,i,j) = 0.0
193  if (parame(i,6) /= 0.0 .AND. parame(j,6) /= 0.0) then
194  do m = 1, 4
195  z3_m = REAL(m) * z3**(m-1) * z3_r(k)
196  eij = (parame(i,3)*parame(j,3))**0.5
197  idd2_r(k,i,j) = idd2_r(k,i,j)+ ( ddp2(i,j,m) + eij/t * ddp4(i,j,m) ) * z3_m
198  end do
199  do o = 1, ncomp
200  idd3_r(k,i,j,o) = 0.0
201  if (parame(o,6) /= 0.0) then
202  do m = 1, 4
203  idd3_r(k,i,j,o)=idd3_r(k,i,j,o) + ddp3(i,j,o,m)*REAL(m)*z3**(m-1)*z3_r(k)
204  end do
205  end if
206  end do
207  end if
208  end do
209  end do
210 
211  fdd2_r(k) = 0.0
212  fdd3_r(k) = 0.0
213  do i = 1, ncomp
214  fdd2_r(k) = fdd2_r(k) + 2.0*rhoi(i)*xijfa(i,k) * idd2(i,k)
215  do j = 1, ncomp
216  if (parame(i,6) /= 0.0 .AND. parame(j,6) /= 0.0) then
217  fdd2_r(k) = fdd2_r(k) + rhoi(i)*rhoi(j)*xijfa(i,j) * idd2_r(k,i,j)
218  fdd3_r(k) = fdd3_r(k) + 3.0 * rhoi(i)*rhoi(j)*xijkfa(i,j,k) * idd3(i,j,k)
219  do o = 1, ncomp
220  if (parame(o,6) /= 0.0) then
221  fdd3_r(k) =fdd3_r(k) + rhoi(i)*rhoi(j)*rhoi(o)*xijkfa(i,j,o) * idd3_r(k,i,j,o)
222  end if
223  end do
224  end if
225  end do
226  end do
227 
228  if (fdd2 < -1.e-50 .AND. fdd3 /= 0.0 .AND. fdd2_r(k) /= 0.0 .AND. fdd3_r(k) /= 0.0) then
229 
230  fdd_r(k) = fdd2* (fdd2*fdd2_r(k) - 2.0*fdd3*fdd2_r(k) + fdd2*fdd3_r(k)) / (fdd2-fdd3)**2
231 
232  end if
233 
234  end do
235 
236  !-----------------------------------------------------------------------------
237  ! second derivatives
238  !-----------------------------------------------------------------------------
239 
240  do k = 1, ncomp
241  do l = 1, ncomp
242 
243  do i = 1, ncomp
244  do j = 1, ncomp
245  idd2_rkrl(i,j) = 0.0
246  if (parame(i,6) /= 0.0 .AND. parame(j,6) /= 0.0) then
247  do m = 2, 4
248  z3_m = REAL(m-1) * REAL(m) * z3**(m-2) * z3_r(k) * z3_r(l)
249  eij = (parame(i,3)*parame(j,3))**0.5
250  idd2_rkrl(i,j) = idd2_rkrl(i,j)+ ( ddp2(i,j,m) + eij/t * ddp4(i,j,m) ) * z3_m
251  end do
252  do o = 1, ncomp
253  idd3_rkrl(i,j,o) = 0.0
254  if (parame(o,6) /= 0.0) then
255  do m = 2, 4
256  idd3_rkrl(i,j,o)=idd3_rkrl(i,j,o) + ddp3(i,j,o,m)*REAL(m-1)*REAL(m)*z3**(m-2)*z3_r(k)*z3_r(l)
257  end do
258  end if
259  end do
260  end if
261  end do
262  end do
263 
264  fdd2_rkrl = 2.0 * xijfa(k,l) * idd2(k,l)
265  fdd3_rkrl = 0.0
266  do i = 1, ncomp
267  fdd2_rkrl = fdd2_rkrl + 2.0 * rhoi(i) * ( xijfa(i,k) * idd2_r(l,i,k) + xijfa(i,l) * idd2_r(k,i,l))
268  fdd3_rkrl = fdd3_rkrl + 6.0 * rhoi(i) * xijkfa(i,k,l) * idd3(i,k,l)
269  do j = 1, ncomp
270  if (parame(i,6) /= 0.0 .AND. parame(j,6) /= 0.0) then
271  fdd2_rkrl = fdd2_rkrl + rhoi(i)*rhoi(j)*xijfa(i,j) * idd2_rkrl(i,j)
272  fdd3_rkrl = fdd3_rkrl + 3.0 * rhoi(i)*rhoi(j)* ( xijkfa(i,j,k) * idd3_r(l,i,j,k) &
273  + xijkfa(i,j,l) * idd3_r(k,i,j,l) )
274  do o = 1, ncomp
275  if (parame(o,6) /= 0.0) then
276  fdd3_rkrl =fdd3_rkrl + rhoi(i)*rhoi(j)*rhoi(o)*xijkfa(i,j,o) * idd3_rkrl(i,j,o)
277  end if
278  end do
279  end if
280  end do
281  end do
282 
283 
284  add_rr(k,l) = ( 2.0*fdd2*fdd2_r(l)*( fdd2_r(k) + fdd3_r(k) ) + fdd2*fdd2*( fdd2_rkrl + fdd3_rkrl ) &
285  -2.0*( fdd2_r(k)*fdd2_r(l)*fdd3 + fdd2*fdd3_r(l)*fdd2_r(k) + fdd2*fdd3*fdd2_rkrl ) ) &
286  / ( fdd2 - fdd3 )**2 &
287  + fdd_r(k) * 2.0 * ( fdd3_r(l) - fdd2_r(l) ) / ( fdd2 - fdd3 )
288 
289  end do
290  end do
291 
292  deallocate( xijfa, xijkfa, idd2, idd2_r, idd3, idd3_r, fdd2_r, fdd3_r, fdd_r, idd2_rkrl, idd3_rkrl )
293 
294 end subroutine a_rr_dd_gross_vrabec
295 
296 
297 
298 
299 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
300 !
301 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
302 
303 subroutine a_rr_qq_gross( n_comp, Aqq_rr )
304 
305  USE parameters, ONLY: pi, kbol
306  USE eos_variables, ONLY: nc, ncomp, uij, parame, mseg, sig_ij, rho, eta, x, t, dhs
307  USE eos_constants, ONLY: qqp2, qqp3, qqp4
308  IMPLICIT NONE
309 
310  !-----------------------------------------------------------------------------
311  integer, intent(in) :: n_comp
312  real, dimension(n_comp,n_comp), INTENT(IN OUT) :: aqq_rr
313 
314  !-----------------------------------------------------------------------------
315  integer :: i, j, k, l, m, o
316 
317  real :: factor2, factor3, eij
318  real :: z3, z3_m
319  real :: fqq2, fqq3
320  real, allocatable, dimension(:,:) :: xijfa
321  real, allocatable, dimension(:,:,:) :: xijkfa
322 
323  real, dimension(nc) :: qq2, z3_r
324  real, dimension(nc) :: rhoi
325  real, allocatable, dimension(:,:) :: iqq2
326  real, allocatable, dimension(:,:,:) :: iqq3
327  real, allocatable, dimension(:,:,:) :: iqq2_r
328  real, allocatable, dimension(:,:,:,:) :: iqq3_r
329  real, allocatable, dimension(:) :: fqq2_r, fqq3_r, fqq_r
330  real, allocatable, dimension(:,:) :: iqq2_rkrl
331  real, allocatable, dimension(:,:,:) :: iqq3_rkrl
332  real :: fqq2_rkrl, fqq3_rkrl
333  !-----------------------------------------------------------------------------
334 
335  z3 = eta
336  do i = 1, ncomp
337  if ( uij(i,i) == 0.0 ) write (*,*) 'A_rr_QQ_GROSS: do not use dimensionless units'
338  if ( uij(i,i) == 0.0 ) stop
339  qq2(i) = (parame(i,7))**2 *1.e-69 / (uij(i,i)*kbol*mseg(i)*sig_ij(i,i)**5 *1.e-50)
340  z3_r(i) = pi/6.0 * mseg(i) * dhs(i)**3
341  end do
342 
343  !-----------------------------------------------------------------------------
344  ! some basic quantities
345  !-----------------------------------------------------------------------------
346 
347  allocate( xijfa(ncomp,ncomp), xijkfa(ncomp,ncomp,ncomp) )
348  allocate( iqq2(ncomp,ncomp), iqq2_r(ncomp,ncomp,ncomp) )
349  allocate( iqq3(ncomp,ncomp,ncomp), iqq3_r(ncomp,ncomp,ncomp,ncomp) )
350  allocate( fqq2_r(ncomp), fqq3_r(ncomp), fqq_r(ncomp) )
351  allocate( iqq2_rkrl(ncomp,ncomp), iqq3_rkrl(ncomp,ncomp,ncomp) )
352 
353  rhoi( 1:ncomp ) = x(1:ncomp ) * rho
354 
355  factor2= -9.0/16.0*pi
356  factor3= 9.0/16.0*pi**2
357 
358  do i = 1, ncomp
359  do j = 1, ncomp
360  xijfa(i,j) = factor2* uij(i,i)*qq2(i)*sig_ij(i,i)**5 /t &
361  *uij(j,j)*qq2(j)*sig_ij(j,j)**5 /t / sig_ij(i,j)**7
362  do l = 1, ncomp
363  xijkfa(i,j,l)= factor3*uij(i,i)/t*qq2(i)*sig_ij(i,i)**5 &
364  *uij(j,j)/t*qq2(j)*sig_ij(j,j)**5 &
365  *uij(l,l)/t*qq2(l)*sig_ij(l,l)**5 &
366  / ( sig_ij(i,j)*sig_ij(i,l)*sig_ij(j,l) )**3
367  end do
368  end do
369  end do
370 
371  do i = 1, ncomp
372  do j = 1, ncomp
373  iqq2(i,j) = 0.0
374  if (parame(i,7) /= 0.0 .AND. parame(j,7) /= 0.0) then
375  do m = 0, 4
376  eij = (parame(i,3)*parame(j,3))**0.5
377  iqq2(i,j) = iqq2(i,j) + ( qqp2(i,j,m) + eij/t * qqp4(i,j,m) ) * z3**m
378  end do
379  do l = 1, ncomp
380  iqq3(i,j,l) = 0.0
381  if (parame(l,7) /= 0.0) then
382  do m = 0, 4
383  iqq3(i,j,l) = iqq3(i,j,l) + qqp3(i,j,l,m)*z3**m
384  end do
385  end if
386  end do
387  end if
388  end do
389  end do
390 
391  fqq2 = 0.0
392  fqq3 = 0.0
393  do i = 1, ncomp
394  do j = 1, ncomp
395  if (parame(i,7) /= 0.0 .AND. parame(j,7) /= 0.0) then
396  fqq2 = fqq2 + rhoi(i)*rhoi(j)*xijfa(i,j) * iqq2(i,j)
397  do l=1,ncomp
398  if (parame(l,7) /= 0.0) then
399  fqq3 = fqq3 + rhoi(i)*rhoi(j)*rhoi(l)*xijkfa(i,j,l) * iqq3(i,j,l)
400  end if
401  end do
402  end if
403  end do
404  end do
405 
406 
407  !-----------------------------------------------------------------------------
408  ! some first derivatives
409  !-----------------------------------------------------------------------------
410 
411  do k = 1, ncomp
412 
413  do i = 1, ncomp
414  do j = 1, ncomp
415  iqq2_r(k,i,j) = 0.0
416  if (parame(i,7) /= 0.0 .AND. parame(j,7) /= 0.0) then
417  do m = 1, 4
418  z3_m = REAL(m) * z3**(m-1) * z3_r(k)
419  eij = (parame(i,3)*parame(j,3))**0.5
420  iqq2_r(k,i,j) = iqq2_r(k,i,j)+ ( qqp2(i,j,m) + eij/t * qqp4(i,j,m) ) * z3_m
421  end do
422  do o = 1, ncomp
423  iqq3_r(k,i,j,o) = 0.0
424  if (parame(o,7) /= 0.0) then
425  do m = 1, 4
426  iqq3_r(k,i,j,o)=iqq3_r(k,i,j,o) + qqp3(i,j,o,m)*REAL(m)*z3**(m-1)*z3_r(k)
427  end do
428  end if
429  end do
430  end if
431  end do
432  end do
433 
434  fqq2_r(k) = 0.0
435  fqq3_r(k) = 0.0
436  do i = 1, ncomp
437  fqq2_r(k) = fqq2_r(k) + 2.0*rhoi(i)*xijfa(i,k) * iqq2(i,k)
438  do j = 1, ncomp
439  if (parame(i,7) /= 0.0 .AND. parame(j,7) /= 0.0) then
440  fqq2_r(k) = fqq2_r(k) + rhoi(i)*rhoi(j)*xijfa(i,j) * iqq2_r(k,i,j)
441  fqq3_r(k) = fqq3_r(k) + 3.0 * rhoi(i)*rhoi(j)*xijkfa(i,j,k) * iqq3(i,j,k)
442  do o = 1, ncomp
443  if (parame(o,7) /= 0.0) then
444  fqq3_r(k) =fqq3_r(k) + rhoi(i)*rhoi(j)*rhoi(o)*xijkfa(i,j,o) * iqq3_r(k,i,j,o)
445  end if
446  end do
447  end if
448  end do
449  end do
450 
451  if (fqq2 < -1.e-50 .AND. fqq3 /= 0.0 .AND. fqq2_r(k) /= 0.0 .AND. fqq3_r(k) /= 0.0) then
452 
453  fqq_r(k) = fqq2* (fqq2*fqq2_r(k) - 2.0*fqq3*fqq2_r(k) + fqq2*fqq3_r(k)) / (fqq2-fqq3)**2
454 
455  end if
456 
457  end do
458 
459  !-----------------------------------------------------------------------------
460  ! second derivatives
461  !-----------------------------------------------------------------------------
462 
463  do k = 1, ncomp
464  do l = 1, ncomp
465 
466  do i = 1, ncomp
467  do j = 1, ncomp
468  iqq2_rkrl(i,j) = 0.0
469  if (parame(i,7) /= 0.0 .AND. parame(j,7) /= 0.0) then
470  do m = 2, 4
471  z3_m = REAL(m-1) * REAL(m) * z3**(m-2) * z3_r(k) * z3_r(l)
472  eij = (parame(i,3)*parame(j,3))**0.5
473  iqq2_rkrl(i,j) = iqq2_rkrl(i,j)+ ( qqp2(i,j,m) + eij/t * qqp4(i,j,m) ) * z3_m
474  end do
475  do o = 1, ncomp
476  iqq3_rkrl(i,j,o) = 0.0
477  if (parame(o,7) /= 0.0) then
478  do m = 2, 4
479  iqq3_rkrl(i,j,o)=iqq3_rkrl(i,j,o) + qqp3(i,j,o,m)*REAL(m-1)*REAL(m)*z3**(m-2)*z3_r(k)*z3_r(l)
480  end do
481  end if
482  end do
483  end if
484  end do
485  end do
486 
487  fqq2_rkrl = 2.0 * xijfa(k,l) * iqq2(k,l)
488  fqq3_rkrl = 0.0
489  do i = 1, ncomp
490  fqq2_rkrl = fqq2_rkrl + 2.0 * rhoi(i) * ( xijfa(i,k) * iqq2_r(l,i,k) + xijfa(i,l) * iqq2_r(k,i,l))
491  fqq3_rkrl = fqq3_rkrl + 6.0 * rhoi(i) * xijkfa(i,k,l) * iqq3(i,k,l)
492  do j = 1, ncomp
493  if (parame(i,7) /= 0.0 .AND. parame(j,7) /= 0.0) then
494  fqq2_rkrl = fqq2_rkrl + rhoi(i)*rhoi(j)*xijfa(i,j) * iqq2_rkrl(i,j)
495  fqq3_rkrl = fqq3_rkrl + 3.0 * rhoi(i)*rhoi(j)* ( xijkfa(i,j,k) * iqq3_r(l,i,j,k) &
496  + xijkfa(i,j,l) * iqq3_r(k,i,j,l) )
497  do o = 1, ncomp
498  if (parame(o,7) /= 0.0) then
499  fqq3_rkrl =fqq3_rkrl + rhoi(i)*rhoi(j)*rhoi(o)*xijkfa(i,j,o) * iqq3_rkrl(i,j,o)
500  end if
501  end do
502  end if
503  end do
504  end do
505 
506 
507  aqq_rr(k,l) = ( 2.0*fqq2*fqq2_r(l)*( fqq2_r(k) + fqq3_r(k) ) + fqq2*fqq2*( fqq2_rkrl + fqq3_rkrl ) &
508  -2.0*( fqq2_r(k)*fqq2_r(l)*fqq3 + fqq2*fqq3_r(l)*fqq2_r(k) + fqq2*fqq3*fqq2_rkrl ) ) &
509  / ( fqq2 - fqq3 )**2 &
510  + fqq_r(k) * 2.0 * ( fqq3_r(l) - fqq2_r(l) ) / ( fqq2 - fqq3 )
511 
512  end do
513  end do
514 
515  deallocate( xijfa, xijkfa, iqq2, iqq2_r, iqq3, iqq3_r, fqq2_r, fqq3_r, fqq_r, iqq2_rkrl, iqq3_rkrl )
516 
517 end subroutine a_rr_qq_gross
518 
519 !!$ do i=1,ncomp
520 !!$ my2dd(i) = (parame(i,6))**2 *1.E-49 / (uij(i,i)*KBOL*mseg(i)*sig_ij(i,i)**3 *1.E-30)
521 !!$ myfac(i) = parame(i,3)/t*parame(i,2)**4 *my2dd(i)
522 !!$ qq2(i) = (parame(i,7))**2 *1.E-69 / (uij(i,i)*kbol*mseg(i)*sig_ij(i,i)**5 *1.E-50)
523 !!$ q_fac(i) = parame(i,3)/t*parame(i,2)**4 *qq2(i)
524 !!$ end do
525 
526 
527 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
528 !
529 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
530 
531 subroutine a_rr_dq_vrabec_gross( n_comp, Adq_rr )
532 
533  USE parameters, ONLY: pi, kbol
534  USE eos_variables, ONLY: nc, ncomp, uij, parame, mseg, sig_ij, rho, eta, x, t, dhs
535  USE eos_constants, ONLY: dqp2, dqp3, dqp4
536  IMPLICIT NONE
537 
538  !-----------------------------------------------------------------------------
539  integer, intent(in) :: n_comp
540  real, dimension(n_comp,n_comp), INTENT(IN OUT) :: adq_rr
541 
542  !-----------------------------------------------------------------------------
543  integer :: i, j, k, l, m, o
544 
545  real :: factor2, factor3, eij
546  real :: z3, z3_m
547  real :: fdq2, fdq3
548  real, allocatable, dimension(:,:) :: xijfa
549  real, allocatable, dimension(:,:,:) :: xijkfa
550  !real :: xijkf_rk
551 
552  real, dimension(nc) :: z3_r
553  real, dimension(nc) :: my2dd, myfac, qq2, q_fac
554  real, dimension(nc) :: rhoi
555  real, allocatable, dimension(:,:) :: idq2
556  real, allocatable, dimension(:,:,:) :: idq3
557  real, allocatable, dimension(:,:,:) :: idq2_r
558  real, allocatable, dimension(:,:,:,:) :: idq3_r
559  real, allocatable, dimension(:) :: fdq2_r, fdq3_r, fdq_r
560  real, allocatable, dimension(:,:) :: idq2_rkrl
561  real, allocatable, dimension(:,:,:) :: idq3_rkrl
562  real :: fdq2_rkrl, fdq3_rkrl
563  !-----------------------------------------------------------------------------
564 
565  z3 = eta
566  do i = 1, ncomp
567  if ( uij(i,i) == 0.0 ) write (*,*) 'A_rr_DQ_GROSS: do not use dimensionless units'
568  if ( uij(i,i) == 0.0 ) stop
569  my2dd(i) = (parame(i,6))**2 *1.e-49 / (uij(i,i)*kbol*mseg(i)*sig_ij(i,i)**3 *1.e-30)
570  myfac(i) = parame(i,3)/t*parame(i,2)**4 *my2dd(i)
571  qq2(i) = (parame(i,7))**2 *1.e-69 / (uij(i,i)*kbol*mseg(i)*sig_ij(i,i)**5 *1.e-50)
572  q_fac(i) = parame(i,3)/t*parame(i,2)**4 *qq2(i)
573  z3_r(i) = pi/6.0 * mseg(i) * dhs(i)**3
574  end do
575 
576  !-----------------------------------------------------------------------------
577  ! some basic quantities
578  !-----------------------------------------------------------------------------
579 
580  allocate( xijfa(ncomp,ncomp), xijkfa(ncomp,ncomp,ncomp) )
581  allocate( idq2(ncomp,ncomp), idq2_r(ncomp,ncomp,ncomp) )
582  allocate( idq3(ncomp,ncomp,ncomp), idq3_r(ncomp,ncomp,ncomp,ncomp) )
583  allocate( fdq2_r(ncomp), fdq3_r(ncomp), fdq_r(ncomp) )
584  allocate( idq2_rkrl(ncomp,ncomp), idq3_rkrl(ncomp,ncomp,ncomp) )
585 
586  rhoi( 1:ncomp ) = x(1:ncomp ) * rho
587 
588  factor2 = -9.0/4.0*pi
589  factor3 = pi**2
590 
591  do i = 1, ncomp
592  do j = 1, ncomp
593  xijfa(i,j) = factor2* myfac(i) * q_fac(j) /sig_ij(i,j)**5 ! caution: asymmetric ...(i,j) /= ...(j,i)
594  do l = 1, ncomp
595  xijkfa(i,j,l)= factor3*( myfac(i)*q_fac(j)*myfac(l) + myfac(i)*q_fac(j)*q_fac(l)*1.193735 ) &
596  / (sig_ij(i,j)*sig_ij(i,l)*sig_ij(j,l))**2
597  end do
598  end do
599  end do
600 
601  do i = 1, ncomp
602  do j = 1, ncomp
603  idq2(i,j) = 0.0
604  do m = 0, 4
605  eij = (parame(i,3)*parame(j,3))**0.5
606  idq2(i,j) = idq2(i,j) + ( dqp2(i,j,m) + eij/t * dqp4(i,j,m) ) * z3**m
607  end do
608  do l = 1, ncomp
609  idq3(i,j,l) = 0.0
610  do m = 0, 4
611  idq3(i,j,l) = idq3(i,j,l) + dqp3(i,j,l,m)*z3**m
612  end do
613  end do
614  end do
615  end do
616 
617  fdq2 = 0.0
618  fdq3 = 0.0
619  do i = 1, ncomp
620  do j = 1, ncomp
621  fdq2 = fdq2 + rhoi(i)*rhoi(j)*xijfa(i,j) * idq2(i,j)
622  do l=1,ncomp
623  fdq3 = fdq3 + rhoi(i)*rhoi(j)*rhoi(l)*xijkfa(i,j,l) * idq3(i,j,l)
624  end do
625  end do
626  end do
627 
628 
629  !-----------------------------------------------------------------------------
630  ! some first derivatives
631  !-----------------------------------------------------------------------------
632 
633  do k = 1, ncomp
634 
635  do i = 1, ncomp
636  do j = 1, ncomp
637  idq2_r(k,i,j) = 0.0
638  do m = 1, 4
639  z3_m = REAL(m) * z3**(m-1) * z3_r(k)
640  eij = (parame(i,3)*parame(j,3))**0.5
641  idq2_r(k,i,j) = idq2_r(k,i,j)+ ( dqp2(i,j,m) + eij/t * dqp4(i,j,m) ) * z3_m
642  end do
643  do o = 1, ncomp
644  idq3_r(k,i,j,o) = 0.0
645  do m = 1, 4
646  idq3_r(k,i,j,o)=idq3_r(k,i,j,o) + dqp3(i,j,o,m)*REAL(m)*z3**(m-1)*z3_r(k)
647  end do
648  end do
649  end do
650  end do
651 
652  fdq2_r(k) = 0.0
653  fdq3_r(k) = 0.0
654  do i = 1, ncomp
655  !xijf_r(k,i) = factor2 * ( myfac(k)*q_fac(i) + myfac(i)*q_fac(k) ) / sig_ij(i,k)**5
656  !fdq2_r(k) = fdq2_r(k) +rhoi(i)*xijf_r(k,i) *Idq2(i,k)
657  fdq2_r(k) = fdq2_r(k) +rhoi(i)*( xijfa(i,k)+xijfa(k,i) ) *idq2(i,k)
658  do j = 1, ncomp
659  fdq2_r(k) = fdq2_r(k) + rhoi(i)*rhoi(j)*xijfa(i,j) * idq2_r(k,i,j)
660  !xijkf_rk=rhoi(i)*rhoi(j) &
661  ! *( myfac(i)*q_fac(j)*myfac(k) + myfac(i)*q_fac(k)*myfac(j) + myfac(k)*q_fac(i)*myfac(j) &
662  ! + (myfac(i)*q_fac(j)*q_fac(k) + myfac(i)*q_fac(k)*q_fac(j) &
663  ! + myfac(k)*q_fac(i)*q_fac(j))*1.193735 ) / (sig_ij(i,j)*sig_ij(i,k)*sig_ij(j,k))**2
664  !fdq3_r(k) = fdq3_r(k) + factor3 * xijkf_rk * Idq3(i,j,k)
665  fdq3_r(k) = fdq3_r(k) + rhoi(i)*rhoi(j)* (xijkfa(i,j,k)+xijkfa(i,k,j)+xijkfa(k,i,j)) *idq3(i,j,k)
666  do o = 1, ncomp
667  fdq3_r(k) =fdq3_r(k) + rhoi(i)*rhoi(j)*rhoi(o)*xijkfa(i,j,o) * idq3_r(k,i,j,o)
668  end do
669  end do
670  end do
671 
672  if (fdq2 < -1.e-50 .AND. fdq3 /= 0.0 .AND. fdq2_r(k) /= 0.0 .AND. fdq3_r(k) /= 0.0) then
673 
674  fdq_r(k) = fdq2* (fdq2*fdq2_r(k) - 2.0*fdq3*fdq2_r(k) + fdq2*fdq3_r(k)) / (fdq2-fdq3)**2
675 
676  end if
677 
678  end do
679 
680  !-----------------------------------------------------------------------------
681  ! second derivatives
682  !-----------------------------------------------------------------------------
683 
684  do k = 1, ncomp
685  do l = 1, ncomp
686 
687  do i = 1, ncomp
688  do j = 1, ncomp
689  idq2_rkrl(i,j) = 0.0
690  do m = 2, 4
691  z3_m = REAL(m-1) * REAL(m) * z3**(m-2) * z3_r(k) * z3_r(l)
692  eij = (parame(i,3)*parame(j,3))**0.5
693  idq2_rkrl(i,j) = idq2_rkrl(i,j)+ ( dqp2(i,j,m) + eij/t * dqp4(i,j,m) ) * z3_m
694  end do
695  do o = 1, ncomp
696  idq3_rkrl(i,j,o) = 0.0
697  do m = 2, 4
698  idq3_rkrl(i,j,o)=idq3_rkrl(i,j,o) + dqp3(i,j,o,m)*REAL(m-1)*REAL(m)*z3**(m-2)*z3_r(k)*z3_r(l)
699  end do
700  end do
701  end do
702  end do
703 
704 
705  !fdq2_rkrl = factor2 * ( myfac(k)*q_fac(l) + myfac(l)*q_fac(k) ) /sig_ij(k,l)**5 * Idq2(k,l)
706  fdq2_rkrl = ( xijfa(l,k)+xijfa(k,l) ) * idq2(k,l)
707  fdq3_rkrl = 0.0
708  do i = 1, ncomp
709  ! fdq2_rkrl = fdq2_rkrl + rhoi(i) * ( xijf_r(k,i) * Idq2_r(l,i,k) + xijf_r(l,i) * Idq2_r(k,i,l))
710  fdq2_rkrl = fdq2_rkrl + rhoi(i) * ( ( xijfa(i,k)+xijfa(k,i) ) * idq2_r(l,i,k) &
711  + ( xijfa(i,l)+xijfa(l,i) ) * idq2_r(k,i,l) )
712  fdq3_rkrl = fdq3_rkrl + rhoi(i) * ( xijkfa(i,l,k)+xijkfa(i,k,l)+xijkfa(l,i,k) &
713  +xijkfa(k,i,l)+xijkfa(l,k,i)+xijkfa(k,l,i) ) * idq3(i,k,l)
714  do j = 1, ncomp
715  fdq2_rkrl = fdq2_rkrl + rhoi(i)*rhoi(j)*xijfa(i,j) * idq2_rkrl(i,j)
716  fdq3_rkrl = fdq3_rkrl + rhoi(i)*rhoi(j)*( (xijkfa(i,j,k)+xijkfa(i,k,j)+xijkfa(k,i,j)) *idq3_r(l,i,j,k) &
717  + (xijkfa(i,j,l)+xijkfa(i,l,j)+xijkfa(l,i,j)) *idq3_r(k,i,j,l) )
718  do o = 1, ncomp
719  fdq3_rkrl =fdq3_rkrl + rhoi(i)*rhoi(j)*rhoi(o)*xijkfa(i,j,o) * idq3_rkrl(i,j,o)
720  end do
721  end do
722  end do
723 
724  adq_rr(k,l) = ( 2.0*fdq2*fdq2_r(l)*( fdq2_r(k) + fdq3_r(k) ) + fdq2*fdq2*( fdq2_rkrl + fdq3_rkrl ) &
725  -2.0*( fdq2_r(k)*fdq2_r(l)*fdq3 + fdq2*fdq3_r(l)*fdq2_r(k) + fdq2*fdq3*fdq2_rkrl ) ) &
726  / ( fdq2 - fdq3 )**2 &
727  + fdq_r(k) * 2.0 * ( fdq3_r(l) - fdq2_r(l) ) / ( fdq2 - fdq3 )
728 
729  end do
730  end do
731 
732  deallocate( xijfa, xijkfa, idq2, idq2_r, idq3, idq3_r, fdq2_r, fdq3_r, fdq_r, idq2_rkrl, idq3_rkrl )
733 
734 end subroutine a_rr_dq_vrabec_gross
parameters
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains constant...
Definition: modules.f90:6
eos_constants
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:200
eos_variables
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:120