1C> \ingroup nwxc
2C> @{
3C>
4C> \file nwxcm_x_sogga.F
5C> The nwxcm_x_sogga functional
6C>
7C> @}
8C>
9C> \ingroup nwxc_priv
10C> @{
11C>
12C> \brief Evaluate the nwxcm_x_sogga functional [1]
13C>
14C> \f{eqnarray*}{
15C>   {\it t_1} &=& {\it param}\left(1\right)\\\\
16C>   {\it t_2} &=& {\it param}\left(7\right)\\\\
17C>   {\it t_3} &=& {\it param}\left(2\right)\\\\
18C>   {\it t_4} &=& {{0.003680288926083876\,
19C>    \sigma_{\alpha\alpha}}\over{\rho_\alpha^{{{8}\over{3}}}}}\\\\
20C>   {\it t_5} &=& 1.0-{{1.0}\over{{\it t_4}+1.0}}\\\\
21C>   {\it t_6} &=& {\it param}\left(3\right)\\\\
22C>   {\it t_7} &=& {\it param}\left(4\right)\\\\
23C>   {\it t_8} &=& {\it param}\left(5\right)\\\\
24C>   {\it t_9} &=& {\it param}\left(6\right)\\\\
25C>   {\it t_{10}} &=& {\it param}\left(9\right)\\\\
26C>   {\it t_{11}} &=& 1.0-{{1}\over{e^{{\it t_4}}}}\\\\
27C>   {\it t_{12}} &=& {\it param}\left(10\right)\\\\
28C>   {\it t_{13}} &=& {\it param}\left(11\right)\\\\
29C>   {\it t_{14}} &=& {\it param}\left(12\right)\\\\
30C>   {\it t_{15}} &=& {\it param}\left(8\right)\\\\
31C>   {\it t_{16}} &=& {{0.003680288926083876\,
32C>    \sigma_{\beta\beta}}\over{\rho_\beta^{{{8}\over{3}}}}}\\\\
33C>   {\it t_{17}} &=& 1.0-{{1.0}\over{{\it t_{16}}+1.0}}\\\\
34C>   {\it t_{18}} &=& 1.0-{{1}\over{e^{{\it t_{16}}}}}\\\\
35C>   {\it t_{19}} &=& {{0.003680288926083876\,\sigma_{ss}}
36C>    \over{\rho_s^{{{8}\over{3}}}}}\\\\
37C>   {\it t_{20}} &=& 1.0-{{1.0}\over{{\it t_{19}}+1.0}}\\\\
38C>   {\it t_{21}} &=& 1.0-{{1}\over{e^{{\it t_{19}}}}}\\\\
39C>   f &=& -0.9305257363491002\,\rho_\beta^{{{4}\over{3}}}\,
40C>    \left({\it t_{15}}\,{\it t_{18}}+{\it t_{14}}\,{
41C>    \it t_{18}}^{5.0}+{\it t_{13}}\,{\it t_{18}}^{4.0}+{
42C>    \it t_{12}}\,{\it t_{18}}^{3.0}+{\it t_{10}}\,{\it t_{18}}^{2.0}
43C>    +{\it t_9}\,{\it t_{17}}^{5.0}+{\it t_8}\,{\it t_{17}}^{4.0}
44C>    +{\it t_7}\,{\it t_{17}}^{3.0}+{\it t_6}\,{\it t_{17}}^{2.0}
45C>    +{\it t_3}\,{\it t_{17}}+{\it t_2}+{\it t_1}\right)
46C>    -0.9305257363491002\,\rho_\alpha^{{{4}\over{3}}}\,\left({
47C>    \it t_{15}}\,{\it t_{11}}+{\it t_{14}}\,{\it t_{11}}^{5.0}+{
48C>    \it t_{13}}\,{\it t_{11}}^{4.0}+{\it t_{12}}\,{\it t_{11}}^{3.0}
49C>    +{\it t_{10}}\,{\it t_{11}}^{2.0}+{\it t_9}\,{\it t_5}^{5.0}
50C>    +{\it t_8}\,{\it t_5}^{4.0}+{\it t_7}\,{\it t_5}^{3.0}+{
51C>    \it t_6}\,{\it t_5}^{2.0}+{\it t_3}\,{\it t_5}+{\it t_2}+{
52C>    \it t_1}\right)\\\\
53C>   g &=& 0\\\\
54C>   G &=& -0.9305257363491002\,\rho_s^{{{4}\over{3}}}\,\left({
55C>    \it t_{15}}\,{\it t_{21}}+{\it t_{14}}\,{\it t_{21}}^{5.0}+{
56C>    \it t_{13}}\,{\it t_{21}}^{4.0}+{\it t_{12}}\,{\it t_{21}}^{3.0}
57C>    +{\it t_{10}}\,{\it t_{21}}^{2.0}+{\it t_9}\,{\it t_{20}}^{5.0}
58C>    +{\it t_8}\,{\it t_{20}}^{4.0}+{\it t_7}\,{\it t_{20}}^{3.0}
59C>    +{\it t_6}\,{\it t_{20}}^{2.0}+{\it t_3}\,{\it t_{20}}+{\it t_2}
60C>    +{\it t_1}\right)\\\\
61C> \f}
62C>
63C> Code generated with Maxima 5.34.0 [2]
64C> driven by autoxc [3].
65C>
66C> ### References ###
67C>
68C> [1] Y Zhao, DG Truhlar, J.Chem.Phys. 128, 184109 (2008)  , DOI:
69C> <a href="https://doi.org/10.1063/1.2912068 ">
70C> 10.1063/1.2912068 </a>
71C>
72C> [2] Maxima, a computer algebra system,
73C> <a href="http://maxima.sourceforge.net/">
74C> http://maxima.sourceforge.net/</a>
75C>
76C> [3] autoxc, revision 27097 2015-05-08
77C>
78      subroutine nwxcm_x_sogga(param,tol_rho,ipol,nq,wght,
79     +rho,rgamma,fnc,Amat,Cmat)
80c $Id: $
81#ifdef NWXC_QUAD_PREC
82      implicit real(kind=selected_real_kind(30))(a-h,o-z),integer(i-n)
83      integer, parameter :: rk=selected_real_kind(30)
84#else
85      implicit real(kind=selected_real_kind(15))(a-h,o-z),integer(i-n)
86      integer, parameter :: rk=selected_real_kind(15)
87#endif
88#include "nwxc_param.fh"
89      double precision param(*)     !< [Input] Parameters of functional
90      double precision tol_rho      !< [Input] The lower limit on the density
91      integer ipol                  !< [Input] The number of spin channels
92      integer nq                    !< [Input] The number of points
93      double precision wght         !< [Input] The weight of the functional
94      double precision rho(nq,*)    !< [Input] The density
95      double precision rgamma(nq,*) !< [Input] The norm of the density
96                                    !< gradients
97      double precision fnc(nq)      !< [Output] The value of the functional
98c
99c     Sampling Matrices for the XC Kernel
100c
101      double precision Amat(nq,*)   !< [Output] The derivative wrt rho
102      double precision Cmat(nq,*)   !< [Output] The derivative wrt rgamma
103      integer iq
104      double precision tmp
105      double precision rhoa,rhob
106      double precision gammaaa,gammaab,gammabb
107      double precision taua,taub
108      double precision nwxcm_heaviside
109      external         nwxcm_heaviside
110CDIR$ NOVECTOR
111      do iq = 1, nq
112        if (ipol.eq.1) then
113          rhoa    = 0.5d0*rho(iq,R_T)
114          gammaaa = 0.25d0*rgamma(iq,G_TT)
115          if (rhoa.gt.tol_rho) then
116            t1 = rhoa**1.3333333333333333d+0
117            t2 = param(3)
118            t3 = 1/rhoa**2.6666666666666666d+0
119            t4 = 3.6802889260838756d-3*gammaaa*t3
120            t5 = t4+1.0d+0
121            t6 = 1.0d+0-1.0d+0/t5
122            t7 = t6**2.0d+0
123            t8 = param(4)
124            t9 = t6**3.0d+0
125            t10 = param(5)
126            t11 = t6**4.0d+0
127            t12 = param(6)
128            t13 = param(2)
129            t14 = param(9)
130            t15 = exp(-t4)
131            t16 = 1.0d+0-t15
132            t17 = t16**2.0d+0
133            t18 = param(10)
134            t19 = t16**3.0d+0
135            t20 = param(11)
136            t21 = t16**4.0d+0
137            t22 = param(12)
138            t23 = param(8)
139            t24 = t8*t9+t2*t7+t12*t6**5.0d+0+t13*t6+t16*t23+t16**5.0d+0*
140     1         t22+t20*t21+t18*t19+t14*t17+t10*t11+param(7)+param(1)
141            t25 = 1/t5**2
142            t26 = 1/rhoa**3.6666666666666664d+0
143            fnc(iq) = fnc(iq)-1.8610514726982003d+0*t1*t24*wght
144            Amat(iq,D1_RA) = (-9.305257363491002d-1*t1*(-3.9256415211561
145     1         34d-2*gammaaa*t10*t25*t26*t9-2.9442311408671007d-2*gammaa
146     2         a*t25*t26*t7*t8-1.962820760578067d-2*gammaaa*t2*t25*t26*t
147     3         6-9.814103802890335d-3*gammaaa*t13*t25*t26-4.907051901445
148     4         1675d-2*gammaaa*t11*t12*t25*t26-9.814103802890335d-3*gamm
149     5         aaa*t15*t23*t26-4.9070519014451675d-2*gammaaa*t15*t21*t22
150     6         *t26-3.925641521156134d-2*gammaaa*t15*t19*t20*t26-2.94423
151     7         11408671007d-2*gammaaa*t15*t17*t18*t26-1.962820760578067d
152     8         -2*gammaaa*t14*t15*t16*t26)-1.2407009817988002d+0*rhoa**3
153     9         .333333333333333d-1*t24)*wght+Amat(iq,D1_RA)
154            Cmat(iq,D1_GAA) = Cmat(iq,D1_GAA)-9.305257363491002d-1*t1*(t
155     1         25*t3*(1.4721155704335503d-2*t10*t9+1.1040866778251626d-2
156     2         *t7*t8+7.360577852167751d-3*t2*t6+3.6802889260838756d-3*t
157     3         13+1.8401444630419378d-2*t11*t12)+(3.6802889260838756d-3*
158     4         t15*t23+1.8401444630419378d-2*t15*t21*t22+1.4721155704335
159     5         503d-2*t15*t19*t20+1.1040866778251626d-2*t15*t17*t18+7.36
160     6         0577852167751d-3*t14*t15*t16)*t3)*wght
161            Cmat(iq,D1_GAB) = Cmat(iq,D1_GAB)
162          endif ! rhoa.gt.tol_rho
163        else  ! ipol.eq.1
164          rhoa    = rho(iq,R_A)
165          rhob    = rho(iq,R_B)
166          gammaaa = rgamma(iq,G_AA)
167          gammaab = rgamma(iq,G_AB)
168          gammabb = rgamma(iq,G_BB)
169          if (rhoa.gt.tol_rho.and.rhob.gt.tol_rho) then
170            t1 = rhoa**1.3333333333333333d+0
171            t2 = param(1)
172            t3 = param(7)
173            t4 = param(3)
174            t5 = 1/rhoa**2.6666666666666666d+0
175            t6 = 3.6802889260838756d-3*gammaaa*t5
176            t7 = t6+1.0d+0
177            t8 = 1.0d+0-1.0d+0/t7
178            t9 = t8**2.0d+0
179            t10 = param(4)
180            t11 = t8**3.0d+0
181            t12 = param(5)
182            t13 = t8**4.0d+0
183            t14 = param(6)
184            t15 = param(2)
185            t16 = param(9)
186            t17 = exp(-t6)
187            t18 = 1.0d+0-t17
188            t19 = t18**2.0d+0
189            t20 = param(10)
190            t21 = t18**3.0d+0
191            t22 = param(11)
192            t23 = t18**4.0d+0
193            t24 = param(12)
194            t25 = param(8)
195            t26 = t4*t9+t14*t8**5.0d+0+t15*t8+t3+t18*t25+t18**5.0d+0*t24
196     1         +t22*t23+t20*t21+t2+t16*t19+t12*t13+t10*t11
197            t27 = rhob**1.3333333333333333d+0
198            t28 = 1/rhob**2.6666666666666666d+0
199            t29 = 3.6802889260838756d-3*gammabb*t28
200            t30 = t29+1.0d+0
201            t31 = 1.0d+0-1.0d+0/t30
202            t32 = t31**2.0d+0
203            t33 = t31**3.0d+0
204            t34 = t31**4.0d+0
205            t35 = exp(-t29)
206            t36 = 1.0d+0-t35
207            t37 = t36**2.0d+0
208            t38 = t36**3.0d+0
209            t39 = t36**4.0d+0
210            t40 = t32*t4+t22*t39+t20*t38+t16*t37+t24*t36**5.0d+0+t25*t36
211     1         +t12*t34+t10*t33+t14*t31**5.0d+0+t15*t31+t3+t2
212            t41 = 1/t7**2
213            t42 = 1/rhoa**3.6666666666666664d+0
214            t43 = 1/t30**2
215            t44 = 1/rhob**3.6666666666666664d+0
216            t45 = 3.6802889260838756d-3*t15
217            fnc(iq) = (-9.305257363491002d-1*t27*t40-9.305257363491002d-
218     1         1*t1*t26)*wght+fnc(iq)
219            Amat(iq,D1_RA) = (-9.305257363491002d-1*t1*(-2.9442311408671
220     1         007d-2*gammaaa*t10*t41*t42*t9-1.962820760578067d-2*gammaa
221     2         a*t4*t41*t42*t8-9.814103802890335d-3*gammaaa*t15*t41*t42-
222     3         4.9070519014451675d-2*gammaaa*t13*t14*t41*t42-3.925641521
223     4         156134d-2*gammaaa*t11*t12*t41*t42-9.814103802890335d-3*ga
224     5         mmaaa*t17*t25*t42-4.9070519014451675d-2*gammaaa*t17*t23*t
225     6         24*t42-3.925641521156134d-2*gammaaa*t17*t21*t22*t42-2.944
226     7         2311408671007d-2*gammaaa*t17*t19*t20*t42-1.96282076057806
227     8         7d-2*gammaaa*t16*t17*t18*t42)-1.2407009817988002d+0*rhoa*
228     9         *3.333333333333333d-1*t26)*wght+Amat(iq,D1_RA)
229            Amat(iq,D1_RB) = (-9.305257363491002d-1*t27*(-1.962820760578
230     1         067d-2*gammabb*t31*t4*t43*t44-4.9070519014451675d-2*gamma
231     2         bb*t14*t34*t43*t44-3.925641521156134d-2*gammabb*t12*t33*t
232     3         43*t44-2.9442311408671007d-2*gammabb*t10*t32*t43*t44-9.81
233     4         4103802890335d-3*gammabb*t15*t43*t44-4.9070519014451675d-
234     5         2*gammabb*t24*t35*t39*t44-3.925641521156134d-2*gammabb*t2
235     6         2*t35*t38*t44-2.9442311408671007d-2*gammabb*t20*t35*t37*t
236     7         44-1.962820760578067d-2*gammabb*t16*t35*t36*t44-9.8141038
237     8         02890335d-3*gammabb*t25*t35*t44)-1.2407009817988002d+0*rh
238     9         ob**3.333333333333333d-1*t40)*wght+Amat(iq,D1_RB)
239            Cmat(iq,D1_GAA) = Cmat(iq,D1_GAA)-9.305257363491002d-1*t1*(t
240     1         41*t5*(1.1040866778251626d-2*t10*t9+7.360577852167751d-3*
241     2         t4*t8+t45+1.8401444630419378d-2*t13*t14+1.472115570433550
242     3         3d-2*t11*t12)+(3.6802889260838756d-3*t17*t25+1.8401444630
243     4         419378d-2*t17*t23*t24+1.4721155704335503d-2*t17*t21*t22+1
244     5         .1040866778251626d-2*t17*t19*t20+7.360577852167751d-3*t16
245     6         *t17*t18)*t5)*wght
246            Cmat(iq,D1_GAB) = Cmat(iq,D1_GAB)
247            Cmat(iq,D1_GBB) = Cmat(iq,D1_GBB)-9.305257363491002d-1*t27*(
248     1         t28*t43*(t45+7.360577852167751d-3*t31*t4+1.84014446304193
249     2         78d-2*t14*t34+1.4721155704335503d-2*t12*t33+1.10408667782
250     3         51626d-2*t10*t32)+t28*(1.8401444630419378d-2*t24*t35*t39+
251     4         1.4721155704335503d-2*t22*t35*t38+1.1040866778251626d-2*t
252     5         20*t35*t37+7.360577852167751d-3*t16*t35*t36+3.68028892608
253     6         38756d-3*t25*t35))*wght
254          elseif (rhoa.gt.tol_rho.and.rhob.le.tol_rho) then
255            t1 = rhoa**1.3333333333333333d+0
256            t2 = param(3)
257            t3 = 1/rhoa**2.6666666666666666d+0
258            t4 = 3.6802889260838756d-3*gammaaa*t3
259            t5 = t4+1.0d+0
260            t6 = 1.0d+0-1.0d+0/t5
261            t7 = t6**2.0d+0
262            t8 = param(4)
263            t9 = t6**3.0d+0
264            t10 = param(5)
265            t11 = t6**4.0d+0
266            t12 = param(6)
267            t13 = param(2)
268            t14 = param(9)
269            t15 = exp(-t4)
270            t16 = 1.0d+0-t15
271            t17 = t16**2.0d+0
272            t18 = param(10)
273            t19 = t16**3.0d+0
274            t20 = param(11)
275            t21 = t16**4.0d+0
276            t22 = param(12)
277            t23 = param(8)
278            t24 = t8*t9+t2*t7+t12*t6**5.0d+0+t13*t6+t16*t23+t16**5.0d+0*
279     1         t22+t20*t21+t18*t19+t14*t17+t10*t11+param(7)+param(1)
280            t25 = 1/t5**2
281            t26 = 1/rhoa**3.6666666666666664d+0
282            fnc(iq) = fnc(iq)-9.305257363491002d-1*t1*t24*wght
283            Amat(iq,D1_RA) = -9.305257363491002d-1*t1*(-3.92564152115613
284     1         4d-2*gammaaa*t10*t25*t26*t9-2.9442311408671007d-2*gammaaa
285     2         *t25*t26*t7*t8-1.962820760578067d-2*gammaaa*t2*t25*t26*t6
286     3         -9.814103802890335d-3*gammaaa*t13*t25*t26-4.9070519014451
287     4         675d-2*gammaaa*t11*t12*t25*t26-9.814103802890335d-3*gamma
288     5         aa*t15*t23*t26-4.9070519014451675d-2*gammaaa*t15*t21*t22*
289     6         t26-3.925641521156134d-2*gammaaa*t15*t19*t20*t26-2.944231
290     7         1408671007d-2*gammaaa*t15*t17*t18*t26-1.962820760578067d-
291     8         2*gammaaa*t14*t15*t16*t26)*wght-1.2407009817988002d+0*rho
292     9         a**3.333333333333333d-1*t24*wght+Amat(iq,D1_RA)
293            Cmat(iq,D1_GAA) = Cmat(iq,D1_GAA)-9.305257363491002d-1*t1*(t
294     1         25*t3*(1.4721155704335503d-2*t10*t9+1.1040866778251626d-2
295     2         *t7*t8+7.360577852167751d-3*t2*t6+3.6802889260838756d-3*t
296     3         13+1.8401444630419378d-2*t11*t12)+(3.6802889260838756d-3*
297     4         t15*t23+1.8401444630419378d-2*t15*t21*t22+1.4721155704335
298     5         503d-2*t15*t19*t20+1.1040866778251626d-2*t15*t17*t18+7.36
299     6         0577852167751d-3*t14*t15*t16)*t3)*wght
300          elseif (rhoa.le.tol_rho.and.rhob.gt.tol_rho) then
301            t1 = rhob**1.3333333333333333d+0
302            t2 = param(3)
303            t3 = 1/rhob**2.6666666666666666d+0
304            t4 = 3.6802889260838756d-3*gammabb*t3
305            t5 = t4+1.0d+0
306            t6 = 1.0d+0-1.0d+0/t5
307            t7 = t6**2.0d+0
308            t8 = param(4)
309            t9 = t6**3.0d+0
310            t10 = param(5)
311            t11 = t6**4.0d+0
312            t12 = param(6)
313            t13 = param(2)
314            t14 = param(9)
315            t15 = exp(-t4)
316            t16 = 1.0d+0-t15
317            t17 = t16**2.0d+0
318            t18 = param(10)
319            t19 = t16**3.0d+0
320            t20 = param(11)
321            t21 = t16**4.0d+0
322            t22 = param(12)
323            t23 = param(8)
324            t24 = t8*t9+t2*t7+t12*t6**5.0d+0+t13*t6+t16*t23+t16**5.0d+0*
325     1         t22+t20*t21+t18*t19+t14*t17+t10*t11+param(7)+param(1)
326            t25 = 1/t5**2
327            t26 = 1/rhob**3.6666666666666664d+0
328            fnc(iq) = fnc(iq)-9.305257363491002d-1*t1*t24*wght
329            Amat(iq,D1_RB) = -9.305257363491002d-1*t1*(-3.92564152115613
330     1         4d-2*gammabb*t10*t25*t26*t9-2.9442311408671007d-2*gammabb
331     2         *t25*t26*t7*t8-1.962820760578067d-2*gammabb*t2*t25*t26*t6
332     3         -9.814103802890335d-3*gammabb*t13*t25*t26-4.9070519014451
333     4         675d-2*gammabb*t11*t12*t25*t26-9.814103802890335d-3*gamma
334     5         bb*t15*t23*t26-4.9070519014451675d-2*gammabb*t15*t21*t22*
335     6         t26-3.925641521156134d-2*gammabb*t15*t19*t20*t26-2.944231
336     7         1408671007d-2*gammabb*t15*t17*t18*t26-1.962820760578067d-
337     8         2*gammabb*t14*t15*t16*t26)*wght-1.2407009817988002d+0*rho
338     9         b**3.333333333333333d-1*t24*wght+Amat(iq,D1_RB)
339            Cmat(iq,D1_GBB) = Cmat(iq,D1_GBB)-9.305257363491002d-1*t1*(t
340     1         25*t3*(1.4721155704335503d-2*t10*t9+1.1040866778251626d-2
341     2         *t7*t8+7.360577852167751d-3*t2*t6+3.6802889260838756d-3*t
342     3         13+1.8401444630419378d-2*t11*t12)+(3.6802889260838756d-3*
343     4         t15*t23+1.8401444630419378d-2*t15*t21*t22+1.4721155704335
344     5         503d-2*t15*t19*t20+1.1040866778251626d-2*t15*t17*t18+7.36
345     6         0577852167751d-3*t14*t15*t16)*t3)*wght
346          endif ! rhoa.gt.tol_rho.and.rhob.gt.tol_rho
347        endif ! ipol.eq.1
348      enddo ! iq
349      end
350C>
351C> \brief Evaluate the nwxcm_x_sogga functional [1]
352C>
353C> \f{eqnarray*}{
354C>   {\it t_1} &=& {\it param}\left(1\right)\\\\
355C>   {\it t_2} &=& {\it param}\left(7\right)\\\\
356C>   {\it t_3} &=& {\it param}\left(2\right)\\\\
357C>   {\it t_4} &=& {{0.003680288926083876\,
358C>    \sigma_{\alpha\alpha}}\over{\rho_\alpha^{{{8}\over{3}}}}}\\\\
359C>   {\it t_5} &=& 1.0-{{1.0}\over{{\it t_4}+1.0}}\\\\
360C>   {\it t_6} &=& {\it param}\left(3\right)\\\\
361C>   {\it t_7} &=& {\it param}\left(4\right)\\\\
362C>   {\it t_8} &=& {\it param}\left(5\right)\\\\
363C>   {\it t_9} &=& {\it param}\left(6\right)\\\\
364C>   {\it t_{10}} &=& {\it param}\left(9\right)\\\\
365C>   {\it t_{11}} &=& 1.0-{{1}\over{e^{{\it t_4}}}}\\\\
366C>   {\it t_{12}} &=& {\it param}\left(10\right)\\\\
367C>   {\it t_{13}} &=& {\it param}\left(11\right)\\\\
368C>   {\it t_{14}} &=& {\it param}\left(12\right)\\\\
369C>   {\it t_{15}} &=& {\it param}\left(8\right)\\\\
370C>   {\it t_{16}} &=& {{0.003680288926083876\,
371C>    \sigma_{\beta\beta}}\over{\rho_\beta^{{{8}\over{3}}}}}\\\\
372C>   {\it t_{17}} &=& 1.0-{{1.0}\over{{\it t_{16}}+1.0}}\\\\
373C>   {\it t_{18}} &=& 1.0-{{1}\over{e^{{\it t_{16}}}}}\\\\
374C>   {\it t_{19}} &=& {{0.003680288926083876\,\sigma_{ss}}
375C>    \over{\rho_s^{{{8}\over{3}}}}}\\\\
376C>   {\it t_{20}} &=& 1.0-{{1.0}\over{{\it t_{19}}+1.0}}\\\\
377C>   {\it t_{21}} &=& 1.0-{{1}\over{e^{{\it t_{19}}}}}\\\\
378C>   f &=& -0.9305257363491002\,\rho_\beta^{{{4}\over{3}}}\,
379C>    \left({\it t_{15}}\,{\it t_{18}}+{\it t_{14}}\,{
380C>    \it t_{18}}^{5.0}+{\it t_{13}}\,{\it t_{18}}^{4.0}+{
381C>    \it t_{12}}\,{\it t_{18}}^{3.0}+{\it t_{10}}\,{\it t_{18}}^{2.0}
382C>    +{\it t_9}\,{\it t_{17}}^{5.0}+{\it t_8}\,{\it t_{17}}^{4.0}
383C>    +{\it t_7}\,{\it t_{17}}^{3.0}+{\it t_6}\,{\it t_{17}}^{2.0}
384C>    +{\it t_3}\,{\it t_{17}}+{\it t_2}+{\it t_1}\right)
385C>    -0.9305257363491002\,\rho_\alpha^{{{4}\over{3}}}\,\left({
386C>    \it t_{15}}\,{\it t_{11}}+{\it t_{14}}\,{\it t_{11}}^{5.0}+{
387C>    \it t_{13}}\,{\it t_{11}}^{4.0}+{\it t_{12}}\,{\it t_{11}}^{3.0}
388C>    +{\it t_{10}}\,{\it t_{11}}^{2.0}+{\it t_9}\,{\it t_5}^{5.0}
389C>    +{\it t_8}\,{\it t_5}^{4.0}+{\it t_7}\,{\it t_5}^{3.0}+{
390C>    \it t_6}\,{\it t_5}^{2.0}+{\it t_3}\,{\it t_5}+{\it t_2}+{
391C>    \it t_1}\right)\\\\
392C>   g &=& 0\\\\
393C>   G &=& -0.9305257363491002\,\rho_s^{{{4}\over{3}}}\,\left({
394C>    \it t_{15}}\,{\it t_{21}}+{\it t_{14}}\,{\it t_{21}}^{5.0}+{
395C>    \it t_{13}}\,{\it t_{21}}^{4.0}+{\it t_{12}}\,{\it t_{21}}^{3.0}
396C>    +{\it t_{10}}\,{\it t_{21}}^{2.0}+{\it t_9}\,{\it t_{20}}^{5.0}
397C>    +{\it t_8}\,{\it t_{20}}^{4.0}+{\it t_7}\,{\it t_{20}}^{3.0}
398C>    +{\it t_6}\,{\it t_{20}}^{2.0}+{\it t_3}\,{\it t_{20}}+{\it t_2}
399C>    +{\it t_1}\right)\\\\
400C> \f}
401C>
402C> Code generated with Maxima 5.34.0 [2]
403C> driven by autoxc [3].
404C>
405C> ### References ###
406C>
407C> [1] Y Zhao, DG Truhlar, J.Chem.Phys. 128, 184109 (2008)  , DOI:
408C> <a href="https://doi.org/10.1063/1.2912068 ">
409C> 10.1063/1.2912068 </a>
410C>
411C> [2] Maxima, a computer algebra system,
412C> <a href="http://maxima.sourceforge.net/">
413C> http://maxima.sourceforge.net/</a>
414C>
415C> [3] autoxc, revision 27097 2015-05-08
416C>
417      subroutine nwxcm_x_sogga_d2(param,tol_rho,ipol,nq,wght,
418     +rho,rgamma,fnc,Amat,Amat2,Cmat,Cmat2)
419c $Id: $
420#ifdef NWXC_QUAD_PREC
421      implicit real(kind=selected_real_kind(30))(a-h,o-z),integer(i-n)
422      integer, parameter :: rk=selected_real_kind(30)
423#else
424      implicit real(kind=selected_real_kind(15))(a-h,o-z),integer(i-n)
425      integer, parameter :: rk=selected_real_kind(15)
426#endif
427#include "nwxc_param.fh"
428      double precision param(*)     !< [Input] Parameters of functional
429      double precision tol_rho      !< [Input] The lower limit on the density
430      integer ipol                  !< [Input] The number of spin channels
431      integer nq                    !< [Input] The number of points
432      double precision wght         !< [Input] The weight of the functional
433      double precision rho(nq,*)    !< [Input] The density
434      double precision rgamma(nq,*) !< [Input] The norm of the density
435                                    !< gradients
436      double precision fnc(nq)      !< [Output] The value of the functional
437c
438c     Sampling Matrices for the XC Kernel
439c
440      double precision Amat(nq,*)   !< [Output] The derivative wrt rho
441      double precision Cmat(nq,*)   !< [Output] The derivative wrt rgamma
442c
443c     Sampling Matrices for the XC Kernel
444c
445      double precision Amat2(nq,*)  !< [Output] The 2nd derivative wrt rho
446      double precision Cmat2(nq,*)  !< [Output] The 2nd derivative wrt rgamma
447                                    !< and possibly rho
448      integer iq
449      double precision tmp
450      double precision rhoa,rhob
451      double precision gammaaa,gammaab,gammabb
452      double precision taua,taub
453      double precision nwxcm_heaviside
454      external         nwxcm_heaviside
455CDIR$ NOVECTOR
456      do iq = 1, nq
457        if (ipol.eq.1) then
458          rhoa    = 0.5d0*rho(iq,R_T)
459          gammaaa = 0.25d0*rgamma(iq,G_TT)
460          if (rhoa.gt.tol_rho) then
461            t1 = rhoa**1.3333333333333333d+0
462            t2 = param(3)
463            t3 = 1/rhoa**2.6666666666666666d+0
464            t4 = 3.6802889260838756d-3*gammaaa*t3
465            t5 = t4+1.0d+0
466            t6 = 1.0d+0-1.0d+0/t5
467            t7 = t6**2.0d+0
468            t8 = param(4)
469            t9 = t6**3.0d+0
470            t10 = param(5)
471            t11 = t6**4.0d+0
472            t12 = param(6)
473            t13 = param(2)
474            t14 = param(9)
475            t15 = exp(-t4)
476            t16 = 1.0d+0-t15
477            t17 = t16**2.0d+0
478            t18 = param(10)
479            t19 = t16**3.0d+0
480            t20 = param(11)
481            t21 = t16**4.0d+0
482            t22 = param(12)
483            t23 = param(8)
484            t24 = t8*t9+t2*t7+t12*t6**5.0d+0+t13*t6+t16*t23+t16**5.0d+0*
485     1         t22+t20*t21+t18*t19+t14*t17+t10*t11+param(7)+param(1)
486            t25 = rhoa**3.333333333333333d-1
487            t26 = 1/t5**2
488            t27 = 1/rhoa**3.6666666666666664d+0
489            t28 = -3.925641521156134d-2*gammaaa*t10*t26*t27*t9-2.9442311
490     1         408671007d-2*gammaaa*t26*t27*t7*t8-1.962820760578067d-2*g
491     2         ammaaa*t2*t26*t27*t6-9.814103802890335d-3*gammaaa*t13*t26
492     3         *t27-4.9070519014451675d-2*gammaaa*t11*t12*t26*t27-9.8141
493     4         03802890335d-3*gammaaa*t15*t23*t27-4.9070519014451675d-2*
494     5         gammaaa*t15*t21*t22*t27-3.925641521156134d-2*gammaaa*t15*
495     6         t19*t20*t27-2.9442311408671007d-2*gammaaa*t15*t17*t18*t27
496     7         -1.962820760578067d-2*gammaaa*t14*t15*t16*t27
497            t29 = 1.4721155704335503d-2*t10*t9+1.1040866778251626d-2*t7*
498     1         t8+7.360577852167751d-3*t2*t6+3.6802889260838756d-3*t13+1
499     2         .8401444630419378d-2*t11*t12
500            t30 = gammaaa**2
501            t31 = 1/t5**4
502            t32 = 1/rhoa**7.333333333333333d+0
503            t33 = 1/t5**3
504            t34 = 1/rhoa**4.666666666666667d+0
505            t35 = exp(-7.360577852167751d-3*gammaaa*t3)
506            t36 = 1/rhoa**6.333333333333333d+0
507            fnc(iq) = fnc(iq)-1.8610514726982003d+0*t1*t24*wght
508            Amat(iq,D1_RA) = (-9.305257363491002d-1*t1*t28-1.24070098179
509     1         88002d+0*t24*t25)*wght+Amat(iq,D1_RA)
510            Cmat(iq,D1_GAA) = Cmat(iq,D1_GAA)-9.305257363491002d-1*t1*(t
511     1         26*t29*t3+(3.6802889260838756d-3*t15*t23+1.84014446304193
512     2         78d-2*t15*t21*t22+1.4721155704335503d-2*t15*t19*t20+1.104
513     3         0866778251626d-2*t15*t17*t18+7.360577852167751d-3*t14*t15
514     4         *t16)*t3)*wght
515            Cmat(iq,D1_GAB) = Cmat(iq,D1_GAB)
516            Amat2(iq,D2_RA_RA) = (-9.305257363491002d-1*t1*(1.4394018910
517     1         905823d-1*gammaaa*t10*t26*t34*t9-7.705330676312523d-4*t10
518     2         *t30*t32*t33*t9+1.9263326690781307d-3*t12*t30*t31*t32*t9+
519     3         1.0795514183179368d-1*gammaaa*t26*t34*t7*t8-5.77899800723
520     4         4393d-4*t30*t32*t33*t7*t8+5.778998007234393d-4*t30*t31*t3
521     5         2*t6*t8+1.1557996014468785d-3*t10*t30*t31*t32*t7+7.197009
522     6         455452912d-2*gammaaa*t2*t26*t34*t6-3.8526653381562614d-4*
523     7         t2*t30*t32*t33*t6+1.9263326690781307d-3*t19*t22*t30*t32*t
524     8         35+1.1557996014468785d-3*t17*t20*t30*t32*t35+5.7789980072
525     9         34393d-4*t16*t18*t30*t32*t35+1.9263326690781307d-4*t14*t3
526     :         0*t32*t35+3.598504727726456d-2*gammaaa*t13*t26*t34+1.7992
527     ;         52363863228d-1*gammaaa*t11*t12*t26*t34+3.598504727726456d
528     <         -2*gammaaa*t15*t23*t34+1.799252363863228d-1*gammaaa*t15*t
529     =         21*t22*t34+1.4394018910905823d-1*gammaaa*t15*t19*t20*t34+
530     >         1.0795514183179368d-1*gammaaa*t15*t17*t18*t34+7.197009455
531     ?         452912d-2*gammaaa*t14*t15*t16*t34-1.9263326690781307d-4*t
532     @         13*t30*t32*t33-9.631663345390654d-4*t11*t12*t30*t32*t33+1
533     1         .9263326690781307d-4*t2*t30*t31*t32-9.631663345390653d-5*
534     2         t15*t23*t30*t32-4.815831672695327d-4*t15*t21*t22*t30*t32-
535     3         3.8526653381562614d-4*t15*t19*t20*t30*t32-2.8894990036171
536     4         964d-4*t15*t17*t18*t30*t32-1.9263326690781307d-4*t14*t15*
537     5         t16*t30*t32)-2.4814019635976003d+0*t25*t28-4.135669939329
538     6         334d-1*t24/rhoa**6.666666666666666d-1)*wght+Amat2(iq,D2_R
539     7         A_RA)
540            Amat2(iq,D2_RA_RB) = Amat2(iq,D2_RA_RB)
541            Cmat2(iq,D2_RA_GAA) = (-9.305257363491002d-1*t1*(2.889499003
542     1         6171964d-4*gammaaa*t10*t33*t36*t9-7.22374750904299d-4*gam
543     2         maaa*t12*t31*t36*t9-3.925641521156134d-2*t10*t26*t27*t9+2
544     3         .167124252712897d-4*gammaaa*t33*t36*t7*t8-2.9442311408671
545     4         007d-2*t26*t27*t7*t8-2.167124252712897d-4*gammaaa*t31*t36
546     5         *t6*t8-4.334248505425794d-4*gammaaa*t10*t31*t36*t7+1.4447
547     6         495018085982d-4*gammaaa*t2*t33*t36*t6-1.962820760578067d-
548     7         2*t2*t26*t27*t6-7.22374750904299d-4*gammaaa*t19*t22*t35*t
549     8         36-4.334248505425794d-4*gammaaa*t17*t20*t35*t36-2.1671242
550     9         52712897d-4*gammaaa*t16*t18*t35*t36-7.22374750904299d-5*g
551     :         ammaaa*t14*t35*t36+7.22374750904299d-5*gammaaa*t13*t33*t3
552     ;         6+3.611873754521495d-4*gammaaa*t11*t12*t33*t36-7.22374750
553     <         904299d-5*gammaaa*t2*t31*t36+3.611873754521495d-5*gammaaa
554     =         *t15*t23*t36+1.8059368772607476d-4*gammaaa*t15*t21*t22*t3
555     >         6+1.4447495018085982d-4*gammaaa*t15*t19*t20*t36+1.0835621
556     ?         263564485d-4*gammaaa*t15*t17*t18*t36+7.22374750904299d-5*
557     @         gammaaa*t14*t15*t16*t36-9.814103802890335d-3*t13*t26*t27-
558     1         4.9070519014451675d-2*t11*t12*t26*t27-9.814103802890335d-
559     2         3*t15*t23*t27-4.9070519014451675d-2*t15*t21*t22*t27-3.925
560     3         641521156134d-2*t15*t19*t20*t27-2.9442311408671007d-2*t15
561     4         *t17*t18*t27-1.962820760578067d-2*t14*t15*t16*t27)-1.2407
562     5         009817988002d+0*t25*(1.4721155704335503d-2*t10*t26*t3*t9+
563     6         1.1040866778251626d-2*t26*t3*t7*t8+7.360577852167751d-3*t
564     7         2*t26*t3*t6+3.6802889260838756d-3*t13*t26*t3+1.8401444630
565     8         419378d-2*t11*t12*t26*t3+3.6802889260838756d-3*t15*t23*t3
566     9         +1.8401444630419378d-2*t15*t21*t22*t3+1.4721155704335503d
567     :         -2*t15*t19*t20*t3+1.1040866778251626d-2*t15*t17*t18*t3+7.
568     ;         360577852167751d-3*t14*t15*t16*t3))*wght+Cmat2(iq,D2_RA_G
569     <         AA)
570            Cmat2(iq,D2_RA_GAB) = Cmat2(iq,D2_RA_GAB)
571            Cmat2(iq,D2_RA_GBB) = Cmat2(iq,D2_RA_GBB)
572            Cmat2(iq,D2_GAA_GAA) = Cmat2(iq,D2_GAA_GAA)-9.30525736349100
573     1         2d-1*t1*(t26*t3*(2.708905315891122d-4*t12*t26*t3*t9+8.126
574     2         715947673364d-5*t26*t3*t6*t8+1.625343189534673d-4*t10*t26
575     3         *t3*t7+2.7089053158911214d-5*t2*t26*t3)+t3*(2.70890531589
576     4         1122d-4*t19*t22*t3*t35+1.625343189534673d-4*t17*t20*t3*t3
577     5         5+8.126715947673364d-5*t16*t18*t3*t35+2.7089053158911214d
578     6         -5*t14*t3*t35-1.3544526579455607d-5*t15*t23*t3-6.77226328
579     7         9727805d-5*t15*t21*t22*t3-5.417810631782243d-5*t15*t19*t2
580     8         0*t3-4.063357973836682d-5*t15*t17*t18*t3-2.70890531589112
581     9         14d-5*t14*t15*t16*t3)-7.360577852167751d-3*t29*t33/rhoa**
582     :         5.333333333333333d+0)*wght
583            Cmat2(iq,D2_GAA_GAB) = Cmat2(iq,D2_GAA_GAB)
584            Cmat2(iq,D2_GAA_GBB) = Cmat2(iq,D2_GAA_GBB)
585            Cmat2(iq,D2_GAB_GAB) = Cmat2(iq,D2_GAB_GAB)
586          endif ! rhoa.gt.tol_rho
587        else  ! ipol.eq.1
588          rhoa    = rho(iq,R_A)
589          rhob    = rho(iq,R_B)
590          gammaaa = rgamma(iq,G_AA)
591          gammaab = rgamma(iq,G_AB)
592          gammabb = rgamma(iq,G_BB)
593          if (rhoa.gt.tol_rho.and.rhob.gt.tol_rho) then
594            t1 = rhoa**1.3333333333333333d+0
595            t2 = param(1)
596            t3 = param(7)
597            t4 = param(3)
598            t5 = 1/rhoa**2.6666666666666666d+0
599            t6 = 3.6802889260838756d-3*gammaaa*t5
600            t7 = t6+1.0d+0
601            t8 = 1.0d+0-1.0d+0/t7
602            t9 = t8**2.0d+0
603            t10 = param(4)
604            t11 = t8**3.0d+0
605            t12 = param(5)
606            t13 = t8**4.0d+0
607            t14 = param(6)
608            t15 = param(2)
609            t16 = param(9)
610            t17 = exp(-t6)
611            t18 = 1.0d+0-t17
612            t19 = t18**2.0d+0
613            t20 = param(10)
614            t21 = t18**3.0d+0
615            t22 = param(11)
616            t23 = t18**4.0d+0
617            t24 = param(12)
618            t25 = param(8)
619            t26 = t4*t9+t14*t8**5.0d+0+t15*t8+t3+t18*t25+t18**5.0d+0*t24
620     1         +t22*t23+t20*t21+t2+t16*t19+t12*t13+t10*t11
621            t27 = rhob**1.3333333333333333d+0
622            t28 = 1/rhob**2.6666666666666666d+0
623            t29 = 3.6802889260838756d-3*gammabb*t28
624            t30 = t29+1.0d+0
625            t31 = 1.0d+0-1.0d+0/t30
626            t32 = t31**2.0d+0
627            t33 = t31**3.0d+0
628            t34 = t31**4.0d+0
629            t35 = exp(-t29)
630            t36 = 1.0d+0-t35
631            t37 = t36**2.0d+0
632            t38 = t36**3.0d+0
633            t39 = t36**4.0d+0
634            t40 = t32*t4+t22*t39+t20*t38+t16*t37+t24*t36**5.0d+0+t25*t36
635     1         +t12*t34+t10*t33+t14*t31**5.0d+0+t15*t31+t3+t2
636            t41 = rhoa**3.333333333333333d-1
637            t42 = 1/t7**2
638            t43 = 1/rhoa**3.6666666666666664d+0
639            t44 = -2.9442311408671007d-2*gammaaa*t10*t42*t43*t9-1.962820
640     1         760578067d-2*gammaaa*t4*t42*t43*t8-9.814103802890335d-3*g
641     2         ammaaa*t15*t42*t43-4.9070519014451675d-2*gammaaa*t13*t14*
642     3         t42*t43-3.925641521156134d-2*gammaaa*t11*t12*t42*t43-9.81
643     4         4103802890335d-3*gammaaa*t17*t25*t43-4.9070519014451675d-
644     5         2*gammaaa*t17*t23*t24*t43-3.925641521156134d-2*gammaaa*t1
645     6         7*t21*t22*t43-2.9442311408671007d-2*gammaaa*t17*t19*t20*t
646     7         43-1.962820760578067d-2*gammaaa*t16*t17*t18*t43
647            t45 = rhob**3.333333333333333d-1
648            t46 = 1/t30**2
649            t47 = 1/rhob**3.6666666666666664d+0
650            t48 = -1.962820760578067d-2*gammabb*t31*t4*t46*t47-4.9070519
651     1         014451675d-2*gammabb*t14*t34*t46*t47-3.925641521156134d-2
652     2         *gammabb*t12*t33*t46*t47-2.9442311408671007d-2*gammabb*t1
653     3         0*t32*t46*t47-9.814103802890335d-3*gammabb*t15*t46*t47-4.
654     4         9070519014451675d-2*gammabb*t24*t35*t39*t47-3.92564152115
655     5         6134d-2*gammabb*t22*t35*t38*t47-2.9442311408671007d-2*gam
656     6         mabb*t20*t35*t37*t47-1.962820760578067d-2*gammabb*t16*t35
657     7         *t36*t47-9.814103802890335d-3*gammabb*t25*t35*t47
658            t49 = 3.6802889260838756d-3*t15
659            t50 = 1.1040866778251626d-2*t10*t9+7.360577852167751d-3*t4*t
660     1         8+t49+1.8401444630419378d-2*t13*t14+1.4721155704335503d-2
661     2         *t11*t12
662            t51 = t49+7.360577852167751d-3*t31*t4+1.8401444630419378d-2*
663     1         t14*t34+1.4721155704335503d-2*t12*t33+1.1040866778251626d
664     2         -2*t10*t32
665            t52 = gammaaa**2
666            t53 = 1/t7**4
667            t54 = 1/rhoa**7.333333333333333d+0
668            t55 = 1/t7**3
669            t56 = 1/rhoa**4.666666666666667d+0
670            t57 = exp(-7.360577852167751d-3*gammaaa*t5)
671            t58 = gammabb**2
672            t59 = 1/t30**4
673            t60 = 1/rhob**7.333333333333333d+0
674            t61 = 1/t30**3
675            t62 = 1/rhob**4.666666666666667d+0
676            t63 = exp(-7.360577852167751d-3*gammabb*t28)
677            t64 = 1/rhoa**6.333333333333333d+0
678            t65 = 1/rhob**6.333333333333333d+0
679            fnc(iq) = (-9.305257363491002d-1*t27*t40-9.305257363491002d-
680     1         1*t1*t26)*wght+fnc(iq)
681            Amat(iq,D1_RA) = (-9.305257363491002d-1*t1*t44-1.24070098179
682     1         88002d+0*t26*t41)*wght+Amat(iq,D1_RA)
683            Amat(iq,D1_RB) = (-9.305257363491002d-1*t27*t48-1.2407009817
684     1         988002d+0*t40*t45)*wght+Amat(iq,D1_RB)
685            Cmat(iq,D1_GAA) = Cmat(iq,D1_GAA)-9.305257363491002d-1*t1*(t
686     1         42*t5*t50+(3.6802889260838756d-3*t17*t25+1.84014446304193
687     2         78d-2*t17*t23*t24+1.4721155704335503d-2*t17*t21*t22+1.104
688     3         0866778251626d-2*t17*t19*t20+7.360577852167751d-3*t16*t17
689     4         *t18)*t5)*wght
690            Cmat(iq,D1_GAB) = Cmat(iq,D1_GAB)
691            Cmat(iq,D1_GBB) = Cmat(iq,D1_GBB)-9.305257363491002d-1*t27*(
692     1         t28*t46*t51+t28*(1.8401444630419378d-2*t24*t35*t39+1.4721
693     2         155704335503d-2*t22*t35*t38+1.1040866778251626d-2*t20*t35
694     3         *t37+7.360577852167751d-3*t16*t35*t36+3.6802889260838756d
695     4         -3*t25*t35))*wght
696            Amat2(iq,D2_RA_RA) = (-9.305257363491002d-1*t1*(1.0795514183
697     1         179368d-1*gammaaa*t10*t42*t56*t9-5.778998007234393d-4*t10
698     2         *t52*t54*t55*t9+1.1557996014468785d-3*t12*t52*t53*t54*t9+
699     3         7.197009455452912d-2*gammaaa*t4*t42*t56*t8-3.852665338156
700     4         2614d-4*t4*t52*t54*t55*t8+5.778998007234393d-4*t10*t52*t5
701     5         3*t54*t8+1.9263326690781307d-3*t21*t24*t52*t54*t57+1.1557
702     6         996014468785d-3*t19*t22*t52*t54*t57+5.778998007234393d-4*
703     7         t18*t20*t52*t54*t57+1.9263326690781307d-4*t16*t52*t54*t57
704     8         +3.598504727726456d-2*gammaaa*t15*t42*t56+1.7992523638632
705     9         28d-1*gammaaa*t13*t14*t42*t56+1.4394018910905823d-1*gamma
706     :         aa*t11*t12*t42*t56+3.598504727726456d-2*gammaaa*t17*t25*t
707     ;         56+1.799252363863228d-1*gammaaa*t17*t23*t24*t56+1.4394018
708     <         910905823d-1*gammaaa*t17*t21*t22*t56+1.0795514183179368d-
709     =         1*gammaaa*t17*t19*t20*t56+7.197009455452912d-2*gammaaa*t1
710     >         6*t17*t18*t56-1.9263326690781307d-4*t15*t52*t54*t55-9.631
711     ?         663345390654d-4*t13*t14*t52*t54*t55-7.705330676312523d-4*
712     @         t11*t12*t52*t54*t55+1.9263326690781307d-4*t4*t52*t53*t54+
713     1         1.9263326690781307d-3*t11*t14*t52*t53*t54-9.6316633453906
714     2         53d-5*t17*t25*t52*t54-4.815831672695327d-4*t17*t23*t24*t5
715     3         2*t54-3.8526653381562614d-4*t17*t21*t22*t52*t54-2.8894990
716     4         036171964d-4*t17*t19*t20*t52*t54-1.9263326690781307d-4*t1
717     5         6*t17*t18*t52*t54)-2.4814019635976003d+0*t41*t44-4.135669
718     6         939329334d-1*t26/rhoa**6.666666666666666d-1)*wght+Amat2(i
719     7         q,D2_RA_RA)
720            Amat2(iq,D2_RA_RB) = Amat2(iq,D2_RA_RB)
721            Amat2(iq,D2_RB_RB) = (-9.305257363491002d-1*t27*(1.926332669
722     1         0781307d-3*t24*t38*t58*t60*t63+1.1557996014468785d-3*t22*
723     2         t37*t58*t60*t63+5.778998007234393d-4*t20*t36*t58*t60*t63+
724     3         1.9263326690781307d-4*t16*t58*t60*t63+7.197009455452912d-
725     4         2*gammabb*t31*t4*t46*t62+1.799252363863228d-1*gammabb*t14
726     5         *t34*t46*t62+1.4394018910905823d-1*gammabb*t12*t33*t46*t6
727     6         2+1.0795514183179368d-1*gammabb*t10*t32*t46*t62+3.5985047
728     7         27726456d-2*gammabb*t15*t46*t62+1.799252363863228d-1*gamm
729     8         abb*t24*t35*t39*t62+1.4394018910905823d-1*gammabb*t22*t35
730     9         *t38*t62+1.0795514183179368d-1*gammabb*t20*t35*t37*t62+7.
731     :         197009455452912d-2*gammabb*t16*t35*t36*t62+3.598504727726
732     ;         456d-2*gammabb*t25*t35*t62-3.8526653381562614d-4*t31*t4*t
733     <         58*t60*t61-9.631663345390654d-4*t14*t34*t58*t60*t61-7.705
734     =         330676312523d-4*t12*t33*t58*t60*t61-5.778998007234393d-4*
735     >         t10*t32*t58*t60*t61-1.9263326690781307d-4*t15*t58*t60*t61
736     ?         +1.9263326690781307d-4*t4*t58*t59*t60+1.9263326690781307d
737     @         -3*t14*t33*t58*t59*t60+1.1557996014468785d-3*t12*t32*t58*
738     1         t59*t60+5.778998007234393d-4*t10*t31*t58*t59*t60-4.815831
739     2         672695327d-4*t24*t35*t39*t58*t60-3.8526653381562614d-4*t2
740     3         2*t35*t38*t58*t60-2.8894990036171964d-4*t20*t35*t37*t58*t
741     4         60-1.9263326690781307d-4*t16*t35*t36*t58*t60-9.6316633453
742     5         90653d-5*t25*t35*t58*t60)-2.4814019635976003d+0*t45*t48-4
743     6         .135669939329334d-1*t40/rhob**6.666666666666666d-1)*wght+
744     7         Amat2(iq,D2_RB_RB)
745            Cmat2(iq,D2_RA_GAA) = (-9.305257363491002d-1*t1*(2.167124252
746     1         712897d-4*gammaaa*t10*t55*t64*t9-4.334248505425794d-4*gam
747     2         maaa*t12*t53*t64*t9-2.9442311408671007d-2*t10*t42*t43*t9+
748     3         1.4447495018085982d-4*gammaaa*t4*t55*t64*t8-2.16712425271
749     4         2897d-4*gammaaa*t10*t53*t64*t8-1.962820760578067d-2*t4*t4
750     5         2*t43*t8-7.22374750904299d-4*gammaaa*t21*t24*t57*t64-4.33
751     6         4248505425794d-4*gammaaa*t19*t22*t57*t64-2.16712425271289
752     7         7d-4*gammaaa*t18*t20*t57*t64-7.22374750904299d-5*gammaaa*
753     8         t16*t57*t64+7.22374750904299d-5*gammaaa*t15*t55*t64+3.611
754     9         873754521495d-4*gammaaa*t13*t14*t55*t64+2.889499003617196
755     :         4d-4*gammaaa*t11*t12*t55*t64-7.22374750904299d-5*gammaaa*
756     ;         t4*t53*t64-7.22374750904299d-4*gammaaa*t11*t14*t53*t64+3.
757     <         611873754521495d-5*gammaaa*t17*t25*t64+1.8059368772607476
758     =         d-4*gammaaa*t17*t23*t24*t64+1.4447495018085982d-4*gammaaa
759     >         *t17*t21*t22*t64+1.0835621263564485d-4*gammaaa*t17*t19*t2
760     ?         0*t64+7.22374750904299d-5*gammaaa*t16*t17*t18*t64-9.81410
761     @         3802890335d-3*t15*t42*t43-4.9070519014451675d-2*t13*t14*t
762     1         42*t43-3.925641521156134d-2*t11*t12*t42*t43-9.81410380289
763     2         0335d-3*t17*t25*t43-4.9070519014451675d-2*t17*t23*t24*t43
764     3         -3.925641521156134d-2*t17*t21*t22*t43-2.9442311408671007d
765     4         -2*t17*t19*t20*t43-1.962820760578067d-2*t16*t17*t18*t43)-
766     5         1.2407009817988002d+0*t41*(1.1040866778251626d-2*t10*t42*
767     6         t5*t9+7.360577852167751d-3*t4*t42*t5*t8+3.680288926083875
768     7         6d-3*t15*t42*t5+1.8401444630419378d-2*t13*t14*t42*t5+1.47
769     8         21155704335503d-2*t11*t12*t42*t5+3.6802889260838756d-3*t1
770     9         7*t25*t5+1.8401444630419378d-2*t17*t23*t24*t5+1.472115570
771     :         4335503d-2*t17*t21*t22*t5+1.1040866778251626d-2*t17*t19*t
772     ;         20*t5+7.360577852167751d-3*t16*t17*t18*t5))*wght+Cmat2(iq
773     <         ,D2_RA_GAA)
774            Cmat2(iq,D2_RA_GAB) = Cmat2(iq,D2_RA_GAB)
775            Cmat2(iq,D2_RA_GBB) = Cmat2(iq,D2_RA_GBB)
776            Cmat2(iq,D2_RB_GAA) = Cmat2(iq,D2_RB_GAA)
777            Cmat2(iq,D2_RB_GAB) = Cmat2(iq,D2_RB_GAB)
778            Cmat2(iq,D2_RB_GBB) = (-9.305257363491002d-1*t27*(-7.2237475
779     1         0904299d-4*gammabb*t24*t38*t63*t65-4.334248505425794d-4*g
780     2         ammabb*t22*t37*t63*t65-2.167124252712897d-4*gammabb*t20*t
781     3         36*t63*t65-7.22374750904299d-5*gammabb*t16*t63*t65+1.4447
782     4         495018085982d-4*gammabb*t31*t4*t61*t65+3.611873754521495d
783     5         -4*gammabb*t14*t34*t61*t65+2.8894990036171964d-4*gammabb*
784     6         t12*t33*t61*t65+2.167124252712897d-4*gammabb*t10*t32*t61*
785     7         t65+7.22374750904299d-5*gammabb*t15*t61*t65-7.22374750904
786     8         299d-5*gammabb*t4*t59*t65-7.22374750904299d-4*gammabb*t14
787     9         *t33*t59*t65-4.334248505425794d-4*gammabb*t12*t32*t59*t65
788     :         -2.167124252712897d-4*gammabb*t10*t31*t59*t65+1.805936877
789     ;         2607476d-4*gammabb*t24*t35*t39*t65+1.4447495018085982d-4*
790     <         gammabb*t22*t35*t38*t65+1.0835621263564485d-4*gammabb*t20
791     =         *t35*t37*t65+7.22374750904299d-5*gammabb*t16*t35*t36*t65+
792     >         3.611873754521495d-5*gammabb*t25*t35*t65-1.96282076057806
793     ?         7d-2*t31*t4*t46*t47-4.9070519014451675d-2*t14*t34*t46*t47
794     @         -3.925641521156134d-2*t12*t33*t46*t47-2.9442311408671007d
795     1         -2*t10*t32*t46*t47-9.814103802890335d-3*t15*t46*t47-4.907
796     2         0519014451675d-2*t24*t35*t39*t47-3.925641521156134d-2*t22
797     3         *t35*t38*t47-2.9442311408671007d-2*t20*t35*t37*t47-1.9628
798     4         20760578067d-2*t16*t35*t36*t47-9.814103802890335d-3*t25*t
799     5         35*t47)-1.2407009817988002d+0*t45*(7.360577852167751d-3*t
800     6         28*t31*t4*t46+1.8401444630419378d-2*t14*t28*t34*t46+1.472
801     7         1155704335503d-2*t12*t28*t33*t46+1.1040866778251626d-2*t1
802     8         0*t28*t32*t46+3.6802889260838756d-3*t15*t28*t46+1.8401444
803     9         630419378d-2*t24*t28*t35*t39+1.4721155704335503d-2*t22*t2
804     :         8*t35*t38+1.1040866778251626d-2*t20*t28*t35*t37+7.3605778
805     ;         52167751d-3*t16*t28*t35*t36+3.6802889260838756d-3*t25*t28
806     <         *t35))*wght+Cmat2(iq,D2_RB_GBB)
807            Cmat2(iq,D2_GAA_GAA) = Cmat2(iq,D2_GAA_GAA)-9.30525736349100
808     1         2d-1*t1*(t42*t5*(1.625343189534673d-4*t12*t42*t5*t9+8.126
809     2         715947673364d-5*t10*t42*t5*t8+2.7089053158911214d-5*t4*t4
810     3         2*t5+2.708905315891122d-4*t11*t14*t42*t5)+t5*(2.708905315
811     4         891122d-4*t21*t24*t5*t57+1.625343189534673d-4*t19*t22*t5*
812     5         t57+8.126715947673364d-5*t18*t20*t5*t57+2.708905315891121
813     6         4d-5*t16*t5*t57-1.3544526579455607d-5*t17*t25*t5-6.772263
814     7         289727805d-5*t17*t23*t24*t5-5.417810631782243d-5*t17*t21*
815     8         t22*t5-4.063357973836682d-5*t17*t19*t20*t5-2.708905315891
816     9         1214d-5*t16*t17*t18*t5)-7.360577852167751d-3*t50*t55/rhoa
817     :         **5.333333333333333d+0)*wght
818            Cmat2(iq,D2_GAA_GAB) = Cmat2(iq,D2_GAA_GAB)
819            Cmat2(iq,D2_GAA_GBB) = Cmat2(iq,D2_GAA_GBB)
820            Cmat2(iq,D2_GAB_GAB) = Cmat2(iq,D2_GAB_GAB)
821            Cmat2(iq,D2_GAB_GBB) = Cmat2(iq,D2_GAB_GBB)
822            Cmat2(iq,D2_GBB_GBB) = Cmat2(iq,D2_GBB_GBB)-9.30525736349100
823     1         2d-1*t27*(t28*(2.708905315891122d-4*t24*t28*t38*t63+1.625
824     2         343189534673d-4*t22*t28*t37*t63+8.126715947673364d-5*t20*
825     3         t28*t36*t63+2.7089053158911214d-5*t16*t28*t63-6.772263289
826     4         727805d-5*t24*t28*t35*t39-5.417810631782243d-5*t22*t28*t3
827     5         5*t38-4.063357973836682d-5*t20*t28*t35*t37-2.708905315891
828     6         1214d-5*t16*t28*t35*t36-1.3544526579455607d-5*t25*t28*t35
829     7         )-7.360577852167751d-3*t51*t61/rhob**5.333333333333333d+0
830     8         +t28*t46*(2.7089053158911214d-5*t28*t4*t46+2.708905315891
831     9         122d-4*t14*t28*t33*t46+1.625343189534673d-4*t12*t28*t32*t
832     :         46+8.126715947673364d-5*t10*t28*t31*t46))*wght
833          elseif (rhoa.gt.tol_rho.and.rhob.le.tol_rho) then
834            t1 = rhoa**1.3333333333333333d+0
835            t2 = param(3)
836            t3 = 1/rhoa**2.6666666666666666d+0
837            t4 = 3.6802889260838756d-3*gammaaa*t3
838            t5 = t4+1.0d+0
839            t6 = 1.0d+0-1.0d+0/t5
840            t7 = t6**2.0d+0
841            t8 = param(4)
842            t9 = t6**3.0d+0
843            t10 = param(5)
844            t11 = t6**4.0d+0
845            t12 = param(6)
846            t13 = param(2)
847            t14 = param(9)
848            t15 = exp(-t4)
849            t16 = 1.0d+0-t15
850            t17 = t16**2.0d+0
851            t18 = param(10)
852            t19 = t16**3.0d+0
853            t20 = param(11)
854            t21 = t16**4.0d+0
855            t22 = param(12)
856            t23 = param(8)
857            t24 = t8*t9+t2*t7+t12*t6**5.0d+0+t13*t6+t16*t23+t16**5.0d+0*
858     1         t22+t20*t21+t18*t19+t14*t17+t10*t11+param(7)+param(1)
859            t25 = rhoa**3.333333333333333d-1
860            t26 = 1/t5**2
861            t27 = 1/rhoa**3.6666666666666664d+0
862            t28 = -3.925641521156134d-2*gammaaa*t10*t26*t27*t9-2.9442311
863     1         408671007d-2*gammaaa*t26*t27*t7*t8-1.962820760578067d-2*g
864     2         ammaaa*t2*t26*t27*t6-9.814103802890335d-3*gammaaa*t13*t26
865     3         *t27-4.9070519014451675d-2*gammaaa*t11*t12*t26*t27-9.8141
866     4         03802890335d-3*gammaaa*t15*t23*t27-4.9070519014451675d-2*
867     5         gammaaa*t15*t21*t22*t27-3.925641521156134d-2*gammaaa*t15*
868     6         t19*t20*t27-2.9442311408671007d-2*gammaaa*t15*t17*t18*t27
869     7         -1.962820760578067d-2*gammaaa*t14*t15*t16*t27
870            t29 = 1.4721155704335503d-2*t10*t9+1.1040866778251626d-2*t7*
871     1         t8+7.360577852167751d-3*t2*t6+3.6802889260838756d-3*t13+1
872     2         .8401444630419378d-2*t11*t12
873            t30 = gammaaa**2
874            t31 = 1/t5**4
875            t32 = 1/rhoa**7.333333333333333d+0
876            t33 = 1/t5**3
877            t34 = 1/rhoa**4.666666666666667d+0
878            t35 = exp(-7.360577852167751d-3*gammaaa*t3)
879            t36 = 1/rhoa**6.333333333333333d+0
880            fnc(iq) = fnc(iq)-9.305257363491002d-1*t1*t24*wght
881            Amat(iq,D1_RA) = -9.305257363491002d-1*t1*t28*wght-1.2407009
882     1         817988002d+0*t24*t25*wght+Amat(iq,D1_RA)
883            Cmat(iq,D1_GAA) = Cmat(iq,D1_GAA)-9.305257363491002d-1*t1*(t
884     1         26*t29*t3+(3.6802889260838756d-3*t15*t23+1.84014446304193
885     2         78d-2*t15*t21*t22+1.4721155704335503d-2*t15*t19*t20+1.104
886     3         0866778251626d-2*t15*t17*t18+7.360577852167751d-3*t14*t15
887     4         *t16)*t3)*wght
888            Amat2(iq,D2_RA_RA) = -9.305257363491002d-1*t1*(1.43940189109
889     1         05823d-1*gammaaa*t10*t26*t34*t9-7.705330676312523d-4*t10*
890     2         t30*t32*t33*t9+1.9263326690781307d-3*t12*t30*t31*t32*t9+1
891     3         .0795514183179368d-1*gammaaa*t26*t34*t7*t8-5.778998007234
892     4         393d-4*t30*t32*t33*t7*t8+5.778998007234393d-4*t30*t31*t32
893     5         *t6*t8+1.1557996014468785d-3*t10*t30*t31*t32*t7+7.1970094
894     6         55452912d-2*gammaaa*t2*t26*t34*t6-3.8526653381562614d-4*t
895     7         2*t30*t32*t33*t6+1.9263326690781307d-3*t19*t22*t30*t32*t3
896     8         5+1.1557996014468785d-3*t17*t20*t30*t32*t35+5.77899800723
897     9         4393d-4*t16*t18*t30*t32*t35+1.9263326690781307d-4*t14*t30
898     :         *t32*t35+3.598504727726456d-2*gammaaa*t13*t26*t34+1.79925
899     ;         2363863228d-1*gammaaa*t11*t12*t26*t34+3.598504727726456d-
900     <         2*gammaaa*t15*t23*t34+1.799252363863228d-1*gammaaa*t15*t2
901     =         1*t22*t34+1.4394018910905823d-1*gammaaa*t15*t19*t20*t34+1
902     >         .0795514183179368d-1*gammaaa*t15*t17*t18*t34+7.1970094554
903     ?         52912d-2*gammaaa*t14*t15*t16*t34-1.9263326690781307d-4*t1
904     @         3*t30*t32*t33-9.631663345390654d-4*t11*t12*t30*t32*t33+1.
905     1         9263326690781307d-4*t2*t30*t31*t32-9.631663345390653d-5*t
906     2         15*t23*t30*t32-4.815831672695327d-4*t15*t21*t22*t30*t32-3
907     3         .8526653381562614d-4*t15*t19*t20*t30*t32-2.88949900361719
908     4         64d-4*t15*t17*t18*t30*t32-1.9263326690781307d-4*t14*t15*t
909     5         16*t30*t32)*wght-2.4814019635976003d+0*t25*t28*wght-4.135
910     6         669939329334d-1*t24*wght/rhoa**6.666666666666666d-1+Amat2
911     7         (iq,D2_RA_RA)
912            Cmat2(iq,D2_RA_GAA) = -9.305257363491002d-1*t1*(2.8894990036
913     1         171964d-4*gammaaa*t10*t33*t36*t9-7.22374750904299d-4*gamm
914     2         aaa*t12*t31*t36*t9-3.925641521156134d-2*t10*t26*t27*t9+2.
915     3         167124252712897d-4*gammaaa*t33*t36*t7*t8-2.94423114086710
916     4         07d-2*t26*t27*t7*t8-2.167124252712897d-4*gammaaa*t31*t36*
917     5         t6*t8-4.334248505425794d-4*gammaaa*t10*t31*t36*t7+1.44474
918     6         95018085982d-4*gammaaa*t2*t33*t36*t6-1.962820760578067d-2
919     7         *t2*t26*t27*t6-7.22374750904299d-4*gammaaa*t19*t22*t35*t3
920     8         6-4.334248505425794d-4*gammaaa*t17*t20*t35*t36-2.16712425
921     9         2712897d-4*gammaaa*t16*t18*t35*t36-7.22374750904299d-5*ga
922     :         mmaaa*t14*t35*t36+7.22374750904299d-5*gammaaa*t13*t33*t36
923     ;         +3.611873754521495d-4*gammaaa*t11*t12*t33*t36-7.223747509
924     <         04299d-5*gammaaa*t2*t31*t36+3.611873754521495d-5*gammaaa*
925     =         t15*t23*t36+1.8059368772607476d-4*gammaaa*t15*t21*t22*t36
926     >         +1.4447495018085982d-4*gammaaa*t15*t19*t20*t36+1.08356212
927     ?         63564485d-4*gammaaa*t15*t17*t18*t36+7.22374750904299d-5*g
928     @         ammaaa*t14*t15*t16*t36-9.814103802890335d-3*t13*t26*t27-4
929     1         .9070519014451675d-2*t11*t12*t26*t27-9.814103802890335d-3
930     2         *t15*t23*t27-4.9070519014451675d-2*t15*t21*t22*t27-3.9256
931     3         41521156134d-2*t15*t19*t20*t27-2.9442311408671007d-2*t15*
932     4         t17*t18*t27-1.962820760578067d-2*t14*t15*t16*t27)*wght-1.
933     5         2407009817988002d+0*t25*(1.4721155704335503d-2*t10*t26*t3
934     6         *t9+1.1040866778251626d-2*t26*t3*t7*t8+7.360577852167751d
935     7         -3*t2*t26*t3*t6+3.6802889260838756d-3*t13*t26*t3+1.840144
936     8         4630419378d-2*t11*t12*t26*t3+3.6802889260838756d-3*t15*t2
937     9         3*t3+1.8401444630419378d-2*t15*t21*t22*t3+1.4721155704335
938     :         503d-2*t15*t19*t20*t3+1.1040866778251626d-2*t15*t17*t18*t
939     ;         3+7.360577852167751d-3*t14*t15*t16*t3)*wght+Cmat2(iq,D2_R
940     <         A_GAA)
941            Cmat2(iq,D2_GAA_GAA) = Cmat2(iq,D2_GAA_GAA)-9.30525736349100
942     1         2d-1*t1*(t26*t3*(2.708905315891122d-4*t12*t26*t3*t9+8.126
943     2         715947673364d-5*t26*t3*t6*t8+1.625343189534673d-4*t10*t26
944     3         *t3*t7+2.7089053158911214d-5*t2*t26*t3)+t3*(2.70890531589
945     4         1122d-4*t19*t22*t3*t35+1.625343189534673d-4*t17*t20*t3*t3
946     5         5+8.126715947673364d-5*t16*t18*t3*t35+2.7089053158911214d
947     6         -5*t14*t3*t35-1.3544526579455607d-5*t15*t23*t3-6.77226328
948     7         9727805d-5*t15*t21*t22*t3-5.417810631782243d-5*t15*t19*t2
949     8         0*t3-4.063357973836682d-5*t15*t17*t18*t3-2.70890531589112
950     9         14d-5*t14*t15*t16*t3)-7.360577852167751d-3*t29*t33/rhoa**
951     :         5.333333333333333d+0)*wght
952          elseif (rhoa.le.tol_rho.and.rhob.gt.tol_rho) then
953            t1 = rhob**1.3333333333333333d+0
954            t2 = param(3)
955            t3 = 1/rhob**2.6666666666666666d+0
956            t4 = 3.6802889260838756d-3*gammabb*t3
957            t5 = t4+1.0d+0
958            t6 = 1.0d+0-1.0d+0/t5
959            t7 = t6**2.0d+0
960            t8 = param(4)
961            t9 = t6**3.0d+0
962            t10 = param(5)
963            t11 = t6**4.0d+0
964            t12 = param(6)
965            t13 = param(2)
966            t14 = param(9)
967            t15 = exp(-t4)
968            t16 = 1.0d+0-t15
969            t17 = t16**2.0d+0
970            t18 = param(10)
971            t19 = t16**3.0d+0
972            t20 = param(11)
973            t21 = t16**4.0d+0
974            t22 = param(12)
975            t23 = param(8)
976            t24 = t8*t9+t2*t7+t12*t6**5.0d+0+t13*t6+t16*t23+t16**5.0d+0*
977     1         t22+t20*t21+t18*t19+t14*t17+t10*t11+param(7)+param(1)
978            t25 = rhob**3.333333333333333d-1
979            t26 = 1/t5**2
980            t27 = 1/rhob**3.6666666666666664d+0
981            t28 = -3.925641521156134d-2*gammabb*t10*t26*t27*t9-2.9442311
982     1         408671007d-2*gammabb*t26*t27*t7*t8-1.962820760578067d-2*g
983     2         ammabb*t2*t26*t27*t6-9.814103802890335d-3*gammabb*t13*t26
984     3         *t27-4.9070519014451675d-2*gammabb*t11*t12*t26*t27-9.8141
985     4         03802890335d-3*gammabb*t15*t23*t27-4.9070519014451675d-2*
986     5         gammabb*t15*t21*t22*t27-3.925641521156134d-2*gammabb*t15*
987     6         t19*t20*t27-2.9442311408671007d-2*gammabb*t15*t17*t18*t27
988     7         -1.962820760578067d-2*gammabb*t14*t15*t16*t27
989            t29 = 1.4721155704335503d-2*t10*t9+1.1040866778251626d-2*t7*
990     1         t8+7.360577852167751d-3*t2*t6+3.6802889260838756d-3*t13+1
991     2         .8401444630419378d-2*t11*t12
992            t30 = gammabb**2
993            t31 = 1/t5**4
994            t32 = 1/rhob**7.333333333333333d+0
995            t33 = 1/t5**3
996            t34 = 1/rhob**4.666666666666667d+0
997            t35 = exp(-7.360577852167751d-3*gammabb*t3)
998            t36 = 1/rhob**6.333333333333333d+0
999            fnc(iq) = fnc(iq)-9.305257363491002d-1*t1*t24*wght
1000            Amat(iq,D1_RB) = -9.305257363491002d-1*t1*t28*wght-1.2407009
1001     1         817988002d+0*t24*t25*wght+Amat(iq,D1_RB)
1002            Cmat(iq,D1_GBB) = Cmat(iq,D1_GBB)-9.305257363491002d-1*t1*(t
1003     1         26*t29*t3+(3.6802889260838756d-3*t15*t23+1.84014446304193
1004     2         78d-2*t15*t21*t22+1.4721155704335503d-2*t15*t19*t20+1.104
1005     3         0866778251626d-2*t15*t17*t18+7.360577852167751d-3*t14*t15
1006     4         *t16)*t3)*wght
1007            Amat2(iq,D2_RB_RB) = -9.305257363491002d-1*t1*(1.43940189109
1008     1         05823d-1*gammabb*t10*t26*t34*t9-7.705330676312523d-4*t10*
1009     2         t30*t32*t33*t9+1.9263326690781307d-3*t12*t30*t31*t32*t9+1
1010     3         .0795514183179368d-1*gammabb*t26*t34*t7*t8-5.778998007234
1011     4         393d-4*t30*t32*t33*t7*t8+5.778998007234393d-4*t30*t31*t32
1012     5         *t6*t8+1.1557996014468785d-3*t10*t30*t31*t32*t7+7.1970094
1013     6         55452912d-2*gammabb*t2*t26*t34*t6-3.8526653381562614d-4*t
1014     7         2*t30*t32*t33*t6+1.9263326690781307d-3*t19*t22*t30*t32*t3
1015     8         5+1.1557996014468785d-3*t17*t20*t30*t32*t35+5.77899800723
1016     9         4393d-4*t16*t18*t30*t32*t35+1.9263326690781307d-4*t14*t30
1017     :         *t32*t35+3.598504727726456d-2*gammabb*t13*t26*t34+1.79925
1018     ;         2363863228d-1*gammabb*t11*t12*t26*t34+3.598504727726456d-
1019     <         2*gammabb*t15*t23*t34+1.799252363863228d-1*gammabb*t15*t2
1020     =         1*t22*t34+1.4394018910905823d-1*gammabb*t15*t19*t20*t34+1
1021     >         .0795514183179368d-1*gammabb*t15*t17*t18*t34+7.1970094554
1022     ?         52912d-2*gammabb*t14*t15*t16*t34-1.9263326690781307d-4*t1
1023     @         3*t30*t32*t33-9.631663345390654d-4*t11*t12*t30*t32*t33+1.
1024     1         9263326690781307d-4*t2*t30*t31*t32-9.631663345390653d-5*t
1025     2         15*t23*t30*t32-4.815831672695327d-4*t15*t21*t22*t30*t32-3
1026     3         .8526653381562614d-4*t15*t19*t20*t30*t32-2.88949900361719
1027     4         64d-4*t15*t17*t18*t30*t32-1.9263326690781307d-4*t14*t15*t
1028     5         16*t30*t32)*wght-2.4814019635976003d+0*t25*t28*wght-4.135
1029     6         669939329334d-1*t24*wght/rhob**6.666666666666666d-1+Amat2
1030     7         (iq,D2_RB_RB)
1031            Cmat2(iq,D2_RB_GBB) = -9.305257363491002d-1*t1*(2.8894990036
1032     1         171964d-4*gammabb*t10*t33*t36*t9-7.22374750904299d-4*gamm
1033     2         abb*t12*t31*t36*t9-3.925641521156134d-2*t10*t26*t27*t9+2.
1034     3         167124252712897d-4*gammabb*t33*t36*t7*t8-2.94423114086710
1035     4         07d-2*t26*t27*t7*t8-2.167124252712897d-4*gammabb*t31*t36*
1036     5         t6*t8-4.334248505425794d-4*gammabb*t10*t31*t36*t7+1.44474
1037     6         95018085982d-4*gammabb*t2*t33*t36*t6-1.962820760578067d-2
1038     7         *t2*t26*t27*t6-7.22374750904299d-4*gammabb*t19*t22*t35*t3
1039     8         6-4.334248505425794d-4*gammabb*t17*t20*t35*t36-2.16712425
1040     9         2712897d-4*gammabb*t16*t18*t35*t36-7.22374750904299d-5*ga
1041     :         mmabb*t14*t35*t36+7.22374750904299d-5*gammabb*t13*t33*t36
1042     ;         +3.611873754521495d-4*gammabb*t11*t12*t33*t36-7.223747509
1043     <         04299d-5*gammabb*t2*t31*t36+3.611873754521495d-5*gammabb*
1044     =         t15*t23*t36+1.8059368772607476d-4*gammabb*t15*t21*t22*t36
1045     >         +1.4447495018085982d-4*gammabb*t15*t19*t20*t36+1.08356212
1046     ?         63564485d-4*gammabb*t15*t17*t18*t36+7.22374750904299d-5*g
1047     @         ammabb*t14*t15*t16*t36-9.814103802890335d-3*t13*t26*t27-4
1048     1         .9070519014451675d-2*t11*t12*t26*t27-9.814103802890335d-3
1049     2         *t15*t23*t27-4.9070519014451675d-2*t15*t21*t22*t27-3.9256
1050     3         41521156134d-2*t15*t19*t20*t27-2.9442311408671007d-2*t15*
1051     4         t17*t18*t27-1.962820760578067d-2*t14*t15*t16*t27)*wght-1.
1052     5         2407009817988002d+0*t25*(1.4721155704335503d-2*t10*t26*t3
1053     6         *t9+1.1040866778251626d-2*t26*t3*t7*t8+7.360577852167751d
1054     7         -3*t2*t26*t3*t6+3.6802889260838756d-3*t13*t26*t3+1.840144
1055     8         4630419378d-2*t11*t12*t26*t3+3.6802889260838756d-3*t15*t2
1056     9         3*t3+1.8401444630419378d-2*t15*t21*t22*t3+1.4721155704335
1057     :         503d-2*t15*t19*t20*t3+1.1040866778251626d-2*t15*t17*t18*t
1058     ;         3+7.360577852167751d-3*t14*t15*t16*t3)*wght+Cmat2(iq,D2_R
1059     <         B_GBB)
1060            Cmat2(iq,D2_GBB_GBB) = Cmat2(iq,D2_GBB_GBB)-9.30525736349100
1061     1         2d-1*t1*(t26*t3*(2.708905315891122d-4*t12*t26*t3*t9+8.126
1062     2         715947673364d-5*t26*t3*t6*t8+1.625343189534673d-4*t10*t26
1063     3         *t3*t7+2.7089053158911214d-5*t2*t26*t3)+t3*(2.70890531589
1064     4         1122d-4*t19*t22*t3*t35+1.625343189534673d-4*t17*t20*t3*t3
1065     5         5+8.126715947673364d-5*t16*t18*t3*t35+2.7089053158911214d
1066     6         -5*t14*t3*t35-1.3544526579455607d-5*t15*t23*t3-6.77226328
1067     7         9727805d-5*t15*t21*t22*t3-5.417810631782243d-5*t15*t19*t2
1068     8         0*t3-4.063357973836682d-5*t15*t17*t18*t3-2.70890531589112
1069     9         14d-5*t14*t15*t16*t3)-7.360577852167751d-3*t29*t33/rhob**
1070     :         5.333333333333333d+0)*wght
1071          endif ! rhoa.gt.tol_rho.and.rhob.gt.tol_rho
1072        endif ! ipol.eq.1
1073      enddo ! iq
1074      end
1075C>
1076C> \brief Evaluate the nwxcm_x_sogga functional [1]
1077C>
1078C> \f{eqnarray*}{
1079C>   {\it t_1} &=& {\it param}\left(1\right)\\\\
1080C>   {\it t_2} &=& {\it param}\left(7\right)\\\\
1081C>   {\it t_3} &=& {\it param}\left(2\right)\\\\
1082C>   {\it t_4} &=& {{0.003680288926083876\,
1083C>    \sigma_{\alpha\alpha}}\over{\rho_\alpha^{{{8}\over{3}}}}}\\\\
1084C>   {\it t_5} &=& 1.0-{{1.0}\over{{\it t_4}+1.0}}\\\\
1085C>   {\it t_6} &=& {\it param}\left(3\right)\\\\
1086C>   {\it t_7} &=& {\it param}\left(4\right)\\\\
1087C>   {\it t_8} &=& {\it param}\left(5\right)\\\\
1088C>   {\it t_9} &=& {\it param}\left(6\right)\\\\
1089C>   {\it t_{10}} &=& {\it param}\left(9\right)\\\\
1090C>   {\it t_{11}} &=& 1.0-{{1}\over{e^{{\it t_4}}}}\\\\
1091C>   {\it t_{12}} &=& {\it param}\left(10\right)\\\\
1092C>   {\it t_{13}} &=& {\it param}\left(11\right)\\\\
1093C>   {\it t_{14}} &=& {\it param}\left(12\right)\\\\
1094C>   {\it t_{15}} &=& {\it param}\left(8\right)\\\\
1095C>   {\it t_{16}} &=& {{0.003680288926083876\,
1096C>    \sigma_{\beta\beta}}\over{\rho_\beta^{{{8}\over{3}}}}}\\\\
1097C>   {\it t_{17}} &=& 1.0-{{1.0}\over{{\it t_{16}}+1.0}}\\\\
1098C>   {\it t_{18}} &=& 1.0-{{1}\over{e^{{\it t_{16}}}}}\\\\
1099C>   {\it t_{19}} &=& {{0.003680288926083876\,\sigma_{ss}}
1100C>    \over{\rho_s^{{{8}\over{3}}}}}\\\\
1101C>   {\it t_{20}} &=& 1.0-{{1.0}\over{{\it t_{19}}+1.0}}\\\\
1102C>   {\it t_{21}} &=& 1.0-{{1}\over{e^{{\it t_{19}}}}}\\\\
1103C>   f &=& -0.9305257363491002\,\rho_\beta^{{{4}\over{3}}}\,
1104C>    \left({\it t_{15}}\,{\it t_{18}}+{\it t_{14}}\,{
1105C>    \it t_{18}}^{5.0}+{\it t_{13}}\,{\it t_{18}}^{4.0}+{
1106C>    \it t_{12}}\,{\it t_{18}}^{3.0}+{\it t_{10}}\,{\it t_{18}}^{2.0}
1107C>    +{\it t_9}\,{\it t_{17}}^{5.0}+{\it t_8}\,{\it t_{17}}^{4.0}
1108C>    +{\it t_7}\,{\it t_{17}}^{3.0}+{\it t_6}\,{\it t_{17}}^{2.0}
1109C>    +{\it t_3}\,{\it t_{17}}+{\it t_2}+{\it t_1}\right)
1110C>    -0.9305257363491002\,\rho_\alpha^{{{4}\over{3}}}\,\left({
1111C>    \it t_{15}}\,{\it t_{11}}+{\it t_{14}}\,{\it t_{11}}^{5.0}+{
1112C>    \it t_{13}}\,{\it t_{11}}^{4.0}+{\it t_{12}}\,{\it t_{11}}^{3.0}
1113C>    +{\it t_{10}}\,{\it t_{11}}^{2.0}+{\it t_9}\,{\it t_5}^{5.0}
1114C>    +{\it t_8}\,{\it t_5}^{4.0}+{\it t_7}\,{\it t_5}^{3.0}+{
1115C>    \it t_6}\,{\it t_5}^{2.0}+{\it t_3}\,{\it t_5}+{\it t_2}+{
1116C>    \it t_1}\right)\\\\
1117C>   g &=& 0\\\\
1118C>   G &=& -0.9305257363491002\,\rho_s^{{{4}\over{3}}}\,\left({
1119C>    \it t_{15}}\,{\it t_{21}}+{\it t_{14}}\,{\it t_{21}}^{5.0}+{
1120C>    \it t_{13}}\,{\it t_{21}}^{4.0}+{\it t_{12}}\,{\it t_{21}}^{3.0}
1121C>    +{\it t_{10}}\,{\it t_{21}}^{2.0}+{\it t_9}\,{\it t_{20}}^{5.0}
1122C>    +{\it t_8}\,{\it t_{20}}^{4.0}+{\it t_7}\,{\it t_{20}}^{3.0}
1123C>    +{\it t_6}\,{\it t_{20}}^{2.0}+{\it t_3}\,{\it t_{20}}+{\it t_2}
1124C>    +{\it t_1}\right)\\\\
1125C> \f}
1126C>
1127C> Code generated with Maxima 5.34.0 [2]
1128C> driven by autoxc [3].
1129C>
1130C> ### References ###
1131C>
1132C> [1] Y Zhao, DG Truhlar, J.Chem.Phys. 128, 184109 (2008)  , DOI:
1133C> <a href="https://doi.org/10.1063/1.2912068 ">
1134C> 10.1063/1.2912068 </a>
1135C>
1136C> [2] Maxima, a computer algebra system,
1137C> <a href="http://maxima.sourceforge.net/">
1138C> http://maxima.sourceforge.net/</a>
1139C>
1140C> [3] autoxc, revision 27097 2015-05-08
1141C>
1142      subroutine nwxcm_x_sogga_d3(param,tol_rho,ipol,nq,wght,
1143     +rho,rgamma,fnc,Amat,Amat2,Amat3,
1144     +Cmat,Cmat2,Cmat3)
1145c $Id: $
1146#ifdef NWXC_QUAD_PREC
1147      implicit real(kind=selected_real_kind(30))(a-h,o-z),integer(i-n)
1148      integer, parameter :: rk=selected_real_kind(30)
1149#else
1150      implicit real(kind=selected_real_kind(15))(a-h,o-z),integer(i-n)
1151      integer, parameter :: rk=selected_real_kind(15)
1152#endif
1153#include "nwxc_param.fh"
1154      double precision param(*)     !< [Input] Parameters of functional
1155      double precision tol_rho      !< [Input] The lower limit on the density
1156      integer ipol                  !< [Input] The number of spin channels
1157      integer nq                    !< [Input] The number of points
1158      double precision wght         !< [Input] The weight of the functional
1159      double precision rho(nq,*)    !< [Input] The density
1160      double precision rgamma(nq,*) !< [Input] The norm of the density
1161                                    !< gradients
1162      double precision fnc(nq)      !< [Output] The value of the functional
1163c
1164c     Sampling Matrices for the XC Kernel
1165c
1166      double precision Amat(nq,*)   !< [Output] The derivative wrt rho
1167      double precision Cmat(nq,*)   !< [Output] The derivative wrt rgamma
1168c
1169c     Sampling Matrices for the XC Kernel
1170c
1171      double precision Amat2(nq,*)  !< [Output] The 2nd derivative wrt rho
1172      double precision Cmat2(nq,*)  !< [Output] The 2nd derivative wrt rgamma
1173                                    !< and possibly rho
1174c
1175c     Sampling Matrices for the XC Kernel
1176c
1177      double precision Amat3(nq,*)  !< [Output] The 3rd derivative wrt rho
1178      double precision Cmat3(nq,*)  !< [Output] The 3rd derivative wrt rgamma
1179                                    !< and possibly rho
1180      integer iq
1181      double precision tmp
1182      double precision rhoa,rhob
1183      double precision gammaaa,gammaab,gammabb
1184      double precision taua,taub
1185      double precision nwxcm_heaviside
1186      external         nwxcm_heaviside
1187CDIR$ NOVECTOR
1188      do iq = 1, nq
1189        if (ipol.eq.1) then
1190          rhoa    = 0.5d0*rho(iq,R_T)
1191          gammaaa = 0.25d0*rgamma(iq,G_TT)
1192          if (rhoa.gt.tol_rho) then
1193            t1 = rhoa**1.3333333333333333d+0
1194            t2 = param(3)
1195            t3 = 1/rhoa**2.6666666666666666d+0
1196            t4 = 3.6802889260838756d-3*gammaaa*t3
1197            t5 = t4+1.0d+0
1198            t6 = 1.0d+0-1.0d+0/t5
1199            t7 = t6**2.0d+0
1200            t8 = param(4)
1201            t9 = t6**3.0d+0
1202            t10 = param(5)
1203            t11 = t6**4.0d+0
1204            t12 = param(6)
1205            t13 = param(2)
1206            t14 = param(9)
1207            t15 = exp(-t4)
1208            t16 = 1.0d+0-t15
1209            t17 = t16**2.0d+0
1210            t18 = param(10)
1211            t19 = t16**3.0d+0
1212            t20 = param(11)
1213            t21 = t16**4.0d+0
1214            t22 = param(12)
1215            t23 = param(8)
1216            t24 = t8*t9+t2*t7+t12*t6**5.0d+0+t13*t6+t16*t23+t16**5.0d+0*
1217     1         t22+t20*t21+t18*t19+t14*t17+t10*t11+param(7)+param(1)
1218            t25 = rhoa**3.333333333333333d-1
1219            t26 = 1/t5**2
1220            t27 = 1/rhoa**3.6666666666666664d+0
1221            t28 = -3.925641521156134d-2*gammaaa*t10*t26*t27*t9-2.9442311
1222     1         408671007d-2*gammaaa*t26*t27*t7*t8-1.962820760578067d-2*g
1223     2         ammaaa*t2*t26*t27*t6-9.814103802890335d-3*gammaaa*t13*t26
1224     3         *t27-4.9070519014451675d-2*gammaaa*t11*t12*t26*t27-9.8141
1225     4         03802890335d-3*gammaaa*t15*t23*t27-4.9070519014451675d-2*
1226     5         gammaaa*t15*t21*t22*t27-3.925641521156134d-2*gammaaa*t15*
1227     6         t19*t20*t27-2.9442311408671007d-2*gammaaa*t15*t17*t18*t27
1228     7         -1.962820760578067d-2*gammaaa*t14*t15*t16*t27
1229            t29 = 1.4721155704335503d-2*t10*t9+1.1040866778251626d-2*t7*
1230     1         t8+7.360577852167751d-3*t2*t6+3.6802889260838756d-3*t13+1
1231     2         .8401444630419378d-2*t11*t12
1232            t30 = 1/rhoa**6.666666666666666d-1
1233            t31 = gammaaa**2
1234            t32 = 1/t5**4
1235            t33 = 1/rhoa**7.333333333333333d+0
1236            t34 = 1/t5**3
1237            t35 = 1/rhoa**4.666666666666667d+0
1238            t36 = exp(-7.360577852167751d-3*gammaaa*t3)
1239            t37 = 1.4394018910905823d-1*gammaaa*t10*t26*t35*t9-7.7053306
1240     1         76312523d-4*t10*t31*t33*t34*t9+1.9263326690781307d-3*t12*
1241     2         t31*t32*t33*t9+1.0795514183179368d-1*gammaaa*t26*t35*t7*t
1242     3         8-5.778998007234393d-4*t31*t33*t34*t7*t8+5.77899800723439
1243     4         3d-4*t31*t32*t33*t6*t8+1.1557996014468785d-3*t10*t31*t32*
1244     5         t33*t7+7.197009455452912d-2*gammaaa*t2*t26*t35*t6-3.85266
1245     6         53381562614d-4*t2*t31*t33*t34*t6+1.9263326690781307d-3*t1
1246     7         9*t22*t31*t33*t36+1.1557996014468785d-3*t17*t20*t31*t33*t
1247     8         36+5.778998007234393d-4*t16*t18*t31*t33*t36+1.92633266907
1248     9         81307d-4*t14*t31*t33*t36+3.598504727726456d-2*gammaaa*t13
1249     :         *t26*t35+1.799252363863228d-1*gammaaa*t11*t12*t26*t35+3.5
1250     ;         98504727726456d-2*gammaaa*t15*t23*t35+1.799252363863228d-
1251     <         1*gammaaa*t15*t21*t22*t35+1.4394018910905823d-1*gammaaa*t
1252     =         15*t19*t20*t35+1.0795514183179368d-1*gammaaa*t15*t17*t18*
1253     >         t35+7.197009455452912d-2*gammaaa*t14*t15*t16*t35-1.926332
1254     ?         6690781307d-4*t13*t31*t33*t34-9.631663345390654d-4*t11*t1
1255     @         2*t31*t33*t34+1.9263326690781307d-4*t2*t31*t32*t33-9.6316
1256     1         63345390653d-5*t15*t23*t31*t33-4.815831672695327d-4*t15*t
1257     2         21*t22*t31*t33-3.8526653381562614d-4*t15*t19*t20*t31*t33-
1258     3         2.8894990036171964d-4*t15*t17*t18*t31*t33-1.9263326690781
1259     4         307d-4*t14*t15*t16*t31*t33
1260            t38 = 1/rhoa**6.333333333333333d+0
1261            t39 = 2.8894990036171964d-4*gammaaa*t10*t34*t38*t9-7.2237475
1262     1         0904299d-4*gammaaa*t12*t32*t38*t9-3.925641521156134d-2*t1
1263     2         0*t26*t27*t9+2.167124252712897d-4*gammaaa*t34*t38*t7*t8-2
1264     3         .9442311408671007d-2*t26*t27*t7*t8-2.167124252712897d-4*g
1265     4         ammaaa*t32*t38*t6*t8-4.334248505425794d-4*gammaaa*t10*t32
1266     5         *t38*t7+1.4447495018085982d-4*gammaaa*t2*t34*t38*t6-1.962
1267     6         820760578067d-2*t2*t26*t27*t6-7.22374750904299d-4*gammaaa
1268     7         *t19*t22*t36*t38-4.334248505425794d-4*gammaaa*t17*t20*t36
1269     8         *t38-2.167124252712897d-4*gammaaa*t16*t18*t36*t38-7.22374
1270     9         750904299d-5*gammaaa*t14*t36*t38+7.22374750904299d-5*gamm
1271     :         aaa*t13*t34*t38+3.611873754521495d-4*gammaaa*t11*t12*t34*
1272     ;         t38-7.22374750904299d-5*gammaaa*t2*t32*t38+3.611873754521
1273     <         495d-5*gammaaa*t15*t23*t38+1.8059368772607476d-4*gammaaa*
1274     =         t15*t21*t22*t38+1.4447495018085982d-4*gammaaa*t15*t19*t20
1275     >         *t38+1.0835621263564485d-4*gammaaa*t15*t17*t18*t38+7.2237
1276     ?         4750904299d-5*gammaaa*t14*t15*t16*t38-9.814103802890335d-
1277     @         3*t13*t26*t27-4.9070519014451675d-2*t11*t12*t26*t27-9.814
1278     1         103802890335d-3*t15*t23*t27-4.9070519014451675d-2*t15*t21
1279     2         *t22*t27-3.925641521156134d-2*t15*t19*t20*t27-2.944231140
1280     3         8671007d-2*t15*t17*t18*t27-1.962820760578067d-2*t14*t15*t
1281     4         16*t27
1282            t40 = 1.4721155704335503d-2*t10*t26*t3*t9+1.1040866778251626
1283     1         d-2*t26*t3*t7*t8+7.360577852167751d-3*t2*t26*t3*t6+3.6802
1284     2         889260838756d-3*t13*t26*t3+1.8401444630419378d-2*t11*t12*
1285     3         t26*t3+3.6802889260838756d-3*t15*t23*t3+1.840144463041937
1286     4         8d-2*t15*t21*t22*t3+1.4721155704335503d-2*t15*t19*t20*t3+
1287     5         1.1040866778251626d-2*t15*t17*t18*t3+7.360577852167751d-3
1288     6         *t14*t15*t16*t3
1289            t41 = 1/rhoa**5.333333333333333d+0
1290            t42 = 2.708905315891122d-4*t12*t26*t3*t9+8.126715947673364d-
1291     1         5*t26*t3*t6*t8+1.625343189534673d-4*t10*t26*t3*t7+2.70890
1292     2         53158911214d-5*t2*t26*t3
1293            t43 = gammaaa**3
1294            t44 = 1/t5**6
1295            t45 = 1/rhoa**11
1296            t46 = 1/t5**5
1297            t47 = 1/rhoa**8.333333333333334d+0
1298            t48 = 1/rhoa**5.666666666666667d+0
1299            t49 = exp(-1.1040866778251626d-2*gammaaa*t3)
1300            t50 = 1/rhoa**10
1301            t51 = 1/rhoa**9
1302            fnc(iq) = fnc(iq)-1.8610514726982003d+0*t1*t24*wght
1303            Amat(iq,D1_RA) = (-9.305257363491002d-1*t1*t28-1.24070098179
1304     1         88002d+0*t24*t25)*wght+Amat(iq,D1_RA)
1305            Cmat(iq,D1_GAA) = Cmat(iq,D1_GAA)-9.305257363491002d-1*t1*(t
1306     1         26*t29*t3+(3.6802889260838756d-3*t15*t23+1.84014446304193
1307     2         78d-2*t15*t21*t22+1.4721155704335503d-2*t15*t19*t20+1.104
1308     3         0866778251626d-2*t15*t17*t18+7.360577852167751d-3*t14*t15
1309     4         *t16)*t3)*wght
1310            Cmat(iq,D1_GAB) = Cmat(iq,D1_GAB)
1311            Amat2(iq,D2_RA_RA) = (-9.305257363491002d-1*t1*t37-4.1356699
1312     1         39329334d-1*t24*t30-2.4814019635976003d+0*t25*t28)*wght+A
1313     2         mat2(iq,D2_RA_RA)
1314            Amat2(iq,D2_RA_RB) = Amat2(iq,D2_RA_RB)
1315            Cmat2(iq,D2_RA_GAA) = (-1.2407009817988002d+0*t25*t40-9.3052
1316     1         57363491002d-1*t1*t39)*wght+Cmat2(iq,D2_RA_GAA)
1317            Cmat2(iq,D2_RA_GAB) = Cmat2(iq,D2_RA_GAB)
1318            Cmat2(iq,D2_RA_GBB) = Cmat2(iq,D2_RA_GBB)
1319            Cmat2(iq,D2_GAA_GAA) = Cmat2(iq,D2_GAA_GAA)-9.30525736349100
1320     1         2d-1*t1*(t26*t3*t42-7.360577852167751d-3*t29*t34*t41+t3*(
1321     2         2.708905315891122d-4*t19*t22*t3*t36+1.625343189534673d-4*
1322     3         t17*t20*t3*t36+8.126715947673364d-5*t16*t18*t3*t36+2.7089
1323     4         053158911214d-5*t14*t3*t36-1.3544526579455607d-5*t15*t23*
1324     5         t3-6.772263289727805d-5*t15*t21*t22*t3-5.417810631782243d
1325     6         -5*t15*t19*t20*t3-4.063357973836682d-5*t15*t17*t18*t3-2.7
1326     7         089053158911214d-5*t14*t15*t16*t3))*wght
1327            Cmat2(iq,D2_GAA_GAB) = Cmat2(iq,D2_GAA_GAB)
1328            Cmat2(iq,D2_GAA_GBB) = Cmat2(iq,D2_GAA_GBB)
1329            Cmat2(iq,D2_GAB_GAB) = Cmat2(iq,D2_GAB_GAB)
1330            Amat3(iq,D3_RA_RA_RA) = (-9.305257363491002d-1*t1*(-6.717208
1331     1         825089385d-1*gammaaa*t10*t26*t48*t9+8.475863743943775d-3*
1332     2         t10*t31*t34*t47*t9-2.1189659359859436d-2*t12*t31*t32*t47*
1333     3         t9+1.1343137263938945d-4*t12*t43*t45*t46*t9-2.26862745278
1334     4         77887d-5*t10*t32*t43*t45*t9-5.037906618817038d-1*gammaaa*
1335     5         t26*t48*t7*t8+6.356897807957831d-3*t31*t34*t47*t7*t8-1.70
1336     6         14705895908416d-5*t32*t43*t45*t7*t8-6.356897807957831d-3*
1337     7         t31*t32*t47*t6*t8+3.4029411791816827d-5*t43*t45*t46*t6*t8
1338     8         -5.671568631969471d-6*t43*t44*t45*t8-1.2713795615915663d-
1339     9         2*t10*t31*t32*t47*t7+6.805882358363365d-5*t10*t43*t45*t46
1340     :         *t7-5.671568631969472d-5*t12*t43*t44*t45*t7-3.35860441254
1341     ;         46923d-1*gammaaa*t2*t26*t48*t6+4.237931871971887d-3*t2*t3
1342     <         1*t34*t47*t6-2.2686274527877887d-5*t10*t43*t44*t45*t6-1.1
1343     =         343137263938943d-5*t2*t32*t43*t45*t6-5.671568631969472d-5
1344     >         *t17*t22*t43*t45*t49-2.2686274527877887d-5*t16*t20*t43*t4
1345     ?         5*t49-5.671568631969471d-6*t18*t43*t45*t49-1.679302206272
1346     @         3462d-1*gammaaa*t13*t26*t48-8.396511031361731d-1*gammaaa*
1347     1         t11*t12*t26*t48-1.6793022062723462d-1*gammaaa*t15*t23*t48
1348     2         -8.396511031361731d-1*gammaaa*t15*t21*t22*t48-6.717208825
1349     3         089385d-1*gammaaa*t15*t19*t20*t48-5.037906618817038d-1*ga
1350     4         mmaaa*t15*t17*t18*t48-3.3586044125446923d-1*gammaaa*t14*t
1351     5         15*t16*t48-2.1189659359859436d-2*t19*t22*t31*t36*t47-1.27
1352     6         13795615915663d-2*t17*t20*t31*t36*t47-6.356897807957831d-
1353     7         3*t16*t18*t31*t36*t47-2.1189659359859436d-3*t14*t31*t36*t
1354     8         47+2.1189659359859436d-3*t13*t31*t34*t47+1.05948296799297
1355     9         18d-2*t11*t12*t31*t34*t47-2.1189659359859436d-3*t2*t31*t3
1356     :         2*t47+1.0594829679929718d-3*t15*t23*t31*t47+5.29741483996
1357     ;         4859d-3*t15*t21*t22*t31*t47+4.237931871971887d-3*t15*t19*
1358     <         t20*t31*t47+3.1784489039789154d-3*t15*t17*t18*t31*t47+2.1
1359     =         189659359859436d-3*t14*t15*t16*t31*t47+1.1343137263938943
1360     >         d-5*t2*t43*t45*t46+5.671568631969472d-5*t19*t22*t36*t43*t
1361     ?         45+3.4029411791816827d-5*t17*t20*t36*t43*t45+1.7014705895
1362     @         908414d-5*t16*t18*t36*t43*t45+5.671568631969471d-6*t14*t3
1363     1         6*t43*t45-5.671568631969471d-6*t13*t32*t43*t45-2.83578431
1364     2         5984736d-5*t11*t12*t32*t43*t45-9.452614386615786d-7*t15*t
1365     3         23*t43*t45-4.726307193307893d-6*t15*t21*t22*t43*t45-3.781
1366     4         0457546463144d-6*t15*t19*t20*t43*t45-2.8357843159847357d-
1367     5         6*t15*t17*t18*t43*t45-1.8905228773231572d-6*t14*t15*t16*t
1368     6         43*t45)-3.7221029453964005d+0*t25*t37-1.2407009817988002d
1369     7         +0*t28*t30+2.7571132928862224d-1*t24/rhoa**1.666666666666
1370     8         6669d+0)*wght+Amat3(iq,D3_RA_RA_RA)
1371            Amat3(iq,D3_RA_RA_RB) = Amat3(iq,D3_RA_RA_RB)
1372            Cmat3(iq,D3_RA_RA_GAA) = (-9.305257363491002d-1*t1*(-4.25367
1373     1         64739771044d-5*t12*t31*t46*t50*t9+8.507352947954206d-6*t1
1374     2         0*t31*t32*t50*t9+1.4394018910905823d-1*t10*t26*t35*t9-2.6
1375     3         00549103255476d-3*gammaaa*t10*t33*t34*t9+6.50137275813869
1376     4         d-3*gammaaa*t12*t32*t33*t9+6.380514710965656d-6*t31*t32*t
1377     5         50*t7*t8+1.0795514183179368d-1*t26*t35*t7*t8-1.9504118274
1378     6         416074d-3*gammaaa*t33*t34*t7*t8-1.2761029421931314d-5*t31
1379     7         *t46*t50*t6*t8+1.9504118274416074d-3*gammaaa*t32*t33*t6*t
1380     8         8+2.1268382369885522d-6*t31*t44*t50*t8-2.552205884386262d
1381     9         -5*t10*t31*t46*t50*t7+2.1268382369885516d-5*t12*t31*t44*t
1382     :         50*t7+3.900823654883215d-3*gammaaa*t10*t32*t33*t7+8.50735
1383     ;         2947954209d-6*t10*t31*t44*t50*t6+4.253676473977103d-6*t2*
1384     <         t31*t32*t50*t6+7.197009455452912d-2*t2*t26*t35*t6-1.30027
1385     =         4551627738d-3*gammaaa*t2*t33*t34*t6+2.1268382369885516d-5
1386     >         *t17*t22*t31*t49*t50+8.507352947954209d-6*t16*t20*t31*t49
1387     ?         *t50+2.1268382369885522d-6*t18*t31*t49*t50-4.253676473977
1388     @         103d-6*t2*t31*t46*t50-2.1268382369885522d-5*t19*t22*t31*t
1389     1         36*t50-1.276102942193131d-5*t17*t20*t31*t36*t50-6.3805147
1390     2         10965657d-6*t16*t18*t31*t36*t50-2.1268382369885516d-6*t14
1391     3         *t31*t36*t50+2.1268382369885516d-6*t13*t31*t32*t50+1.0634
1392     4         191184942758d-5*t11*t12*t31*t32*t50+3.54473039498092d-7*t
1393     5         15*t23*t31*t50+1.77236519749046d-6*t15*t21*t22*t31*t50+1.
1394     6         4178921579923678d-6*t15*t19*t20*t31*t50+1.063419118494276
1395     7         1d-6*t15*t17*t18*t31*t50+7.08946078996184d-7*t14*t15*t16*
1396     8         t31*t50+6.50137275813869d-3*gammaaa*t19*t22*t33*t36+3.900
1397     9         823654883215d-3*gammaaa*t17*t20*t33*t36+1.950411827441607
1398     :         4d-3*gammaaa*t16*t18*t33*t36+6.501372758138692d-4*gammaaa
1399     ;         *t14*t33*t36+3.598504727726456d-2*t13*t26*t35+1.799252363
1400     <         863228d-1*t11*t12*t26*t35+3.598504727726456d-2*t15*t23*t3
1401     =         5+1.799252363863228d-1*t15*t21*t22*t35+1.4394018910905823
1402     >         d-1*t15*t19*t20*t35+1.0795514183179368d-1*t15*t17*t18*t35
1403     ?         +7.197009455452912d-2*t14*t15*t16*t35-6.501372758138692d-
1404     @         4*gammaaa*t13*t33*t34-3.250686379069345d-3*gammaaa*t11*t1
1405     1         2*t33*t34+6.501372758138692d-4*gammaaa*t2*t32*t33-3.25068
1406     2         6379069346d-4*gammaaa*t15*t23*t33-1.6253431895346726d-3*g
1407     3         ammaaa*t15*t21*t22*t33-1.300274551627738d-3*gammaaa*t15*t
1408     4         19*t20*t33-9.752059137208037d-4*gammaaa*t15*t17*t18*t33-6
1409     5         .501372758138692d-4*gammaaa*t14*t15*t16*t33)-4.1356699393
1410     6         29334d-1*t30*t40-2.4814019635976003d+0*t25*t39)*wght+Cmat
1411     7         3(iq,D3_RA_RA_GAA)
1412            Cmat3(iq,D3_RA_RA_GAB) = Cmat3(iq,D3_RA_RA_GAB)
1413            Cmat3(iq,D3_RA_RA_GBB) = Cmat3(iq,D3_RA_RA_GBB)
1414            Cmat3(iq,D3_RA_RB_GAA) = Cmat3(iq,D3_RA_RB_GAA)
1415            Cmat3(iq,D3_RA_RB_GAB) = Cmat3(iq,D3_RA_RB_GAB)
1416            Cmat3(iq,D3_RA_GAA_GAA) = (-9.305257363491002d-1*t1*(1.59512
1417     1         8677741414d-5*gammaaa*t12*t46*t51*t9-3.190257355482828d-6
1418     2         *gammaaa*t10*t32*t51*t9+5.778998007234393d-4*t10*t34*t38*
1419     3         t9-1.4447495018085982d-3*t12*t32*t38*t9-2.392693016612120
1420     4         7d-6*gammaaa*t32*t51*t7*t8+4.334248505425794d-4*t34*t38*t
1421     5         7*t8+4.785386033224242d-6*gammaaa*t46*t51*t6*t8-4.3342485
1422     6         05425794d-4*t32*t38*t6*t8-7.97564338870707d-7*gammaaa*t44
1423     7         *t51*t8+9.570772066448484d-6*gammaaa*t10*t46*t51*t7-7.975
1424     8         64338870707d-6*gammaaa*t12*t44*t51*t7-8.668497010851588d-
1425     9         4*t10*t32*t38*t7-3.190257355482828d-6*gammaaa*t10*t44*t51
1426     :         *t6-1.595128677741414d-6*gammaaa*t2*t32*t51*t6+2.88949900
1427     ;         36171964d-4*t2*t34*t38*t6-7.97564338870707d-6*gammaaa*t17
1428     <         *t22*t49*t51-3.190257355482828d-6*gammaaa*t16*t20*t49*t51
1429     =         -7.97564338870707d-7*gammaaa*t18*t49*t51+1.59512867774141
1430     >         42d-6*gammaaa*t2*t46*t51+7.97564338870707d-6*gammaaa*t19*
1431     ?         t22*t36*t51+4.785386033224242d-6*gammaaa*t17*t20*t36*t51+
1432     @         2.392693016612121d-6*gammaaa*t16*t18*t36*t51+7.9756433887
1433     1         07071d-7*gammaaa*t14*t36*t51-7.97564338870707d-7*gammaaa*
1434     2         t13*t32*t51-3.987821694353535d-6*gammaaa*t11*t12*t32*t51-
1435     3         1.329273898117845d-7*gammaaa*t15*t23*t51-6.64636949058922
1436     4         6d-7*gammaaa*t15*t21*t22*t51-5.31709559247138d-7*gammaaa*
1437     5         t15*t19*t20*t51-3.987821694353535d-7*gammaaa*t15*t17*t18*
1438     6         t51-2.65854779623569d-7*gammaaa*t14*t15*t16*t51-1.4447495
1439     7         018085982d-3*t19*t22*t36*t38-8.668497010851588d-4*t17*t20
1440     8         *t36*t38-4.334248505425794d-4*t16*t18*t36*t38-1.444749501
1441     9         8085982d-4*t14*t36*t38+1.4447495018085982d-4*t13*t34*t38+
1442     :         7.22374750904299d-4*t11*t12*t34*t38-1.4447495018085982d-4
1443     ;         *t2*t32*t38+7.22374750904299d-5*t15*t23*t38+3.61187375452
1444     <         1495d-4*t15*t21*t22*t38+2.8894990036171964d-4*t15*t19*t20
1445     =         *t38+2.167124252712897d-4*t15*t17*t18*t38+1.4447495018085
1446     >         982d-4*t14*t15*t16*t38)-1.2407009817988002d+0*t25*(-1.083
1447     ?         5621263564485d-4*t10*t34*t41*t9+2.708905315891122d-4*t12*
1448     @         t32*t41*t9-8.126715947673364d-5*t34*t41*t7*t8+8.126715947
1449     1         673364d-5*t32*t41*t6*t8+1.625343189534673d-4*t10*t32*t41*
1450     2         t7-5.417810631782243d-5*t2*t34*t41*t6+2.708905315891122d-
1451     3         4*t19*t22*t36*t41+1.625343189534673d-4*t17*t20*t36*t41+8.
1452     4         126715947673364d-5*t16*t18*t36*t41+2.7089053158911214d-5*
1453     5         t14*t36*t41-2.7089053158911214d-5*t13*t34*t41-1.354452657
1454     6         945561d-4*t11*t12*t34*t41+2.7089053158911214d-5*t2*t32*t4
1455     7         1-1.3544526579455607d-5*t15*t23*t41-6.772263289727805d-5*
1456     8         t15*t21*t22*t41-5.417810631782243d-5*t15*t19*t20*t41-4.06
1457     9         3357973836682d-5*t15*t17*t18*t41-2.7089053158911214d-5*t1
1458     :         4*t15*t16*t41))*wght+Cmat3(iq,D3_RA_GAA_GAA)
1459            Cmat3(iq,D3_RA_GAA_GAB) = Cmat3(iq,D3_RA_GAA_GAB)
1460            Cmat3(iq,D3_RA_GAA_GBB) = Cmat3(iq,D3_RA_GAA_GBB)
1461            Cmat3(iq,D3_RA_GAB_GAB) = Cmat3(iq,D3_RA_GAB_GAB)
1462            Cmat3(iq,D3_RA_GAB_GBB) = Cmat3(iq,D3_RA_GAB_GBB)
1463            Cmat3(iq,D3_RA_GBB_GBB) = Cmat3(iq,D3_RA_GBB_GBB)
1464            Cmat3(iq,D3_GAA_GAA_GAA) = Cmat3(iq,D3_GAA_GAA_GAA)-9.305257
1465     1         363491002d-1*t1*(t26*t3*(-1.9939108471767675d-6*t12*t34*t
1466     2         41*t9-5.981732541530301d-7*t34*t41*t6*t8+2.99086627076515
1467     3         07d-7*t32*t41*t8-1.1963465083060604d-6*t10*t34*t41*t7+2.9
1468     4         90866270765151d-6*t12*t32*t41*t7+1.1963465083060604d-6*t1
1469     5         0*t32*t41*t6-1.9939108471767675d-7*t2*t34*t41)+t3*(2.9908
1470     6         66270765151d-6*t17*t22*t41*t49+1.1963465083060604d-6*t16*
1471     7         t20*t41*t49+2.9908662707651507d-7*t18*t41*t49-2.990866270
1472     8         765151d-6*t19*t22*t36*t41-1.7945197624590903d-6*t17*t20*t
1473     9         36*t41-8.972598812295453d-7*t16*t18*t36*t41-2.99086627076
1474     :         5151d-7*t14*t36*t41+4.9847771179419187d-8*t15*t23*t41+2.4
1475     ;         923885589709593d-7*t15*t21*t22*t41+1.9939108471767675d-7*
1476     <         t15*t19*t20*t41+1.4954331353825753d-7*t15*t17*t18*t41+9.9
1477     =         69554235883837d-8*t14*t15*t16*t41)-1.4721155704335503d-2*
1478     >         t34*t41*t42+8.126715947673364d-5*t29*t32/rhoa**8)*wght
1479            Cmat3(iq,D3_GAA_GAA_GAB) = Cmat3(iq,D3_GAA_GAA_GAB)
1480            Cmat3(iq,D3_GAA_GAA_GBB) = Cmat3(iq,D3_GAA_GAA_GBB)
1481            Cmat3(iq,D3_GAA_GAB_GAB) = Cmat3(iq,D3_GAA_GAB_GAB)
1482            Cmat3(iq,D3_GAA_GAB_GBB) = Cmat3(iq,D3_GAA_GAB_GBB)
1483            Cmat3(iq,D3_GAB_GAB_GAB) = Cmat3(iq,D3_GAB_GAB_GAB)
1484          endif ! rhoa.gt.tol_rho
1485        else  ! ipol.eq.1
1486          rhoa    = rho(iq,R_A)
1487          rhob    = rho(iq,R_B)
1488          gammaaa = rgamma(iq,G_AA)
1489          gammaab = rgamma(iq,G_AB)
1490          gammabb = rgamma(iq,G_BB)
1491          if (rhoa.gt.tol_rho.and.rhob.gt.tol_rho) then
1492            t1 = rhoa**1.3333333333333333d+0
1493            t2 = param(1)
1494            t3 = param(7)
1495            t4 = param(3)
1496            t5 = 1/rhoa**2.6666666666666666d+0
1497            t6 = 3.6802889260838756d-3*gammaaa*t5
1498            t7 = t6+1.0d+0
1499            t8 = 1.0d+0-1.0d+0/t7
1500            t9 = t8**2.0d+0
1501            t10 = param(4)
1502            t11 = t8**3.0d+0
1503            t12 = param(5)
1504            t13 = t8**4.0d+0
1505            t14 = param(6)
1506            t15 = param(2)
1507            t16 = param(9)
1508            t17 = exp(-t6)
1509            t18 = 1.0d+0-t17
1510            t19 = t18**2.0d+0
1511            t20 = param(10)
1512            t21 = t18**3.0d+0
1513            t22 = param(11)
1514            t23 = t18**4.0d+0
1515            t24 = param(12)
1516            t25 = param(8)
1517            t26 = t4*t9+t14*t8**5.0d+0+t15*t8+t3+t18*t25+t18**5.0d+0*t24
1518     1         +t22*t23+t20*t21+t2+t16*t19+t12*t13+t10*t11
1519            t27 = rhob**1.3333333333333333d+0
1520            t28 = 1/rhob**2.6666666666666666d+0
1521            t29 = 3.6802889260838756d-3*gammabb*t28
1522            t30 = t29+1.0d+0
1523            t31 = 1.0d+0-1.0d+0/t30
1524            t32 = t31**2.0d+0
1525            t33 = t31**3.0d+0
1526            t34 = t31**4.0d+0
1527            t35 = exp(-t29)
1528            t36 = 1.0d+0-t35
1529            t37 = t36**2.0d+0
1530            t38 = t36**3.0d+0
1531            t39 = t36**4.0d+0
1532            t40 = t32*t4+t22*t39+t20*t38+t16*t37+t24*t36**5.0d+0+t25*t36
1533     1         +t12*t34+t10*t33+t14*t31**5.0d+0+t15*t31+t3+t2
1534            t41 = rhoa**3.333333333333333d-1
1535            t42 = 1/t7**2
1536            t43 = 1/rhoa**3.6666666666666664d+0
1537            t44 = -2.9442311408671007d-2*gammaaa*t10*t42*t43*t9-1.962820
1538     1         760578067d-2*gammaaa*t4*t42*t43*t8-9.814103802890335d-3*g
1539     2         ammaaa*t15*t42*t43-4.9070519014451675d-2*gammaaa*t13*t14*
1540     3         t42*t43-3.925641521156134d-2*gammaaa*t11*t12*t42*t43-9.81
1541     4         4103802890335d-3*gammaaa*t17*t25*t43-4.9070519014451675d-
1542     5         2*gammaaa*t17*t23*t24*t43-3.925641521156134d-2*gammaaa*t1
1543     6         7*t21*t22*t43-2.9442311408671007d-2*gammaaa*t17*t19*t20*t
1544     7         43-1.962820760578067d-2*gammaaa*t16*t17*t18*t43
1545            t45 = rhob**3.333333333333333d-1
1546            t46 = 1/t30**2
1547            t47 = 1/rhob**3.6666666666666664d+0
1548            t48 = -1.962820760578067d-2*gammabb*t31*t4*t46*t47-4.9070519
1549     1         014451675d-2*gammabb*t14*t34*t46*t47-3.925641521156134d-2
1550     2         *gammabb*t12*t33*t46*t47-2.9442311408671007d-2*gammabb*t1
1551     3         0*t32*t46*t47-9.814103802890335d-3*gammabb*t15*t46*t47-4.
1552     4         9070519014451675d-2*gammabb*t24*t35*t39*t47-3.92564152115
1553     5         6134d-2*gammabb*t22*t35*t38*t47-2.9442311408671007d-2*gam
1554     6         mabb*t20*t35*t37*t47-1.962820760578067d-2*gammabb*t16*t35
1555     7         *t36*t47-9.814103802890335d-3*gammabb*t25*t35*t47
1556            t49 = 3.6802889260838756d-3*t15
1557            t50 = 1.1040866778251626d-2*t10*t9+7.360577852167751d-3*t4*t
1558     1         8+t49+1.8401444630419378d-2*t13*t14+1.4721155704335503d-2
1559     2         *t11*t12
1560            t51 = t49+7.360577852167751d-3*t31*t4+1.8401444630419378d-2*
1561     1         t14*t34+1.4721155704335503d-2*t12*t33+1.1040866778251626d
1562     2         -2*t10*t32
1563            t52 = 1/rhoa**6.666666666666666d-1
1564            t53 = gammaaa**2
1565            t54 = 1/t7**4
1566            t55 = 1/rhoa**7.333333333333333d+0
1567            t56 = 1/t7**3
1568            t57 = 1/rhoa**4.666666666666667d+0
1569            t58 = exp(-7.360577852167751d-3*gammaaa*t5)
1570            t59 = 1.0795514183179368d-1*gammaaa*t10*t42*t57*t9-5.7789980
1571     1         07234393d-4*t10*t53*t55*t56*t9+1.1557996014468785d-3*t12*
1572     2         t53*t54*t55*t9+7.197009455452912d-2*gammaaa*t4*t42*t57*t8
1573     3         -3.8526653381562614d-4*t4*t53*t55*t56*t8+5.77899800723439
1574     4         3d-4*t10*t53*t54*t55*t8+1.9263326690781307d-3*t21*t24*t53
1575     5         *t55*t58+1.1557996014468785d-3*t19*t22*t53*t55*t58+5.7789
1576     6         98007234393d-4*t18*t20*t53*t55*t58+1.9263326690781307d-4*
1577     7         t16*t53*t55*t58+3.598504727726456d-2*gammaaa*t15*t42*t57+
1578     8         1.799252363863228d-1*gammaaa*t13*t14*t42*t57+1.4394018910
1579     9         905823d-1*gammaaa*t11*t12*t42*t57+3.598504727726456d-2*ga
1580     :         mmaaa*t17*t25*t57+1.799252363863228d-1*gammaaa*t17*t23*t2
1581     ;         4*t57+1.4394018910905823d-1*gammaaa*t17*t21*t22*t57+1.079
1582     <         5514183179368d-1*gammaaa*t17*t19*t20*t57+7.19700945545291
1583     =         2d-2*gammaaa*t16*t17*t18*t57-1.9263326690781307d-4*t15*t5
1584     >         3*t55*t56-9.631663345390654d-4*t13*t14*t53*t55*t56-7.7053
1585     ?         30676312523d-4*t11*t12*t53*t55*t56+1.9263326690781307d-4*
1586     @         t4*t53*t54*t55+1.9263326690781307d-3*t11*t14*t53*t54*t55-
1587     1         9.631663345390653d-5*t17*t25*t53*t55-4.815831672695327d-4
1588     2         *t17*t23*t24*t53*t55-3.8526653381562614d-4*t17*t21*t22*t5
1589     3         3*t55-2.8894990036171964d-4*t17*t19*t20*t53*t55-1.9263326
1590     4         690781307d-4*t16*t17*t18*t53*t55
1591            t60 = 1/rhob**6.666666666666666d-1
1592            t61 = gammabb**2
1593            t62 = 1/t30**4
1594            t63 = 1/rhob**7.333333333333333d+0
1595            t64 = 1/t30**3
1596            t65 = 1/rhob**4.666666666666667d+0
1597            t66 = exp(-7.360577852167751d-3*gammabb*t28)
1598            t67 = 1.9263326690781307d-3*t24*t38*t61*t63*t66+1.1557996014
1599     1         468785d-3*t22*t37*t61*t63*t66+5.778998007234393d-4*t20*t3
1600     2         6*t61*t63*t66+1.9263326690781307d-4*t16*t61*t63*t66+7.197
1601     3         009455452912d-2*gammabb*t31*t4*t46*t65+1.799252363863228d
1602     4         -1*gammabb*t14*t34*t46*t65+1.4394018910905823d-1*gammabb*
1603     5         t12*t33*t46*t65+1.0795514183179368d-1*gammabb*t10*t32*t46
1604     6         *t65+3.598504727726456d-2*gammabb*t15*t46*t65+1.799252363
1605     7         863228d-1*gammabb*t24*t35*t39*t65+1.4394018910905823d-1*g
1606     8         ammabb*t22*t35*t38*t65+1.0795514183179368d-1*gammabb*t20*
1607     9         t35*t37*t65+7.197009455452912d-2*gammabb*t16*t35*t36*t65+
1608     :         3.598504727726456d-2*gammabb*t25*t35*t65-3.85266533815626
1609     ;         14d-4*t31*t4*t61*t63*t64-9.631663345390654d-4*t14*t34*t61
1610     <         *t63*t64-7.705330676312523d-4*t12*t33*t61*t63*t64-5.77899
1611     =         8007234393d-4*t10*t32*t61*t63*t64-1.9263326690781307d-4*t
1612     >         15*t61*t63*t64+1.9263326690781307d-4*t4*t61*t62*t63+1.926
1613     ?         3326690781307d-3*t14*t33*t61*t62*t63+1.1557996014468785d-
1614     @         3*t12*t32*t61*t62*t63+5.778998007234393d-4*t10*t31*t61*t6
1615     1         2*t63-4.815831672695327d-4*t24*t35*t39*t61*t63-3.85266533
1616     2         81562614d-4*t22*t35*t38*t61*t63-2.8894990036171964d-4*t20
1617     3         *t35*t37*t61*t63-1.9263326690781307d-4*t16*t35*t36*t61*t6
1618     4         3-9.631663345390653d-5*t25*t35*t61*t63
1619            t68 = 1/rhoa**6.333333333333333d+0
1620            t69 = 2.167124252712897d-4*gammaaa*t10*t56*t68*t9-4.33424850
1621     1         5425794d-4*gammaaa*t12*t54*t68*t9-2.9442311408671007d-2*t
1622     2         10*t42*t43*t9+1.4447495018085982d-4*gammaaa*t4*t56*t68*t8
1623     3         -2.167124252712897d-4*gammaaa*t10*t54*t68*t8-1.9628207605
1624     4         78067d-2*t4*t42*t43*t8-7.22374750904299d-4*gammaaa*t21*t2
1625     5         4*t58*t68-4.334248505425794d-4*gammaaa*t19*t22*t58*t68-2.
1626     6         167124252712897d-4*gammaaa*t18*t20*t58*t68-7.223747509042
1627     7         99d-5*gammaaa*t16*t58*t68+7.22374750904299d-5*gammaaa*t15
1628     8         *t56*t68+3.611873754521495d-4*gammaaa*t13*t14*t56*t68+2.8
1629     9         894990036171964d-4*gammaaa*t11*t12*t56*t68-7.223747509042
1630     :         99d-5*gammaaa*t4*t54*t68-7.22374750904299d-4*gammaaa*t11*
1631     ;         t14*t54*t68+3.611873754521495d-5*gammaaa*t17*t25*t68+1.80
1632     <         59368772607476d-4*gammaaa*t17*t23*t24*t68+1.4447495018085
1633     =         982d-4*gammaaa*t17*t21*t22*t68+1.0835621263564485d-4*gamm
1634     >         aaa*t17*t19*t20*t68+7.22374750904299d-5*gammaaa*t16*t17*t
1635     ?         18*t68-9.814103802890335d-3*t15*t42*t43-4.907051901445167
1636     @         5d-2*t13*t14*t42*t43-3.925641521156134d-2*t11*t12*t42*t43
1637     1         -9.814103802890335d-3*t17*t25*t43-4.9070519014451675d-2*t
1638     2         17*t23*t24*t43-3.925641521156134d-2*t17*t21*t22*t43-2.944
1639     3         2311408671007d-2*t17*t19*t20*t43-1.962820760578067d-2*t16
1640     4         *t17*t18*t43
1641            t70 = 1.1040866778251626d-2*t10*t42*t5*t9+7.360577852167751d
1642     1         -3*t4*t42*t5*t8+3.6802889260838756d-3*t15*t42*t5+1.840144
1643     2         4630419378d-2*t13*t14*t42*t5+1.4721155704335503d-2*t11*t1
1644     3         2*t42*t5+3.6802889260838756d-3*t17*t25*t5+1.8401444630419
1645     4         378d-2*t17*t23*t24*t5+1.4721155704335503d-2*t17*t21*t22*t
1646     5         5+1.1040866778251626d-2*t17*t19*t20*t5+7.360577852167751d
1647     6         -3*t16*t17*t18*t5
1648            t71 = 1/rhob**6.333333333333333d+0
1649            t72 = -7.22374750904299d-4*gammabb*t24*t38*t66*t71-4.3342485
1650     1         05425794d-4*gammabb*t22*t37*t66*t71-2.167124252712897d-4*
1651     2         gammabb*t20*t36*t66*t71-7.22374750904299d-5*gammabb*t16*t
1652     3         66*t71+1.4447495018085982d-4*gammabb*t31*t4*t64*t71+3.611
1653     4         873754521495d-4*gammabb*t14*t34*t64*t71+2.889499003617196
1654     5         4d-4*gammabb*t12*t33*t64*t71+2.167124252712897d-4*gammabb
1655     6         *t10*t32*t64*t71+7.22374750904299d-5*gammabb*t15*t64*t71-
1656     7         7.22374750904299d-5*gammabb*t4*t62*t71-7.22374750904299d-
1657     8         4*gammabb*t14*t33*t62*t71-4.334248505425794d-4*gammabb*t1
1658     9         2*t32*t62*t71-2.167124252712897d-4*gammabb*t10*t31*t62*t7
1659     :         1+1.8059368772607476d-4*gammabb*t24*t35*t39*t71+1.4447495
1660     ;         018085982d-4*gammabb*t22*t35*t38*t71+1.0835621263564485d-
1661     <         4*gammabb*t20*t35*t37*t71+7.22374750904299d-5*gammabb*t16
1662     =         *t35*t36*t71+3.611873754521495d-5*gammabb*t25*t35*t71-1.9
1663     >         62820760578067d-2*t31*t4*t46*t47-4.9070519014451675d-2*t1
1664     ?         4*t34*t46*t47-3.925641521156134d-2*t12*t33*t46*t47-2.9442
1665     @         311408671007d-2*t10*t32*t46*t47-9.814103802890335d-3*t15*
1666     1         t46*t47-4.9070519014451675d-2*t24*t35*t39*t47-3.925641521
1667     2         156134d-2*t22*t35*t38*t47-2.9442311408671007d-2*t20*t35*t
1668     3         37*t47-1.962820760578067d-2*t16*t35*t36*t47-9.81410380289
1669     4         0335d-3*t25*t35*t47
1670            t73 = 7.360577852167751d-3*t28*t31*t4*t46+1.8401444630419378
1671     1         d-2*t14*t28*t34*t46+1.4721155704335503d-2*t12*t28*t33*t46
1672     2         +1.1040866778251626d-2*t10*t28*t32*t46+3.6802889260838756
1673     3         d-3*t15*t28*t46+1.8401444630419378d-2*t24*t28*t35*t39+1.4
1674     4         721155704335503d-2*t22*t28*t35*t38+1.1040866778251626d-2*
1675     5         t20*t28*t35*t37+7.360577852167751d-3*t16*t28*t35*t36+3.68
1676     6         02889260838756d-3*t25*t28*t35
1677            t74 = 1/rhoa**5.333333333333333d+0
1678            t75 = 1.625343189534673d-4*t12*t42*t5*t9+8.126715947673364d-
1679     1         5*t10*t42*t5*t8+2.7089053158911214d-5*t4*t42*t5+2.7089053
1680     2         15891122d-4*t11*t14*t42*t5
1681            t76 = 1/rhob**5.333333333333333d+0
1682            t77 = 2.7089053158911214d-5*t28*t4*t46+2.708905315891122d-4*
1683     1         t14*t28*t33*t46+1.625343189534673d-4*t12*t28*t32*t46+8.12
1684     2         6715947673364d-5*t10*t28*t31*t46
1685            t78 = gammaaa**3
1686            t79 = 1/t7**6
1687            t80 = 1/rhoa**11
1688            t81 = 1/t7**5
1689            t82 = 1/rhoa**8.333333333333334d+0
1690            t83 = 1/rhoa**5.666666666666667d+0
1691            t84 = exp(-1.1040866778251626d-2*gammaaa*t5)
1692            t85 = gammabb**3
1693            t86 = 1/t30**6
1694            t87 = 1/rhob**11
1695            t88 = 1/t30**5
1696            t89 = 1/rhob**8.333333333333334d+0
1697            t90 = 1/rhob**5.666666666666667d+0
1698            t91 = exp(-1.1040866778251626d-2*gammabb*t28)
1699            t92 = 1/rhoa**10
1700            t93 = 1/rhob**10
1701            t94 = 1/rhoa**9
1702            t95 = 1/rhob**9
1703            fnc(iq) = (-9.305257363491002d-1*t27*t40-9.305257363491002d-
1704     1         1*t1*t26)*wght+fnc(iq)
1705            Amat(iq,D1_RA) = (-9.305257363491002d-1*t1*t44-1.24070098179
1706     1         88002d+0*t26*t41)*wght+Amat(iq,D1_RA)
1707            Amat(iq,D1_RB) = (-9.305257363491002d-1*t27*t48-1.2407009817
1708     1         988002d+0*t40*t45)*wght+Amat(iq,D1_RB)
1709            Cmat(iq,D1_GAA) = Cmat(iq,D1_GAA)-9.305257363491002d-1*t1*(t
1710     1         42*t5*t50+(3.6802889260838756d-3*t17*t25+1.84014446304193
1711     2         78d-2*t17*t23*t24+1.4721155704335503d-2*t17*t21*t22+1.104
1712     3         0866778251626d-2*t17*t19*t20+7.360577852167751d-3*t16*t17
1713     4         *t18)*t5)*wght
1714            Cmat(iq,D1_GAB) = Cmat(iq,D1_GAB)
1715            Cmat(iq,D1_GBB) = Cmat(iq,D1_GBB)-9.305257363491002d-1*t27*(
1716     1         t28*t46*t51+t28*(1.8401444630419378d-2*t24*t35*t39+1.4721
1717     2         155704335503d-2*t22*t35*t38+1.1040866778251626d-2*t20*t35
1718     3         *t37+7.360577852167751d-3*t16*t35*t36+3.6802889260838756d
1719     4         -3*t25*t35))*wght
1720            Amat2(iq,D2_RA_RA) = (-9.305257363491002d-1*t1*t59-4.1356699
1721     1         39329334d-1*t26*t52-2.4814019635976003d+0*t41*t44)*wght+A
1722     2         mat2(iq,D2_RA_RA)
1723            Amat2(iq,D2_RA_RB) = Amat2(iq,D2_RA_RB)
1724            Amat2(iq,D2_RB_RB) = (-9.305257363491002d-1*t27*t67-4.135669
1725     1         939329334d-1*t40*t60-2.4814019635976003d+0*t45*t48)*wght+
1726     2         Amat2(iq,D2_RB_RB)
1727            Cmat2(iq,D2_RA_GAA) = (-1.2407009817988002d+0*t41*t70-9.3052
1728     1         57363491002d-1*t1*t69)*wght+Cmat2(iq,D2_RA_GAA)
1729            Cmat2(iq,D2_RA_GAB) = Cmat2(iq,D2_RA_GAB)
1730            Cmat2(iq,D2_RA_GBB) = Cmat2(iq,D2_RA_GBB)
1731            Cmat2(iq,D2_RB_GAA) = Cmat2(iq,D2_RB_GAA)
1732            Cmat2(iq,D2_RB_GAB) = Cmat2(iq,D2_RB_GAB)
1733            Cmat2(iq,D2_RB_GBB) = (-1.2407009817988002d+0*t45*t73-9.3052
1734     1         57363491002d-1*t27*t72)*wght+Cmat2(iq,D2_RB_GBB)
1735            Cmat2(iq,D2_GAA_GAA) = Cmat2(iq,D2_GAA_GAA)-9.30525736349100
1736     1         2d-1*t1*(t42*t5*t75-7.360577852167751d-3*t50*t56*t74+t5*(
1737     2         2.708905315891122d-4*t21*t24*t5*t58+1.625343189534673d-4*
1738     3         t19*t22*t5*t58+8.126715947673364d-5*t18*t20*t5*t58+2.7089
1739     4         053158911214d-5*t16*t5*t58-1.3544526579455607d-5*t17*t25*
1740     5         t5-6.772263289727805d-5*t17*t23*t24*t5-5.417810631782243d
1741     6         -5*t17*t21*t22*t5-4.063357973836682d-5*t17*t19*t20*t5-2.7
1742     7         089053158911214d-5*t16*t17*t18*t5))*wght
1743            Cmat2(iq,D2_GAA_GAB) = Cmat2(iq,D2_GAA_GAB)
1744            Cmat2(iq,D2_GAA_GBB) = Cmat2(iq,D2_GAA_GBB)
1745            Cmat2(iq,D2_GAB_GAB) = Cmat2(iq,D2_GAB_GAB)
1746            Cmat2(iq,D2_GAB_GBB) = Cmat2(iq,D2_GAB_GBB)
1747            Cmat2(iq,D2_GBB_GBB) = Cmat2(iq,D2_GBB_GBB)-9.30525736349100
1748     1         2d-1*t27*(t28*t46*t77-7.360577852167751d-3*t51*t64*t76+t2
1749     2         8*(2.708905315891122d-4*t24*t28*t38*t66+1.625343189534673
1750     3         d-4*t22*t28*t37*t66+8.126715947673364d-5*t20*t28*t36*t66+
1751     4         2.7089053158911214d-5*t16*t28*t66-6.772263289727805d-5*t2
1752     5         4*t28*t35*t39-5.417810631782243d-5*t22*t28*t35*t38-4.0633
1753     6         57973836682d-5*t20*t28*t35*t37-2.7089053158911214d-5*t16*
1754     7         t28*t35*t36-1.3544526579455607d-5*t25*t28*t35))*wght
1755            Amat3(iq,D3_RA_RA_RA) = (-9.305257363491002d-1*t1*(-5.037906
1756     1         618817038d-1*gammaaa*t10*t42*t83*t9+6.356897807957831d-3*
1757     2         t10*t53*t56*t82*t9-1.2713795615915663d-2*t12*t53*t54*t82*
1758     3         t9+6.805882358363365d-5*t12*t78*t80*t81*t9-5.671568631969
1759     4         472d-5*t14*t78*t79*t80*t9-1.7014705895908416d-5*t10*t54*t
1760     5         78*t80*t9-5.671568631969472d-5*t19*t24*t78*t80*t84-2.2686
1761     6         274527877887d-5*t18*t22*t78*t80*t84-5.671568631969471d-6*
1762     7         t20*t78*t80*t84-3.3586044125446923d-1*gammaaa*t4*t42*t8*t
1763     8         83-1.6793022062723462d-1*gammaaa*t15*t42*t83-8.3965110313
1764     9         61731d-1*gammaaa*t13*t14*t42*t83-6.717208825089385d-1*gam
1765     :         maaa*t11*t12*t42*t83-1.6793022062723462d-1*gammaaa*t17*t2
1766     ;         5*t83-8.396511031361731d-1*gammaaa*t17*t23*t24*t83-6.7172
1767     <         08825089385d-1*gammaaa*t17*t21*t22*t83-5.037906618817038d
1768     =         -1*gammaaa*t17*t19*t20*t83-3.3586044125446923d-1*gammaaa*
1769     >         t16*t17*t18*t83+4.237931871971887d-3*t4*t53*t56*t8*t82-6.
1770     ?         356897807957831d-3*t10*t53*t54*t8*t82-2.1189659359859436d
1771     @         -2*t21*t24*t53*t58*t82-1.2713795615915663d-2*t19*t22*t53*
1772     1         t58*t82-6.356897807957831d-3*t18*t20*t53*t58*t82-2.118965
1773     2         9359859436d-3*t16*t53*t58*t82+2.1189659359859436d-3*t15*t
1774     3         53*t56*t82+1.0594829679929718d-2*t13*t14*t53*t56*t82+8.47
1775     4         5863743943775d-3*t11*t12*t53*t56*t82-2.1189659359859436d-
1776     5         3*t4*t53*t54*t82-2.1189659359859436d-2*t11*t14*t53*t54*t8
1777     6         2+1.0594829679929718d-3*t17*t25*t53*t82+5.297414839964859
1778     7         d-3*t17*t23*t24*t53*t82+4.237931871971887d-3*t17*t21*t22*
1779     8         t53*t82+3.1784489039789154d-3*t17*t19*t20*t53*t82+2.11896
1780     9         59359859436d-3*t16*t17*t18*t53*t82+3.4029411791816827d-5*
1781     :         t10*t78*t8*t80*t81+1.1343137263938943d-5*t4*t78*t80*t81+1
1782     ;         .1343137263938945d-4*t11*t14*t78*t80*t81-2.26862745278778
1783     <         87d-5*t12*t78*t79*t8*t80-1.1343137263938943d-5*t4*t54*t78
1784     =         *t8*t80-5.671568631969471d-6*t10*t78*t79*t80+5.6715686319
1785     >         69472d-5*t21*t24*t58*t78*t80+3.4029411791816827d-5*t19*t2
1786     ?         2*t58*t78*t80+1.7014705895908414d-5*t18*t20*t58*t78*t80+5
1787     @         .671568631969471d-6*t16*t58*t78*t80-5.671568631969471d-6*
1788     1         t15*t54*t78*t80-2.835784315984736d-5*t13*t14*t54*t78*t80-
1789     2         2.2686274527877887d-5*t11*t12*t54*t78*t80-9.4526143866157
1790     3         86d-7*t17*t25*t78*t80-4.726307193307893d-6*t17*t23*t24*t7
1791     4         8*t80-3.7810457546463144d-6*t17*t21*t22*t78*t80-2.8357843
1792     5         159847357d-6*t17*t19*t20*t78*t80-1.8905228773231572d-6*t1
1793     6         6*t17*t18*t78*t80)-3.7221029453964005d+0*t41*t59-1.240700
1794     7         9817988002d+0*t44*t52+2.7571132928862224d-1*t26/rhoa**1.6
1795     8         666666666666669d+0)*wght+Amat3(iq,D3_RA_RA_RA)
1796            Amat3(iq,D3_RA_RA_RB) = Amat3(iq,D3_RA_RA_RB)
1797            Amat3(iq,D3_RA_RB_RB) = Amat3(iq,D3_RA_RB_RB)
1798            Amat3(iq,D3_RB_RB_RB) = (-9.305257363491002d-1*t27*(-5.67156
1799     1         8631969472d-5*t24*t37*t85*t87*t91-2.2686274527877887d-5*t
1800     2         22*t36*t85*t87*t91-5.671568631969471d-6*t20*t85*t87*t91-3
1801     3         .3586044125446923d-1*gammabb*t31*t4*t46*t90-8.39651103136
1802     4         1731d-1*gammabb*t14*t34*t46*t90-6.717208825089385d-1*gamm
1803     5         abb*t12*t33*t46*t90-5.037906618817038d-1*gammabb*t10*t32*
1804     6         t46*t90-1.6793022062723462d-1*gammabb*t15*t46*t90-8.39651
1805     7         1031361731d-1*gammabb*t24*t35*t39*t90-6.717208825089385d-
1806     8         1*gammabb*t22*t35*t38*t90-5.037906618817038d-1*gammabb*t2
1807     9         0*t35*t37*t90-3.3586044125446923d-1*gammabb*t16*t35*t36*t
1808     :         90-1.6793022062723462d-1*gammabb*t25*t35*t90-2.1189659359
1809     ;         859436d-2*t24*t38*t61*t66*t89-1.2713795615915663d-2*t22*t
1810     <         37*t61*t66*t89-6.356897807957831d-3*t20*t36*t61*t66*t89-2
1811     =         .1189659359859436d-3*t16*t61*t66*t89+4.237931871971887d-3
1812     >         *t31*t4*t61*t64*t89+1.0594829679929718d-2*t14*t34*t61*t64
1813     ?         *t89+8.475863743943775d-3*t12*t33*t61*t64*t89+6.356897807
1814     @         957831d-3*t10*t32*t61*t64*t89+2.1189659359859436d-3*t15*t
1815     1         61*t64*t89-2.1189659359859436d-3*t4*t61*t62*t89-2.1189659
1816     2         359859436d-2*t14*t33*t61*t62*t89-1.2713795615915663d-2*t1
1817     3         2*t32*t61*t62*t89-6.356897807957831d-3*t10*t31*t61*t62*t8
1818     4         9+5.297414839964859d-3*t24*t35*t39*t61*t89+4.237931871971
1819     5         887d-3*t22*t35*t38*t61*t89+3.1784489039789154d-3*t20*t35*
1820     6         t37*t61*t89+2.1189659359859436d-3*t16*t35*t36*t61*t89+1.0
1821     7         594829679929718d-3*t25*t35*t61*t89+1.1343137263938943d-5*
1822     8         t4*t85*t87*t88+1.1343137263938945d-4*t14*t33*t85*t87*t88+
1823     9         6.805882358363365d-5*t12*t32*t85*t87*t88+3.40294117918168
1824     :         27d-5*t10*t31*t85*t87*t88-5.671568631969472d-5*t14*t32*t8
1825     ;         5*t86*t87-2.2686274527877887d-5*t12*t31*t85*t86*t87-5.671
1826     <         568631969471d-6*t10*t85*t86*t87+5.671568631969472d-5*t24*
1827     =         t38*t66*t85*t87+3.4029411791816827d-5*t22*t37*t66*t85*t87
1828     >         +1.7014705895908414d-5*t20*t36*t66*t85*t87+5.671568631969
1829     ?         471d-6*t16*t66*t85*t87-1.1343137263938943d-5*t31*t4*t62*t
1830     @         85*t87-2.835784315984736d-5*t14*t34*t62*t85*t87-2.2686274
1831     1         527877887d-5*t12*t33*t62*t85*t87-1.7014705895908416d-5*t1
1832     2         0*t32*t62*t85*t87-5.671568631969471d-6*t15*t62*t85*t87-4.
1833     3         726307193307893d-6*t24*t35*t39*t85*t87-3.7810457546463144
1834     4         d-6*t22*t35*t38*t85*t87-2.8357843159847357d-6*t20*t35*t37
1835     5         *t85*t87-1.8905228773231572d-6*t16*t35*t36*t85*t87-9.4526
1836     6         14386615786d-7*t25*t35*t85*t87)-3.7221029453964005d+0*t45
1837     7         *t67-1.2407009817988002d+0*t48*t60+2.7571132928862224d-1*
1838     8         t40/rhob**1.6666666666666669d+0)*wght+Amat3(iq,D3_RB_RB_R
1839     9         B)
1840            Cmat3(iq,D3_RA_RA_GAA) = (-9.305257363491002d-1*t1*(-2.55220
1841     1         5884386262d-5*t12*t53*t81*t9*t92+2.1268382369885516d-5*t1
1842     2         4*t53*t79*t9*t92+6.380514710965656d-6*t10*t53*t54*t9*t92+
1843     3         2.1268382369885516d-5*t19*t24*t53*t84*t92+8.5073529479542
1844     4         09d-6*t18*t22*t53*t84*t92+2.1268382369885522d-6*t20*t53*t
1845     5         84*t92-1.2761029421931314d-5*t10*t53*t8*t81*t92-4.2536764
1846     6         73977103d-6*t4*t53*t81*t92-4.2536764739771044d-5*t11*t14*
1847     7         t53*t81*t92+8.507352947954209d-6*t12*t53*t79*t8*t92+4.253
1848     8         676473977103d-6*t4*t53*t54*t8*t92+2.1268382369885522d-6*t
1849     9         10*t53*t79*t92-2.1268382369885522d-5*t21*t24*t53*t58*t92-
1850     :         1.276102942193131d-5*t19*t22*t53*t58*t92-6.38051471096565
1851     ;         7d-6*t18*t20*t53*t58*t92-2.1268382369885516d-6*t16*t53*t5
1852     <         8*t92+2.1268382369885516d-6*t15*t53*t54*t92+1.06341911849
1853     =         42758d-5*t13*t14*t53*t54*t92+8.507352947954206d-6*t11*t12
1854     >         *t53*t54*t92+3.54473039498092d-7*t17*t25*t53*t92+1.772365
1855     ?         19749046d-6*t17*t23*t24*t53*t92+1.4178921579923678d-6*t17
1856     @         *t21*t22*t53*t92+1.0634191184942761d-6*t17*t19*t20*t53*t9
1857     1         2+7.08946078996184d-7*t16*t17*t18*t53*t92+1.0795514183179
1858     2         368d-1*t10*t42*t57*t9-1.9504118274416074d-3*gammaaa*t10*t
1859     3         55*t56*t9+3.900823654883215d-3*gammaaa*t12*t54*t55*t9+7.1
1860     4         97009455452912d-2*t4*t42*t57*t8-1.300274551627738d-3*gamm
1861     5         aaa*t4*t55*t56*t8+1.9504118274416074d-3*gammaaa*t10*t54*t
1862     6         55*t8+6.50137275813869d-3*gammaaa*t21*t24*t55*t58+3.90082
1863     7         3654883215d-3*gammaaa*t19*t22*t55*t58+1.9504118274416074d
1864     8         -3*gammaaa*t18*t20*t55*t58+6.501372758138692d-4*gammaaa*t
1865     9         16*t55*t58+3.598504727726456d-2*t15*t42*t57+1.79925236386
1866     :         3228d-1*t13*t14*t42*t57+1.4394018910905823d-1*t11*t12*t42
1867     ;         *t57+3.598504727726456d-2*t17*t25*t57+1.799252363863228d-
1868     <         1*t17*t23*t24*t57+1.4394018910905823d-1*t17*t21*t22*t57+1
1869     =         .0795514183179368d-1*t17*t19*t20*t57+7.197009455452912d-2
1870     >         *t16*t17*t18*t57-6.501372758138692d-4*gammaaa*t15*t55*t56
1871     ?         -3.250686379069345d-3*gammaaa*t13*t14*t55*t56-2.600549103
1872     @         255476d-3*gammaaa*t11*t12*t55*t56+6.501372758138692d-4*ga
1873     1         mmaaa*t4*t54*t55+6.50137275813869d-3*gammaaa*t11*t14*t54*
1874     2         t55-3.250686379069346d-4*gammaaa*t17*t25*t55-1.6253431895
1875     3         346726d-3*gammaaa*t17*t23*t24*t55-1.300274551627738d-3*ga
1876     4         mmaaa*t17*t21*t22*t55-9.752059137208037d-4*gammaaa*t17*t1
1877     5         9*t20*t55-6.501372758138692d-4*gammaaa*t16*t17*t18*t55)-4
1878     6         .135669939329334d-1*t52*t70-2.4814019635976003d+0*t41*t69
1879     7         )*wght+Cmat3(iq,D3_RA_RA_GAA)
1880            Cmat3(iq,D3_RA_RA_GAB) = Cmat3(iq,D3_RA_RA_GAB)
1881            Cmat3(iq,D3_RA_RA_GBB) = Cmat3(iq,D3_RA_RA_GBB)
1882            Cmat3(iq,D3_RA_RB_GAA) = Cmat3(iq,D3_RA_RB_GAA)
1883            Cmat3(iq,D3_RA_RB_GAB) = Cmat3(iq,D3_RA_RB_GAB)
1884            Cmat3(iq,D3_RA_RB_GBB) = Cmat3(iq,D3_RA_RB_GBB)
1885            Cmat3(iq,D3_RB_RB_GAA) = Cmat3(iq,D3_RB_RB_GAA)
1886            Cmat3(iq,D3_RB_RB_GAB) = Cmat3(iq,D3_RB_RB_GAB)
1887            Cmat3(iq,D3_RB_RB_GBB) = (-9.305257363491002d-1*t27*(2.12683
1888     1         82369885516d-5*t24*t37*t61*t91*t93+8.507352947954209d-6*t
1889     2         22*t36*t61*t91*t93+2.1268382369885522d-6*t20*t61*t91*t93-
1890     3         4.253676473977103d-6*t4*t61*t88*t93-4.2536764739771044d-5
1891     4         *t14*t33*t61*t88*t93-2.552205884386262d-5*t12*t32*t61*t88
1892     5         *t93-1.2761029421931314d-5*t10*t31*t61*t88*t93+2.12683823
1893     6         69885516d-5*t14*t32*t61*t86*t93+8.507352947954209d-6*t12*
1894     7         t31*t61*t86*t93+2.1268382369885522d-6*t10*t61*t86*t93-2.1
1895     8         268382369885522d-5*t24*t38*t61*t66*t93-1.276102942193131d
1896     9         -5*t22*t37*t61*t66*t93-6.380514710965657d-6*t20*t36*t61*t
1897     :         66*t93-2.1268382369885516d-6*t16*t61*t66*t93+4.2536764739
1898     ;         77103d-6*t31*t4*t61*t62*t93+1.0634191184942758d-5*t14*t34
1899     <         *t61*t62*t93+8.507352947954206d-6*t12*t33*t61*t62*t93+6.3
1900     =         80514710965656d-6*t10*t32*t61*t62*t93+2.1268382369885516d
1901     >         -6*t15*t61*t62*t93+1.77236519749046d-6*t24*t35*t39*t61*t9
1902     ?         3+1.4178921579923678d-6*t22*t35*t38*t61*t93+1.06341911849
1903     @         42761d-6*t20*t35*t37*t61*t93+7.08946078996184d-7*t16*t35*
1904     1         t36*t61*t93+3.54473039498092d-7*t25*t35*t61*t93+6.5013727
1905     2         5813869d-3*gammabb*t24*t38*t63*t66+3.900823654883215d-3*g
1906     3         ammabb*t22*t37*t63*t66+1.9504118274416074d-3*gammabb*t20*
1907     4         t36*t63*t66+6.501372758138692d-4*gammabb*t16*t63*t66+7.19
1908     5         7009455452912d-2*t31*t4*t46*t65+1.799252363863228d-1*t14*
1909     6         t34*t46*t65+1.4394018910905823d-1*t12*t33*t46*t65+1.07955
1910     7         14183179368d-1*t10*t32*t46*t65+3.598504727726456d-2*t15*t
1911     8         46*t65+1.799252363863228d-1*t24*t35*t39*t65+1.43940189109
1912     9         05823d-1*t22*t35*t38*t65+1.0795514183179368d-1*t20*t35*t3
1913     :         7*t65+7.197009455452912d-2*t16*t35*t36*t65+3.598504727726
1914     ;         456d-2*t25*t35*t65-1.300274551627738d-3*gammabb*t31*t4*t6
1915     <         3*t64-3.250686379069345d-3*gammabb*t14*t34*t63*t64-2.6005
1916     =         49103255476d-3*gammabb*t12*t33*t63*t64-1.9504118274416074
1917     >         d-3*gammabb*t10*t32*t63*t64-6.501372758138692d-4*gammabb*
1918     ?         t15*t63*t64+6.501372758138692d-4*gammabb*t4*t62*t63+6.501
1919     @         37275813869d-3*gammabb*t14*t33*t62*t63+3.900823654883215d
1920     1         -3*gammabb*t12*t32*t62*t63+1.9504118274416074d-3*gammabb*
1921     2         t10*t31*t62*t63-1.6253431895346726d-3*gammabb*t24*t35*t39
1922     3         *t63-1.300274551627738d-3*gammabb*t22*t35*t38*t63-9.75205
1923     4         9137208037d-4*gammabb*t20*t35*t37*t63-6.501372758138692d-
1924     5         4*gammabb*t16*t35*t36*t63-3.250686379069346d-4*gammabb*t2
1925     6         5*t35*t63)-4.135669939329334d-1*t60*t73-2.481401963597600
1926     7         3d+0*t45*t72)*wght+Cmat3(iq,D3_RB_RB_GBB)
1927            Cmat3(iq,D3_RA_GAA_GAA) = (-9.305257363491002d-1*t1*(9.57077
1928     1         2066448484d-6*gammaaa*t12*t81*t9*t94-7.97564338870707d-6*
1929     2         gammaaa*t14*t79*t9*t94-2.3926930166121207d-6*gammaaa*t10*
1930     3         t54*t9*t94-7.97564338870707d-6*gammaaa*t19*t24*t84*t94-3.
1931     4         190257355482828d-6*gammaaa*t18*t22*t84*t94-7.975643388707
1932     5         07d-7*gammaaa*t20*t84*t94+4.785386033224242d-6*gammaaa*t1
1933     6         0*t8*t81*t94+1.5951286777414142d-6*gammaaa*t4*t81*t94+1.5
1934     7         95128677741414d-5*gammaaa*t11*t14*t81*t94-3.1902573554828
1935     8         28d-6*gammaaa*t12*t79*t8*t94-1.595128677741414d-6*gammaaa
1936     9         *t4*t54*t8*t94-7.97564338870707d-7*gammaaa*t10*t79*t94+7.
1937     :         97564338870707d-6*gammaaa*t21*t24*t58*t94+4.7853860332242
1938     ;         42d-6*gammaaa*t19*t22*t58*t94+2.392693016612121d-6*gammaa
1939     <         a*t18*t20*t58*t94+7.975643388707071d-7*gammaaa*t16*t58*t9
1940     =         4-7.97564338870707d-7*gammaaa*t15*t54*t94-3.9878216943535
1941     >         35d-6*gammaaa*t13*t14*t54*t94-3.190257355482828d-6*gammaa
1942     ?         a*t11*t12*t54*t94-1.329273898117845d-7*gammaaa*t17*t25*t9
1943     @         4-6.646369490589226d-7*gammaaa*t17*t23*t24*t94-5.31709559
1944     1         247138d-7*gammaaa*t17*t21*t22*t94-3.987821694353535d-7*ga
1945     2         mmaaa*t17*t19*t20*t94-2.65854779623569d-7*gammaaa*t16*t17
1946     3         *t18*t94+4.334248505425794d-4*t10*t56*t68*t9-8.6684970108
1947     4         51588d-4*t12*t54*t68*t9+2.8894990036171964d-4*t4*t56*t68*
1948     5         t8-4.334248505425794d-4*t10*t54*t68*t8-1.4447495018085982
1949     6         d-3*t21*t24*t58*t68-8.668497010851588d-4*t19*t22*t58*t68-
1950     7         4.334248505425794d-4*t18*t20*t58*t68-1.4447495018085982d-
1951     8         4*t16*t58*t68+1.4447495018085982d-4*t15*t56*t68+7.2237475
1952     9         0904299d-4*t13*t14*t56*t68+5.778998007234393d-4*t11*t12*t
1953     :         56*t68-1.4447495018085982d-4*t4*t54*t68-1.444749501808598
1954     ;         2d-3*t11*t14*t54*t68+7.22374750904299d-5*t17*t25*t68+3.61
1955     <         1873754521495d-4*t17*t23*t24*t68+2.8894990036171964d-4*t1
1956     =         7*t21*t22*t68+2.167124252712897d-4*t17*t19*t20*t68+1.4447
1957     >         495018085982d-4*t16*t17*t18*t68)-1.2407009817988002d+0*t4
1958     ?         1*(-8.126715947673364d-5*t10*t56*t74*t9+1.625343189534673
1959     @         d-4*t12*t54*t74*t9-5.417810631782243d-5*t4*t56*t74*t8+8.1
1960     1         26715947673364d-5*t10*t54*t74*t8+2.708905315891122d-4*t21
1961     2         *t24*t58*t74+1.625343189534673d-4*t19*t22*t58*t74+8.12671
1962     3         5947673364d-5*t18*t20*t58*t74+2.7089053158911214d-5*t16*t
1963     4         58*t74-2.7089053158911214d-5*t15*t56*t74-1.35445265794556
1964     5         1d-4*t13*t14*t56*t74-1.0835621263564485d-4*t11*t12*t56*t7
1965     6         4+2.7089053158911214d-5*t4*t54*t74+2.708905315891122d-4*t
1966     7         11*t14*t54*t74-1.3544526579455607d-5*t17*t25*t74-6.772263
1967     8         289727805d-5*t17*t23*t24*t74-5.417810631782243d-5*t17*t21
1968     9         *t22*t74-4.063357973836682d-5*t17*t19*t20*t74-2.708905315
1969     :         8911214d-5*t16*t17*t18*t74))*wght+Cmat3(iq,D3_RA_GAA_GAA)
1970            Cmat3(iq,D3_RA_GAA_GAB) = Cmat3(iq,D3_RA_GAA_GAB)
1971            Cmat3(iq,D3_RA_GAA_GBB) = Cmat3(iq,D3_RA_GAA_GBB)
1972            Cmat3(iq,D3_RA_GAB_GAB) = Cmat3(iq,D3_RA_GAB_GAB)
1973            Cmat3(iq,D3_RA_GAB_GBB) = Cmat3(iq,D3_RA_GAB_GBB)
1974            Cmat3(iq,D3_RA_GBB_GBB) = Cmat3(iq,D3_RA_GBB_GBB)
1975            Cmat3(iq,D3_RB_GAA_GAA) = Cmat3(iq,D3_RB_GAA_GAA)
1976            Cmat3(iq,D3_RB_GAA_GAB) = Cmat3(iq,D3_RB_GAA_GAB)
1977            Cmat3(iq,D3_RB_GAA_GBB) = Cmat3(iq,D3_RB_GAA_GBB)
1978            Cmat3(iq,D3_RB_GAB_GAB) = Cmat3(iq,D3_RB_GAB_GAB)
1979            Cmat3(iq,D3_RB_GAB_GBB) = Cmat3(iq,D3_RB_GAB_GBB)
1980            Cmat3(iq,D3_RB_GBB_GBB) = (-9.305257363491002d-1*t27*(-7.975
1981     1         64338870707d-6*gammabb*t24*t37*t91*t95-3.190257355482828d
1982     2         -6*gammabb*t22*t36*t91*t95-7.97564338870707d-7*gammabb*t2
1983     3         0*t91*t95+1.5951286777414142d-6*gammabb*t4*t88*t95+1.5951
1984     4         28677741414d-5*gammabb*t14*t33*t88*t95+9.570772066448484d
1985     5         -6*gammabb*t12*t32*t88*t95+4.785386033224242d-6*gammabb*t
1986     6         10*t31*t88*t95-7.97564338870707d-6*gammabb*t14*t32*t86*t9
1987     7         5-3.190257355482828d-6*gammabb*t12*t31*t86*t95-7.97564338
1988     8         870707d-7*gammabb*t10*t86*t95+7.97564338870707d-6*gammabb
1989     9         *t24*t38*t66*t95+4.785386033224242d-6*gammabb*t22*t37*t66
1990     :         *t95+2.392693016612121d-6*gammabb*t20*t36*t66*t95+7.97564
1991     ;         3388707071d-7*gammabb*t16*t66*t95-1.595128677741414d-6*ga
1992     <         mmabb*t31*t4*t62*t95-3.987821694353535d-6*gammabb*t14*t34
1993     =         *t62*t95-3.190257355482828d-6*gammabb*t12*t33*t62*t95-2.3
1994     >         926930166121207d-6*gammabb*t10*t32*t62*t95-7.975643388707
1995     ?         07d-7*gammabb*t15*t62*t95-6.646369490589226d-7*gammabb*t2
1996     @         4*t35*t39*t95-5.31709559247138d-7*gammabb*t22*t35*t38*t95
1997     1         -3.987821694353535d-7*gammabb*t20*t35*t37*t95-2.658547796
1998     2         23569d-7*gammabb*t16*t35*t36*t95-1.329273898117845d-7*gam
1999     3         mabb*t25*t35*t95-1.4447495018085982d-3*t24*t38*t66*t71-8.
2000     4         668497010851588d-4*t22*t37*t66*t71-4.334248505425794d-4*t
2001     5         20*t36*t66*t71-1.4447495018085982d-4*t16*t66*t71+2.889499
2002     6         0036171964d-4*t31*t4*t64*t71+7.22374750904299d-4*t14*t34*
2003     7         t64*t71+5.778998007234393d-4*t12*t33*t64*t71+4.3342485054
2004     8         25794d-4*t10*t32*t64*t71+1.4447495018085982d-4*t15*t64*t7
2005     9         1-1.4447495018085982d-4*t4*t62*t71-1.4447495018085982d-3*
2006     :         t14*t33*t62*t71-8.668497010851588d-4*t12*t32*t62*t71-4.33
2007     ;         4248505425794d-4*t10*t31*t62*t71+3.611873754521495d-4*t24
2008     <         *t35*t39*t71+2.8894990036171964d-4*t22*t35*t38*t71+2.1671
2009     =         24252712897d-4*t20*t35*t37*t71+1.4447495018085982d-4*t16*
2010     >         t35*t36*t71+7.22374750904299d-5*t25*t35*t71)-1.2407009817
2011     ?         988002d+0*t45*(2.708905315891122d-4*t24*t38*t66*t76+1.625
2012     @         343189534673d-4*t22*t37*t66*t76+8.126715947673364d-5*t20*
2013     1         t36*t66*t76+2.7089053158911214d-5*t16*t66*t76-5.417810631
2014     2         782243d-5*t31*t4*t64*t76-1.354452657945561d-4*t14*t34*t64
2015     3         *t76-1.0835621263564485d-4*t12*t33*t64*t76-8.126715947673
2016     4         364d-5*t10*t32*t64*t76-2.7089053158911214d-5*t15*t64*t76+
2017     5         2.7089053158911214d-5*t4*t62*t76+2.708905315891122d-4*t14
2018     6         *t33*t62*t76+1.625343189534673d-4*t12*t32*t62*t76+8.12671
2019     7         5947673364d-5*t10*t31*t62*t76-6.772263289727805d-5*t24*t3
2020     8         5*t39*t76-5.417810631782243d-5*t22*t35*t38*t76-4.06335797
2021     9         3836682d-5*t20*t35*t37*t76-2.7089053158911214d-5*t16*t35*
2022     :         t36*t76-1.3544526579455607d-5*t25*t35*t76))*wght+Cmat3(iq
2023     ;         ,D3_RB_GBB_GBB)
2024            Cmat3(iq,D3_GAA_GAA_GAA) = Cmat3(iq,D3_GAA_GAA_GAA)-9.305257
2025     1         363491002d-1*t1*(t42*t5*(-1.1963465083060604d-6*t12*t56*t
2026     2         74*t9+2.990866270765151d-6*t14*t54*t74*t9-5.9817325415303
2027     3         01d-7*t10*t56*t74*t8+1.1963465083060604d-6*t12*t54*t74*t8
2028     4         -1.9939108471767675d-7*t4*t56*t74-1.9939108471767675d-6*t
2029     5         11*t14*t56*t74+2.9908662707651507d-7*t10*t54*t74)+t5*(2.9
2030     6         90866270765151d-6*t19*t24*t74*t84+1.1963465083060604d-6*t
2031     7         18*t22*t74*t84+2.9908662707651507d-7*t20*t74*t84-2.990866
2032     8         270765151d-6*t21*t24*t58*t74-1.7945197624590903d-6*t19*t2
2033     9         2*t58*t74-8.972598812295453d-7*t18*t20*t58*t74-2.99086627
2034     :         0765151d-7*t16*t58*t74+4.9847771179419187d-8*t17*t25*t74+
2035     ;         2.4923885589709593d-7*t17*t23*t24*t74+1.9939108471767675d
2036     <         -7*t17*t21*t22*t74+1.4954331353825753d-7*t17*t19*t20*t74+
2037     =         9.969554235883837d-8*t16*t17*t18*t74)-1.4721155704335503d
2038     >         -2*t56*t74*t75+8.126715947673364d-5*t50*t54/rhoa**8)*wght
2039            Cmat3(iq,D3_GAA_GAA_GAB) = Cmat3(iq,D3_GAA_GAA_GAB)
2040            Cmat3(iq,D3_GAA_GAA_GBB) = Cmat3(iq,D3_GAA_GAA_GBB)
2041            Cmat3(iq,D3_GAA_GAB_GAB) = Cmat3(iq,D3_GAA_GAB_GAB)
2042            Cmat3(iq,D3_GAA_GAB_GBB) = Cmat3(iq,D3_GAA_GAB_GBB)
2043            Cmat3(iq,D3_GAA_GBB_GBB) = Cmat3(iq,D3_GAA_GBB_GBB)
2044            Cmat3(iq,D3_GAB_GAB_GAB) = Cmat3(iq,D3_GAB_GAB_GAB)
2045            Cmat3(iq,D3_GAB_GAB_GBB) = Cmat3(iq,D3_GAB_GAB_GBB)
2046            Cmat3(iq,D3_GAB_GBB_GBB) = Cmat3(iq,D3_GAB_GBB_GBB)
2047            Cmat3(iq,D3_GBB_GBB_GBB) = Cmat3(iq,D3_GBB_GBB_GBB)-9.305257
2048     1         363491002d-1*t27*(t28*(2.990866270765151d-6*t24*t37*t76*t
2049     2         91+1.1963465083060604d-6*t22*t36*t76*t91+2.99086627076515
2050     3         07d-7*t20*t76*t91-2.990866270765151d-6*t24*t38*t66*t76-1.
2051     4         7945197624590903d-6*t22*t37*t66*t76-8.972598812295453d-7*
2052     5         t20*t36*t66*t76-2.990866270765151d-7*t16*t66*t76+2.492388
2053     6         5589709593d-7*t24*t35*t39*t76+1.9939108471767675d-7*t22*t
2054     7         35*t38*t76+1.4954331353825753d-7*t20*t35*t37*t76+9.969554
2055     8         235883837d-8*t16*t35*t36*t76+4.9847771179419187d-8*t25*t3
2056     9         5*t76)-1.4721155704335503d-2*t64*t76*t77+t28*t46*(-1.9939
2057     :         108471767675d-7*t4*t64*t76-1.9939108471767675d-6*t14*t33*
2058     ;         t64*t76-1.1963465083060604d-6*t12*t32*t64*t76-5.981732541
2059     <         530301d-7*t10*t31*t64*t76+2.990866270765151d-6*t14*t32*t6
2060     =         2*t76+1.1963465083060604d-6*t12*t31*t62*t76+2.99086627076
2061     >         51507d-7*t10*t62*t76)+8.126715947673364d-5*t51*t62/rhob**
2062     ?         8)*wght
2063          elseif (rhoa.gt.tol_rho.and.rhob.le.tol_rho) then
2064            t1 = rhoa**1.3333333333333333d+0
2065            t2 = param(3)
2066            t3 = 1/rhoa**2.6666666666666666d+0
2067            t4 = 3.6802889260838756d-3*gammaaa*t3
2068            t5 = t4+1.0d+0
2069            t6 = 1.0d+0-1.0d+0/t5
2070            t7 = t6**2.0d+0
2071            t8 = param(4)
2072            t9 = t6**3.0d+0
2073            t10 = param(5)
2074            t11 = t6**4.0d+0
2075            t12 = param(6)
2076            t13 = param(2)
2077            t14 = param(9)
2078            t15 = exp(-t4)
2079            t16 = 1.0d+0-t15
2080            t17 = t16**2.0d+0
2081            t18 = param(10)
2082            t19 = t16**3.0d+0
2083            t20 = param(11)
2084            t21 = t16**4.0d+0
2085            t22 = param(12)
2086            t23 = param(8)
2087            t24 = t8*t9+t2*t7+t12*t6**5.0d+0+t13*t6+t16*t23+t16**5.0d+0*
2088     1         t22+t20*t21+t18*t19+t14*t17+t10*t11+param(7)+param(1)
2089            t25 = rhoa**3.333333333333333d-1
2090            t26 = 1/t5**2
2091            t27 = 1/rhoa**3.6666666666666664d+0
2092            t28 = -3.925641521156134d-2*gammaaa*t10*t26*t27*t9-2.9442311
2093     1         408671007d-2*gammaaa*t26*t27*t7*t8-1.962820760578067d-2*g
2094     2         ammaaa*t2*t26*t27*t6-9.814103802890335d-3*gammaaa*t13*t26
2095     3         *t27-4.9070519014451675d-2*gammaaa*t11*t12*t26*t27-9.8141
2096     4         03802890335d-3*gammaaa*t15*t23*t27-4.9070519014451675d-2*
2097     5         gammaaa*t15*t21*t22*t27-3.925641521156134d-2*gammaaa*t15*
2098     6         t19*t20*t27-2.9442311408671007d-2*gammaaa*t15*t17*t18*t27
2099     7         -1.962820760578067d-2*gammaaa*t14*t15*t16*t27
2100            t29 = 1.4721155704335503d-2*t10*t9+1.1040866778251626d-2*t7*
2101     1         t8+7.360577852167751d-3*t2*t6+3.6802889260838756d-3*t13+1
2102     2         .8401444630419378d-2*t11*t12
2103            t30 = 1/rhoa**6.666666666666666d-1
2104            t31 = gammaaa**2
2105            t32 = 1/t5**4
2106            t33 = 1/rhoa**7.333333333333333d+0
2107            t34 = 1/t5**3
2108            t35 = 1/rhoa**4.666666666666667d+0
2109            t36 = exp(-7.360577852167751d-3*gammaaa*t3)
2110            t37 = 1.4394018910905823d-1*gammaaa*t10*t26*t35*t9-7.7053306
2111     1         76312523d-4*t10*t31*t33*t34*t9+1.9263326690781307d-3*t12*
2112     2         t31*t32*t33*t9+1.0795514183179368d-1*gammaaa*t26*t35*t7*t
2113     3         8-5.778998007234393d-4*t31*t33*t34*t7*t8+5.77899800723439
2114     4         3d-4*t31*t32*t33*t6*t8+1.1557996014468785d-3*t10*t31*t32*
2115     5         t33*t7+7.197009455452912d-2*gammaaa*t2*t26*t35*t6-3.85266
2116     6         53381562614d-4*t2*t31*t33*t34*t6+1.9263326690781307d-3*t1
2117     7         9*t22*t31*t33*t36+1.1557996014468785d-3*t17*t20*t31*t33*t
2118     8         36+5.778998007234393d-4*t16*t18*t31*t33*t36+1.92633266907
2119     9         81307d-4*t14*t31*t33*t36+3.598504727726456d-2*gammaaa*t13
2120     :         *t26*t35+1.799252363863228d-1*gammaaa*t11*t12*t26*t35+3.5
2121     ;         98504727726456d-2*gammaaa*t15*t23*t35+1.799252363863228d-
2122     <         1*gammaaa*t15*t21*t22*t35+1.4394018910905823d-1*gammaaa*t
2123     =         15*t19*t20*t35+1.0795514183179368d-1*gammaaa*t15*t17*t18*
2124     >         t35+7.197009455452912d-2*gammaaa*t14*t15*t16*t35-1.926332
2125     ?         6690781307d-4*t13*t31*t33*t34-9.631663345390654d-4*t11*t1
2126     @         2*t31*t33*t34+1.9263326690781307d-4*t2*t31*t32*t33-9.6316
2127     1         63345390653d-5*t15*t23*t31*t33-4.815831672695327d-4*t15*t
2128     2         21*t22*t31*t33-3.8526653381562614d-4*t15*t19*t20*t31*t33-
2129     3         2.8894990036171964d-4*t15*t17*t18*t31*t33-1.9263326690781
2130     4         307d-4*t14*t15*t16*t31*t33
2131            t38 = 1/rhoa**6.333333333333333d+0
2132            t39 = 2.8894990036171964d-4*gammaaa*t10*t34*t38*t9-7.2237475
2133     1         0904299d-4*gammaaa*t12*t32*t38*t9-3.925641521156134d-2*t1
2134     2         0*t26*t27*t9+2.167124252712897d-4*gammaaa*t34*t38*t7*t8-2
2135     3         .9442311408671007d-2*t26*t27*t7*t8-2.167124252712897d-4*g
2136     4         ammaaa*t32*t38*t6*t8-4.334248505425794d-4*gammaaa*t10*t32
2137     5         *t38*t7+1.4447495018085982d-4*gammaaa*t2*t34*t38*t6-1.962
2138     6         820760578067d-2*t2*t26*t27*t6-7.22374750904299d-4*gammaaa
2139     7         *t19*t22*t36*t38-4.334248505425794d-4*gammaaa*t17*t20*t36
2140     8         *t38-2.167124252712897d-4*gammaaa*t16*t18*t36*t38-7.22374
2141     9         750904299d-5*gammaaa*t14*t36*t38+7.22374750904299d-5*gamm
2142     :         aaa*t13*t34*t38+3.611873754521495d-4*gammaaa*t11*t12*t34*
2143     ;         t38-7.22374750904299d-5*gammaaa*t2*t32*t38+3.611873754521
2144     <         495d-5*gammaaa*t15*t23*t38+1.8059368772607476d-4*gammaaa*
2145     =         t15*t21*t22*t38+1.4447495018085982d-4*gammaaa*t15*t19*t20
2146     >         *t38+1.0835621263564485d-4*gammaaa*t15*t17*t18*t38+7.2237
2147     ?         4750904299d-5*gammaaa*t14*t15*t16*t38-9.814103802890335d-
2148     @         3*t13*t26*t27-4.9070519014451675d-2*t11*t12*t26*t27-9.814
2149     1         103802890335d-3*t15*t23*t27-4.9070519014451675d-2*t15*t21
2150     2         *t22*t27-3.925641521156134d-2*t15*t19*t20*t27-2.944231140
2151     3         8671007d-2*t15*t17*t18*t27-1.962820760578067d-2*t14*t15*t
2152     4         16*t27
2153            t40 = 1.4721155704335503d-2*t10*t26*t3*t9+1.1040866778251626
2154     1         d-2*t26*t3*t7*t8+7.360577852167751d-3*t2*t26*t3*t6+3.6802
2155     2         889260838756d-3*t13*t26*t3+1.8401444630419378d-2*t11*t12*
2156     3         t26*t3+3.6802889260838756d-3*t15*t23*t3+1.840144463041937
2157     4         8d-2*t15*t21*t22*t3+1.4721155704335503d-2*t15*t19*t20*t3+
2158     5         1.1040866778251626d-2*t15*t17*t18*t3+7.360577852167751d-3
2159     6         *t14*t15*t16*t3
2160            t41 = 1/rhoa**5.333333333333333d+0
2161            t42 = 2.708905315891122d-4*t12*t26*t3*t9+8.126715947673364d-
2162     1         5*t26*t3*t6*t8+1.625343189534673d-4*t10*t26*t3*t7+2.70890
2163     2         53158911214d-5*t2*t26*t3
2164            t43 = gammaaa**3
2165            t44 = 1/t5**6
2166            t45 = 1/rhoa**11
2167            t46 = 1/t5**5
2168            t47 = 1/rhoa**8.333333333333334d+0
2169            t48 = 1/rhoa**5.666666666666667d+0
2170            t49 = exp(-1.1040866778251626d-2*gammaaa*t3)
2171            t50 = 1/rhoa**10
2172            t51 = 1/rhoa**9
2173            fnc(iq) = fnc(iq)-9.305257363491002d-1*t1*t24*wght
2174            Amat(iq,D1_RA) = -9.305257363491002d-1*t1*t28*wght-1.2407009
2175     1         817988002d+0*t24*t25*wght+Amat(iq,D1_RA)
2176            Cmat(iq,D1_GAA) = Cmat(iq,D1_GAA)-9.305257363491002d-1*t1*(t
2177     1         26*t29*t3+(3.6802889260838756d-3*t15*t23+1.84014446304193
2178     2         78d-2*t15*t21*t22+1.4721155704335503d-2*t15*t19*t20+1.104
2179     3         0866778251626d-2*t15*t17*t18+7.360577852167751d-3*t14*t15
2180     4         *t16)*t3)*wght
2181            Amat2(iq,D2_RA_RA) = -9.305257363491002d-1*t1*t37*wght-4.135
2182     1         669939329334d-1*t24*t30*wght-2.4814019635976003d+0*t25*t2
2183     2         8*wght+Amat2(iq,D2_RA_RA)
2184            Cmat2(iq,D2_RA_GAA) = -1.2407009817988002d+0*t25*t40*wght-9.
2185     1         305257363491002d-1*t1*t39*wght+Cmat2(iq,D2_RA_GAA)
2186            Cmat2(iq,D2_GAA_GAA) = Cmat2(iq,D2_GAA_GAA)-9.30525736349100
2187     1         2d-1*t1*(t26*t3*t42-7.360577852167751d-3*t29*t34*t41+t3*(
2188     2         2.708905315891122d-4*t19*t22*t3*t36+1.625343189534673d-4*
2189     3         t17*t20*t3*t36+8.126715947673364d-5*t16*t18*t3*t36+2.7089
2190     4         053158911214d-5*t14*t3*t36-1.3544526579455607d-5*t15*t23*
2191     5         t3-6.772263289727805d-5*t15*t21*t22*t3-5.417810631782243d
2192     6         -5*t15*t19*t20*t3-4.063357973836682d-5*t15*t17*t18*t3-2.7
2193     7         089053158911214d-5*t14*t15*t16*t3))*wght
2194            Amat3(iq,D3_RA_RA_RA) = -9.305257363491002d-1*t1*(-6.7172088
2195     1         25089385d-1*gammaaa*t10*t26*t48*t9+8.475863743943775d-3*t
2196     2         10*t31*t34*t47*t9-2.1189659359859436d-2*t12*t31*t32*t47*t
2197     3         9+1.1343137263938945d-4*t12*t43*t45*t46*t9-2.268627452787
2198     4         7887d-5*t10*t32*t43*t45*t9-5.037906618817038d-1*gammaaa*t
2199     5         26*t48*t7*t8+6.356897807957831d-3*t31*t34*t47*t7*t8-1.701
2200     6         4705895908416d-5*t32*t43*t45*t7*t8-6.356897807957831d-3*t
2201     7         31*t32*t47*t6*t8+3.4029411791816827d-5*t43*t45*t46*t6*t8-
2202     8         5.671568631969471d-6*t43*t44*t45*t8-1.2713795615915663d-2
2203     9         *t10*t31*t32*t47*t7+6.805882358363365d-5*t10*t43*t45*t46*
2204     :         t7-5.671568631969472d-5*t12*t43*t44*t45*t7-3.358604412544
2205     ;         6923d-1*gammaaa*t2*t26*t48*t6+4.237931871971887d-3*t2*t31
2206     <         *t34*t47*t6-2.2686274527877887d-5*t10*t43*t44*t45*t6-1.13
2207     =         43137263938943d-5*t2*t32*t43*t45*t6-5.671568631969472d-5*
2208     >         t17*t22*t43*t45*t49-2.2686274527877887d-5*t16*t20*t43*t45
2209     ?         *t49-5.671568631969471d-6*t18*t43*t45*t49-1.6793022062723
2210     @         462d-1*gammaaa*t13*t26*t48-8.396511031361731d-1*gammaaa*t
2211     1         11*t12*t26*t48-1.6793022062723462d-1*gammaaa*t15*t23*t48-
2212     2         8.396511031361731d-1*gammaaa*t15*t21*t22*t48-6.7172088250
2213     3         89385d-1*gammaaa*t15*t19*t20*t48-5.037906618817038d-1*gam
2214     4         maaa*t15*t17*t18*t48-3.3586044125446923d-1*gammaaa*t14*t1
2215     5         5*t16*t48-2.1189659359859436d-2*t19*t22*t31*t36*t47-1.271
2216     6         3795615915663d-2*t17*t20*t31*t36*t47-6.356897807957831d-3
2217     7         *t16*t18*t31*t36*t47-2.1189659359859436d-3*t14*t31*t36*t4
2218     8         7+2.1189659359859436d-3*t13*t31*t34*t47+1.059482967992971
2219     9         8d-2*t11*t12*t31*t34*t47-2.1189659359859436d-3*t2*t31*t32
2220     :         *t47+1.0594829679929718d-3*t15*t23*t31*t47+5.297414839964
2221     ;         859d-3*t15*t21*t22*t31*t47+4.237931871971887d-3*t15*t19*t
2222     <         20*t31*t47+3.1784489039789154d-3*t15*t17*t18*t31*t47+2.11
2223     =         89659359859436d-3*t14*t15*t16*t31*t47+1.1343137263938943d
2224     >         -5*t2*t43*t45*t46+5.671568631969472d-5*t19*t22*t36*t43*t4
2225     ?         5+3.4029411791816827d-5*t17*t20*t36*t43*t45+1.70147058959
2226     @         08414d-5*t16*t18*t36*t43*t45+5.671568631969471d-6*t14*t36
2227     1         *t43*t45-5.671568631969471d-6*t13*t32*t43*t45-2.835784315
2228     2         984736d-5*t11*t12*t32*t43*t45-9.452614386615786d-7*t15*t2
2229     3         3*t43*t45-4.726307193307893d-6*t15*t21*t22*t43*t45-3.7810
2230     4         457546463144d-6*t15*t19*t20*t43*t45-2.8357843159847357d-6
2231     5         *t15*t17*t18*t43*t45-1.8905228773231572d-6*t14*t15*t16*t4
2232     6         3*t45)*wght-3.7221029453964005d+0*t25*t37*wght-1.24070098
2233     7         17988002d+0*t28*t30*wght+2.7571132928862224d-1*t24*wght/r
2234     8         hoa**1.6666666666666669d+0+Amat3(iq,D3_RA_RA_RA)
2235            Cmat3(iq,D3_RA_RA_GAA) = -9.305257363491002d-1*t1*(-4.253676
2236     1         4739771044d-5*t12*t31*t46*t50*t9+8.507352947954206d-6*t10
2237     2         *t31*t32*t50*t9+1.4394018910905823d-1*t10*t26*t35*t9-2.60
2238     3         0549103255476d-3*gammaaa*t10*t33*t34*t9+6.50137275813869d
2239     4         -3*gammaaa*t12*t32*t33*t9+6.380514710965656d-6*t31*t32*t5
2240     5         0*t7*t8+1.0795514183179368d-1*t26*t35*t7*t8-1.95041182744
2241     6         16074d-3*gammaaa*t33*t34*t7*t8-1.2761029421931314d-5*t31*
2242     7         t46*t50*t6*t8+1.9504118274416074d-3*gammaaa*t32*t33*t6*t8
2243     8         +2.1268382369885522d-6*t31*t44*t50*t8-2.552205884386262d-
2244     9         5*t10*t31*t46*t50*t7+2.1268382369885516d-5*t12*t31*t44*t5
2245     :         0*t7+3.900823654883215d-3*gammaaa*t10*t32*t33*t7+8.507352
2246     ;         947954209d-6*t10*t31*t44*t50*t6+4.253676473977103d-6*t2*t
2247     <         31*t32*t50*t6+7.197009455452912d-2*t2*t26*t35*t6-1.300274
2248     =         551627738d-3*gammaaa*t2*t33*t34*t6+2.1268382369885516d-5*
2249     >         t17*t22*t31*t49*t50+8.507352947954209d-6*t16*t20*t31*t49*
2250     ?         t50+2.1268382369885522d-6*t18*t31*t49*t50-4.2536764739771
2251     @         03d-6*t2*t31*t46*t50-2.1268382369885522d-5*t19*t22*t31*t3
2252     1         6*t50-1.276102942193131d-5*t17*t20*t31*t36*t50-6.38051471
2253     2         0965657d-6*t16*t18*t31*t36*t50-2.1268382369885516d-6*t14*
2254     3         t31*t36*t50+2.1268382369885516d-6*t13*t31*t32*t50+1.06341
2255     4         91184942758d-5*t11*t12*t31*t32*t50+3.54473039498092d-7*t1
2256     5         5*t23*t31*t50+1.77236519749046d-6*t15*t21*t22*t31*t50+1.4
2257     6         178921579923678d-6*t15*t19*t20*t31*t50+1.0634191184942761
2258     7         d-6*t15*t17*t18*t31*t50+7.08946078996184d-7*t14*t15*t16*t
2259     8         31*t50+6.50137275813869d-3*gammaaa*t19*t22*t33*t36+3.9008
2260     9         23654883215d-3*gammaaa*t17*t20*t33*t36+1.9504118274416074
2261     :         d-3*gammaaa*t16*t18*t33*t36+6.501372758138692d-4*gammaaa*
2262     ;         t14*t33*t36+3.598504727726456d-2*t13*t26*t35+1.7992523638
2263     <         63228d-1*t11*t12*t26*t35+3.598504727726456d-2*t15*t23*t35
2264     =         +1.799252363863228d-1*t15*t21*t22*t35+1.4394018910905823d
2265     >         -1*t15*t19*t20*t35+1.0795514183179368d-1*t15*t17*t18*t35+
2266     ?         7.197009455452912d-2*t14*t15*t16*t35-6.501372758138692d-4
2267     @         *gammaaa*t13*t33*t34-3.250686379069345d-3*gammaaa*t11*t12
2268     1         *t33*t34+6.501372758138692d-4*gammaaa*t2*t32*t33-3.250686
2269     2         379069346d-4*gammaaa*t15*t23*t33-1.6253431895346726d-3*ga
2270     3         mmaaa*t15*t21*t22*t33-1.300274551627738d-3*gammaaa*t15*t1
2271     4         9*t20*t33-9.752059137208037d-4*gammaaa*t15*t17*t18*t33-6.
2272     5         501372758138692d-4*gammaaa*t14*t15*t16*t33)*wght-4.135669
2273     6         939329334d-1*t30*t40*wght-2.4814019635976003d+0*t25*t39*w
2274     7         ght+Cmat3(iq,D3_RA_RA_GAA)
2275            Cmat3(iq,D3_RA_GAA_GAA) = -9.305257363491002d-1*t1*(1.595128
2276     1         677741414d-5*gammaaa*t12*t46*t51*t9-3.190257355482828d-6*
2277     2         gammaaa*t10*t32*t51*t9+5.778998007234393d-4*t10*t34*t38*t
2278     3         9-1.4447495018085982d-3*t12*t32*t38*t9-2.3926930166121207
2279     4         d-6*gammaaa*t32*t51*t7*t8+4.334248505425794d-4*t34*t38*t7
2280     5         *t8+4.785386033224242d-6*gammaaa*t46*t51*t6*t8-4.33424850
2281     6         5425794d-4*t32*t38*t6*t8-7.97564338870707d-7*gammaaa*t44*
2282     7         t51*t8+9.570772066448484d-6*gammaaa*t10*t46*t51*t7-7.9756
2283     8         4338870707d-6*gammaaa*t12*t44*t51*t7-8.668497010851588d-4
2284     9         *t10*t32*t38*t7-3.190257355482828d-6*gammaaa*t10*t44*t51*
2285     :         t6-1.595128677741414d-6*gammaaa*t2*t32*t51*t6+2.889499003
2286     ;         6171964d-4*t2*t34*t38*t6-7.97564338870707d-6*gammaaa*t17*
2287     <         t22*t49*t51-3.190257355482828d-6*gammaaa*t16*t20*t49*t51-
2288     =         7.97564338870707d-7*gammaaa*t18*t49*t51+1.595128677741414
2289     >         2d-6*gammaaa*t2*t46*t51+7.97564338870707d-6*gammaaa*t19*t
2290     ?         22*t36*t51+4.785386033224242d-6*gammaaa*t17*t20*t36*t51+2
2291     @         .392693016612121d-6*gammaaa*t16*t18*t36*t51+7.97564338870
2292     1         7071d-7*gammaaa*t14*t36*t51-7.97564338870707d-7*gammaaa*t
2293     2         13*t32*t51-3.987821694353535d-6*gammaaa*t11*t12*t32*t51-1
2294     3         .329273898117845d-7*gammaaa*t15*t23*t51-6.646369490589226
2295     4         d-7*gammaaa*t15*t21*t22*t51-5.31709559247138d-7*gammaaa*t
2296     5         15*t19*t20*t51-3.987821694353535d-7*gammaaa*t15*t17*t18*t
2297     6         51-2.65854779623569d-7*gammaaa*t14*t15*t16*t51-1.44474950
2298     7         18085982d-3*t19*t22*t36*t38-8.668497010851588d-4*t17*t20*
2299     8         t36*t38-4.334248505425794d-4*t16*t18*t36*t38-1.4447495018
2300     9         085982d-4*t14*t36*t38+1.4447495018085982d-4*t13*t34*t38+7
2301     :         .22374750904299d-4*t11*t12*t34*t38-1.4447495018085982d-4*
2302     ;         t2*t32*t38+7.22374750904299d-5*t15*t23*t38+3.611873754521
2303     <         495d-4*t15*t21*t22*t38+2.8894990036171964d-4*t15*t19*t20*
2304     =         t38+2.167124252712897d-4*t15*t17*t18*t38+1.44474950180859
2305     >         82d-4*t14*t15*t16*t38)*wght-1.2407009817988002d+0*t25*(-1
2306     ?         .0835621263564485d-4*t10*t34*t41*t9+2.708905315891122d-4*
2307     @         t12*t32*t41*t9-8.126715947673364d-5*t34*t41*t7*t8+8.12671
2308     1         5947673364d-5*t32*t41*t6*t8+1.625343189534673d-4*t10*t32*
2309     2         t41*t7-5.417810631782243d-5*t2*t34*t41*t6+2.7089053158911
2310     3         22d-4*t19*t22*t36*t41+1.625343189534673d-4*t17*t20*t36*t4
2311     4         1+8.126715947673364d-5*t16*t18*t36*t41+2.7089053158911214
2312     5         d-5*t14*t36*t41-2.7089053158911214d-5*t13*t34*t41-1.35445
2313     6         2657945561d-4*t11*t12*t34*t41+2.7089053158911214d-5*t2*t3
2314     7         2*t41-1.3544526579455607d-5*t15*t23*t41-6.772263289727805
2315     8         d-5*t15*t21*t22*t41-5.417810631782243d-5*t15*t19*t20*t41-
2316     9         4.063357973836682d-5*t15*t17*t18*t41-2.7089053158911214d-
2317     :         5*t14*t15*t16*t41)*wght+Cmat3(iq,D3_RA_GAA_GAA)
2318            Cmat3(iq,D3_GAA_GAA_GAA) = Cmat3(iq,D3_GAA_GAA_GAA)-9.305257
2319     1         363491002d-1*t1*(t26*t3*(-1.9939108471767675d-6*t12*t34*t
2320     2         41*t9-5.981732541530301d-7*t34*t41*t6*t8+2.99086627076515
2321     3         07d-7*t32*t41*t8-1.1963465083060604d-6*t10*t34*t41*t7+2.9
2322     4         90866270765151d-6*t12*t32*t41*t7+1.1963465083060604d-6*t1
2323     5         0*t32*t41*t6-1.9939108471767675d-7*t2*t34*t41)+t3*(2.9908
2324     6         66270765151d-6*t17*t22*t41*t49+1.1963465083060604d-6*t16*
2325     7         t20*t41*t49+2.9908662707651507d-7*t18*t41*t49-2.990866270
2326     8         765151d-6*t19*t22*t36*t41-1.7945197624590903d-6*t17*t20*t
2327     9         36*t41-8.972598812295453d-7*t16*t18*t36*t41-2.99086627076
2328     :         5151d-7*t14*t36*t41+4.9847771179419187d-8*t15*t23*t41+2.4
2329     ;         923885589709593d-7*t15*t21*t22*t41+1.9939108471767675d-7*
2330     <         t15*t19*t20*t41+1.4954331353825753d-7*t15*t17*t18*t41+9.9
2331     =         69554235883837d-8*t14*t15*t16*t41)-1.4721155704335503d-2*
2332     >         t34*t41*t42+8.126715947673364d-5*t29*t32/rhoa**8)*wght
2333          elseif (rhoa.le.tol_rho.and.rhob.gt.tol_rho) then
2334            t1 = rhob**1.3333333333333333d+0
2335            t2 = param(3)
2336            t3 = 1/rhob**2.6666666666666666d+0
2337            t4 = 3.6802889260838756d-3*gammabb*t3
2338            t5 = t4+1.0d+0
2339            t6 = 1.0d+0-1.0d+0/t5
2340            t7 = t6**2.0d+0
2341            t8 = param(4)
2342            t9 = t6**3.0d+0
2343            t10 = param(5)
2344            t11 = t6**4.0d+0
2345            t12 = param(6)
2346            t13 = param(2)
2347            t14 = param(9)
2348            t15 = exp(-t4)
2349            t16 = 1.0d+0-t15
2350            t17 = t16**2.0d+0
2351            t18 = param(10)
2352            t19 = t16**3.0d+0
2353            t20 = param(11)
2354            t21 = t16**4.0d+0
2355            t22 = param(12)
2356            t23 = param(8)
2357            t24 = t8*t9+t2*t7+t12*t6**5.0d+0+t13*t6+t16*t23+t16**5.0d+0*
2358     1         t22+t20*t21+t18*t19+t14*t17+t10*t11+param(7)+param(1)
2359            t25 = rhob**3.333333333333333d-1
2360            t26 = 1/t5**2
2361            t27 = 1/rhob**3.6666666666666664d+0
2362            t28 = -3.925641521156134d-2*gammabb*t10*t26*t27*t9-2.9442311
2363     1         408671007d-2*gammabb*t26*t27*t7*t8-1.962820760578067d-2*g
2364     2         ammabb*t2*t26*t27*t6-9.814103802890335d-3*gammabb*t13*t26
2365     3         *t27-4.9070519014451675d-2*gammabb*t11*t12*t26*t27-9.8141
2366     4         03802890335d-3*gammabb*t15*t23*t27-4.9070519014451675d-2*
2367     5         gammabb*t15*t21*t22*t27-3.925641521156134d-2*gammabb*t15*
2368     6         t19*t20*t27-2.9442311408671007d-2*gammabb*t15*t17*t18*t27
2369     7         -1.962820760578067d-2*gammabb*t14*t15*t16*t27
2370            t29 = 1.4721155704335503d-2*t10*t9+1.1040866778251626d-2*t7*
2371     1         t8+7.360577852167751d-3*t2*t6+3.6802889260838756d-3*t13+1
2372     2         .8401444630419378d-2*t11*t12
2373            t30 = 1/rhob**6.666666666666666d-1
2374            t31 = gammabb**2
2375            t32 = 1/t5**4
2376            t33 = 1/rhob**7.333333333333333d+0
2377            t34 = 1/t5**3
2378            t35 = 1/rhob**4.666666666666667d+0
2379            t36 = exp(-7.360577852167751d-3*gammabb*t3)
2380            t37 = 1.4394018910905823d-1*gammabb*t10*t26*t35*t9-7.7053306
2381     1         76312523d-4*t10*t31*t33*t34*t9+1.9263326690781307d-3*t12*
2382     2         t31*t32*t33*t9+1.0795514183179368d-1*gammabb*t26*t35*t7*t
2383     3         8-5.778998007234393d-4*t31*t33*t34*t7*t8+5.77899800723439
2384     4         3d-4*t31*t32*t33*t6*t8+1.1557996014468785d-3*t10*t31*t32*
2385     5         t33*t7+7.197009455452912d-2*gammabb*t2*t26*t35*t6-3.85266
2386     6         53381562614d-4*t2*t31*t33*t34*t6+1.9263326690781307d-3*t1
2387     7         9*t22*t31*t33*t36+1.1557996014468785d-3*t17*t20*t31*t33*t
2388     8         36+5.778998007234393d-4*t16*t18*t31*t33*t36+1.92633266907
2389     9         81307d-4*t14*t31*t33*t36+3.598504727726456d-2*gammabb*t13
2390     :         *t26*t35+1.799252363863228d-1*gammabb*t11*t12*t26*t35+3.5
2391     ;         98504727726456d-2*gammabb*t15*t23*t35+1.799252363863228d-
2392     <         1*gammabb*t15*t21*t22*t35+1.4394018910905823d-1*gammabb*t
2393     =         15*t19*t20*t35+1.0795514183179368d-1*gammabb*t15*t17*t18*
2394     >         t35+7.197009455452912d-2*gammabb*t14*t15*t16*t35-1.926332
2395     ?         6690781307d-4*t13*t31*t33*t34-9.631663345390654d-4*t11*t1
2396     @         2*t31*t33*t34+1.9263326690781307d-4*t2*t31*t32*t33-9.6316
2397     1         63345390653d-5*t15*t23*t31*t33-4.815831672695327d-4*t15*t
2398     2         21*t22*t31*t33-3.8526653381562614d-4*t15*t19*t20*t31*t33-
2399     3         2.8894990036171964d-4*t15*t17*t18*t31*t33-1.9263326690781
2400     4         307d-4*t14*t15*t16*t31*t33
2401            t38 = 1/rhob**6.333333333333333d+0
2402            t39 = 2.8894990036171964d-4*gammabb*t10*t34*t38*t9-7.2237475
2403     1         0904299d-4*gammabb*t12*t32*t38*t9-3.925641521156134d-2*t1
2404     2         0*t26*t27*t9+2.167124252712897d-4*gammabb*t34*t38*t7*t8-2
2405     3         .9442311408671007d-2*t26*t27*t7*t8-2.167124252712897d-4*g
2406     4         ammabb*t32*t38*t6*t8-4.334248505425794d-4*gammabb*t10*t32
2407     5         *t38*t7+1.4447495018085982d-4*gammabb*t2*t34*t38*t6-1.962
2408     6         820760578067d-2*t2*t26*t27*t6-7.22374750904299d-4*gammabb
2409     7         *t19*t22*t36*t38-4.334248505425794d-4*gammabb*t17*t20*t36
2410     8         *t38-2.167124252712897d-4*gammabb*t16*t18*t36*t38-7.22374
2411     9         750904299d-5*gammabb*t14*t36*t38+7.22374750904299d-5*gamm
2412     :         abb*t13*t34*t38+3.611873754521495d-4*gammabb*t11*t12*t34*
2413     ;         t38-7.22374750904299d-5*gammabb*t2*t32*t38+3.611873754521
2414     <         495d-5*gammabb*t15*t23*t38+1.8059368772607476d-4*gammabb*
2415     =         t15*t21*t22*t38+1.4447495018085982d-4*gammabb*t15*t19*t20
2416     >         *t38+1.0835621263564485d-4*gammabb*t15*t17*t18*t38+7.2237
2417     ?         4750904299d-5*gammabb*t14*t15*t16*t38-9.814103802890335d-
2418     @         3*t13*t26*t27-4.9070519014451675d-2*t11*t12*t26*t27-9.814
2419     1         103802890335d-3*t15*t23*t27-4.9070519014451675d-2*t15*t21
2420     2         *t22*t27-3.925641521156134d-2*t15*t19*t20*t27-2.944231140
2421     3         8671007d-2*t15*t17*t18*t27-1.962820760578067d-2*t14*t15*t
2422     4         16*t27
2423            t40 = 1.4721155704335503d-2*t10*t26*t3*t9+1.1040866778251626
2424     1         d-2*t26*t3*t7*t8+7.360577852167751d-3*t2*t26*t3*t6+3.6802
2425     2         889260838756d-3*t13*t26*t3+1.8401444630419378d-2*t11*t12*
2426     3         t26*t3+3.6802889260838756d-3*t15*t23*t3+1.840144463041937
2427     4         8d-2*t15*t21*t22*t3+1.4721155704335503d-2*t15*t19*t20*t3+
2428     5         1.1040866778251626d-2*t15*t17*t18*t3+7.360577852167751d-3
2429     6         *t14*t15*t16*t3
2430            t41 = 1/rhob**5.333333333333333d+0
2431            t42 = 2.708905315891122d-4*t12*t26*t3*t9+8.126715947673364d-
2432     1         5*t26*t3*t6*t8+1.625343189534673d-4*t10*t26*t3*t7+2.70890
2433     2         53158911214d-5*t2*t26*t3
2434            t43 = gammabb**3
2435            t44 = 1/t5**6
2436            t45 = 1/rhob**11
2437            t46 = 1/t5**5
2438            t47 = 1/rhob**8.333333333333334d+0
2439            t48 = 1/rhob**5.666666666666667d+0
2440            t49 = exp(-1.1040866778251626d-2*gammabb*t3)
2441            t50 = 1/rhob**10
2442            t51 = 1/rhob**9
2443            fnc(iq) = fnc(iq)-9.305257363491002d-1*t1*t24*wght
2444            Amat(iq,D1_RB) = -9.305257363491002d-1*t1*t28*wght-1.2407009
2445     1         817988002d+0*t24*t25*wght+Amat(iq,D1_RB)
2446            Cmat(iq,D1_GBB) = Cmat(iq,D1_GBB)-9.305257363491002d-1*t1*(t
2447     1         26*t29*t3+(3.6802889260838756d-3*t15*t23+1.84014446304193
2448     2         78d-2*t15*t21*t22+1.4721155704335503d-2*t15*t19*t20+1.104
2449     3         0866778251626d-2*t15*t17*t18+7.360577852167751d-3*t14*t15
2450     4         *t16)*t3)*wght
2451            Amat2(iq,D2_RB_RB) = -9.305257363491002d-1*t1*t37*wght-4.135
2452     1         669939329334d-1*t24*t30*wght-2.4814019635976003d+0*t25*t2
2453     2         8*wght+Amat2(iq,D2_RB_RB)
2454            Cmat2(iq,D2_RB_GBB) = -1.2407009817988002d+0*t25*t40*wght-9.
2455     1         305257363491002d-1*t1*t39*wght+Cmat2(iq,D2_RB_GBB)
2456            Cmat2(iq,D2_GBB_GBB) = Cmat2(iq,D2_GBB_GBB)-9.30525736349100
2457     1         2d-1*t1*(t26*t3*t42-7.360577852167751d-3*t29*t34*t41+t3*(
2458     2         2.708905315891122d-4*t19*t22*t3*t36+1.625343189534673d-4*
2459     3         t17*t20*t3*t36+8.126715947673364d-5*t16*t18*t3*t36+2.7089
2460     4         053158911214d-5*t14*t3*t36-1.3544526579455607d-5*t15*t23*
2461     5         t3-6.772263289727805d-5*t15*t21*t22*t3-5.417810631782243d
2462     6         -5*t15*t19*t20*t3-4.063357973836682d-5*t15*t17*t18*t3-2.7
2463     7         089053158911214d-5*t14*t15*t16*t3))*wght
2464            Amat3(iq,D3_RB_RB_RB) = -9.305257363491002d-1*t1*(-6.7172088
2465     1         25089385d-1*gammabb*t10*t26*t48*t9+8.475863743943775d-3*t
2466     2         10*t31*t34*t47*t9-2.1189659359859436d-2*t12*t31*t32*t47*t
2467     3         9+1.1343137263938945d-4*t12*t43*t45*t46*t9-2.268627452787
2468     4         7887d-5*t10*t32*t43*t45*t9-5.037906618817038d-1*gammabb*t
2469     5         26*t48*t7*t8+6.356897807957831d-3*t31*t34*t47*t7*t8-1.701
2470     6         4705895908416d-5*t32*t43*t45*t7*t8-6.356897807957831d-3*t
2471     7         31*t32*t47*t6*t8+3.4029411791816827d-5*t43*t45*t46*t6*t8-
2472     8         5.671568631969471d-6*t43*t44*t45*t8-1.2713795615915663d-2
2473     9         *t10*t31*t32*t47*t7+6.805882358363365d-5*t10*t43*t45*t46*
2474     :         t7-5.671568631969472d-5*t12*t43*t44*t45*t7-3.358604412544
2475     ;         6923d-1*gammabb*t2*t26*t48*t6+4.237931871971887d-3*t2*t31
2476     <         *t34*t47*t6-2.2686274527877887d-5*t10*t43*t44*t45*t6-1.13
2477     =         43137263938943d-5*t2*t32*t43*t45*t6-5.671568631969472d-5*
2478     >         t17*t22*t43*t45*t49-2.2686274527877887d-5*t16*t20*t43*t45
2479     ?         *t49-5.671568631969471d-6*t18*t43*t45*t49-1.6793022062723
2480     @         462d-1*gammabb*t13*t26*t48-8.396511031361731d-1*gammabb*t
2481     1         11*t12*t26*t48-1.6793022062723462d-1*gammabb*t15*t23*t48-
2482     2         8.396511031361731d-1*gammabb*t15*t21*t22*t48-6.7172088250
2483     3         89385d-1*gammabb*t15*t19*t20*t48-5.037906618817038d-1*gam
2484     4         mabb*t15*t17*t18*t48-3.3586044125446923d-1*gammabb*t14*t1
2485     5         5*t16*t48-2.1189659359859436d-2*t19*t22*t31*t36*t47-1.271
2486     6         3795615915663d-2*t17*t20*t31*t36*t47-6.356897807957831d-3
2487     7         *t16*t18*t31*t36*t47-2.1189659359859436d-3*t14*t31*t36*t4
2488     8         7+2.1189659359859436d-3*t13*t31*t34*t47+1.059482967992971
2489     9         8d-2*t11*t12*t31*t34*t47-2.1189659359859436d-3*t2*t31*t32
2490     :         *t47+1.0594829679929718d-3*t15*t23*t31*t47+5.297414839964
2491     ;         859d-3*t15*t21*t22*t31*t47+4.237931871971887d-3*t15*t19*t
2492     <         20*t31*t47+3.1784489039789154d-3*t15*t17*t18*t31*t47+2.11
2493     =         89659359859436d-3*t14*t15*t16*t31*t47+1.1343137263938943d
2494     >         -5*t2*t43*t45*t46+5.671568631969472d-5*t19*t22*t36*t43*t4
2495     ?         5+3.4029411791816827d-5*t17*t20*t36*t43*t45+1.70147058959
2496     @         08414d-5*t16*t18*t36*t43*t45+5.671568631969471d-6*t14*t36
2497     1         *t43*t45-5.671568631969471d-6*t13*t32*t43*t45-2.835784315
2498     2         984736d-5*t11*t12*t32*t43*t45-9.452614386615786d-7*t15*t2
2499     3         3*t43*t45-4.726307193307893d-6*t15*t21*t22*t43*t45-3.7810
2500     4         457546463144d-6*t15*t19*t20*t43*t45-2.8357843159847357d-6
2501     5         *t15*t17*t18*t43*t45-1.8905228773231572d-6*t14*t15*t16*t4
2502     6         3*t45)*wght-3.7221029453964005d+0*t25*t37*wght-1.24070098
2503     7         17988002d+0*t28*t30*wght+2.7571132928862224d-1*t24*wght/r
2504     8         hob**1.6666666666666669d+0+Amat3(iq,D3_RB_RB_RB)
2505            Cmat3(iq,D3_RB_RB_GBB) = -9.305257363491002d-1*t1*(-4.253676
2506     1         4739771044d-5*t12*t31*t46*t50*t9+8.507352947954206d-6*t10
2507     2         *t31*t32*t50*t9+1.4394018910905823d-1*t10*t26*t35*t9-2.60
2508     3         0549103255476d-3*gammabb*t10*t33*t34*t9+6.50137275813869d
2509     4         -3*gammabb*t12*t32*t33*t9+6.380514710965656d-6*t31*t32*t5
2510     5         0*t7*t8+1.0795514183179368d-1*t26*t35*t7*t8-1.95041182744
2511     6         16074d-3*gammabb*t33*t34*t7*t8-1.2761029421931314d-5*t31*
2512     7         t46*t50*t6*t8+1.9504118274416074d-3*gammabb*t32*t33*t6*t8
2513     8         +2.1268382369885522d-6*t31*t44*t50*t8-2.552205884386262d-
2514     9         5*t10*t31*t46*t50*t7+2.1268382369885516d-5*t12*t31*t44*t5
2515     :         0*t7+3.900823654883215d-3*gammabb*t10*t32*t33*t7+8.507352
2516     ;         947954209d-6*t10*t31*t44*t50*t6+4.253676473977103d-6*t2*t
2517     <         31*t32*t50*t6+7.197009455452912d-2*t2*t26*t35*t6-1.300274
2518     =         551627738d-3*gammabb*t2*t33*t34*t6+2.1268382369885516d-5*
2519     >         t17*t22*t31*t49*t50+8.507352947954209d-6*t16*t20*t31*t49*
2520     ?         t50+2.1268382369885522d-6*t18*t31*t49*t50-4.2536764739771
2521     @         03d-6*t2*t31*t46*t50-2.1268382369885522d-5*t19*t22*t31*t3
2522     1         6*t50-1.276102942193131d-5*t17*t20*t31*t36*t50-6.38051471
2523     2         0965657d-6*t16*t18*t31*t36*t50-2.1268382369885516d-6*t14*
2524     3         t31*t36*t50+2.1268382369885516d-6*t13*t31*t32*t50+1.06341
2525     4         91184942758d-5*t11*t12*t31*t32*t50+3.54473039498092d-7*t1
2526     5         5*t23*t31*t50+1.77236519749046d-6*t15*t21*t22*t31*t50+1.4
2527     6         178921579923678d-6*t15*t19*t20*t31*t50+1.0634191184942761
2528     7         d-6*t15*t17*t18*t31*t50+7.08946078996184d-7*t14*t15*t16*t
2529     8         31*t50+6.50137275813869d-3*gammabb*t19*t22*t33*t36+3.9008
2530     9         23654883215d-3*gammabb*t17*t20*t33*t36+1.9504118274416074
2531     :         d-3*gammabb*t16*t18*t33*t36+6.501372758138692d-4*gammabb*
2532     ;         t14*t33*t36+3.598504727726456d-2*t13*t26*t35+1.7992523638
2533     <         63228d-1*t11*t12*t26*t35+3.598504727726456d-2*t15*t23*t35
2534     =         +1.799252363863228d-1*t15*t21*t22*t35+1.4394018910905823d
2535     >         -1*t15*t19*t20*t35+1.0795514183179368d-1*t15*t17*t18*t35+
2536     ?         7.197009455452912d-2*t14*t15*t16*t35-6.501372758138692d-4
2537     @         *gammabb*t13*t33*t34-3.250686379069345d-3*gammabb*t11*t12
2538     1         *t33*t34+6.501372758138692d-4*gammabb*t2*t32*t33-3.250686
2539     2         379069346d-4*gammabb*t15*t23*t33-1.6253431895346726d-3*ga
2540     3         mmabb*t15*t21*t22*t33-1.300274551627738d-3*gammabb*t15*t1
2541     4         9*t20*t33-9.752059137208037d-4*gammabb*t15*t17*t18*t33-6.
2542     5         501372758138692d-4*gammabb*t14*t15*t16*t33)*wght-4.135669
2543     6         939329334d-1*t30*t40*wght-2.4814019635976003d+0*t25*t39*w
2544     7         ght+Cmat3(iq,D3_RB_RB_GBB)
2545            Cmat3(iq,D3_RB_GBB_GBB) = -9.305257363491002d-1*t1*(1.595128
2546     1         677741414d-5*gammabb*t12*t46*t51*t9-3.190257355482828d-6*
2547     2         gammabb*t10*t32*t51*t9+5.778998007234393d-4*t10*t34*t38*t
2548     3         9-1.4447495018085982d-3*t12*t32*t38*t9-2.3926930166121207
2549     4         d-6*gammabb*t32*t51*t7*t8+4.334248505425794d-4*t34*t38*t7
2550     5         *t8+4.785386033224242d-6*gammabb*t46*t51*t6*t8-4.33424850
2551     6         5425794d-4*t32*t38*t6*t8-7.97564338870707d-7*gammabb*t44*
2552     7         t51*t8+9.570772066448484d-6*gammabb*t10*t46*t51*t7-7.9756
2553     8         4338870707d-6*gammabb*t12*t44*t51*t7-8.668497010851588d-4
2554     9         *t10*t32*t38*t7-3.190257355482828d-6*gammabb*t10*t44*t51*
2555     :         t6-1.595128677741414d-6*gammabb*t2*t32*t51*t6+2.889499003
2556     ;         6171964d-4*t2*t34*t38*t6-7.97564338870707d-6*gammabb*t17*
2557     <         t22*t49*t51-3.190257355482828d-6*gammabb*t16*t20*t49*t51-
2558     =         7.97564338870707d-7*gammabb*t18*t49*t51+1.595128677741414
2559     >         2d-6*gammabb*t2*t46*t51+7.97564338870707d-6*gammabb*t19*t
2560     ?         22*t36*t51+4.785386033224242d-6*gammabb*t17*t20*t36*t51+2
2561     @         .392693016612121d-6*gammabb*t16*t18*t36*t51+7.97564338870
2562     1         7071d-7*gammabb*t14*t36*t51-7.97564338870707d-7*gammabb*t
2563     2         13*t32*t51-3.987821694353535d-6*gammabb*t11*t12*t32*t51-1
2564     3         .329273898117845d-7*gammabb*t15*t23*t51-6.646369490589226
2565     4         d-7*gammabb*t15*t21*t22*t51-5.31709559247138d-7*gammabb*t
2566     5         15*t19*t20*t51-3.987821694353535d-7*gammabb*t15*t17*t18*t
2567     6         51-2.65854779623569d-7*gammabb*t14*t15*t16*t51-1.44474950
2568     7         18085982d-3*t19*t22*t36*t38-8.668497010851588d-4*t17*t20*
2569     8         t36*t38-4.334248505425794d-4*t16*t18*t36*t38-1.4447495018
2570     9         085982d-4*t14*t36*t38+1.4447495018085982d-4*t13*t34*t38+7
2571     :         .22374750904299d-4*t11*t12*t34*t38-1.4447495018085982d-4*
2572     ;         t2*t32*t38+7.22374750904299d-5*t15*t23*t38+3.611873754521
2573     <         495d-4*t15*t21*t22*t38+2.8894990036171964d-4*t15*t19*t20*
2574     =         t38+2.167124252712897d-4*t15*t17*t18*t38+1.44474950180859
2575     >         82d-4*t14*t15*t16*t38)*wght-1.2407009817988002d+0*t25*(-1
2576     ?         .0835621263564485d-4*t10*t34*t41*t9+2.708905315891122d-4*
2577     @         t12*t32*t41*t9-8.126715947673364d-5*t34*t41*t7*t8+8.12671
2578     1         5947673364d-5*t32*t41*t6*t8+1.625343189534673d-4*t10*t32*
2579     2         t41*t7-5.417810631782243d-5*t2*t34*t41*t6+2.7089053158911
2580     3         22d-4*t19*t22*t36*t41+1.625343189534673d-4*t17*t20*t36*t4
2581     4         1+8.126715947673364d-5*t16*t18*t36*t41+2.7089053158911214
2582     5         d-5*t14*t36*t41-2.7089053158911214d-5*t13*t34*t41-1.35445
2583     6         2657945561d-4*t11*t12*t34*t41+2.7089053158911214d-5*t2*t3
2584     7         2*t41-1.3544526579455607d-5*t15*t23*t41-6.772263289727805
2585     8         d-5*t15*t21*t22*t41-5.417810631782243d-5*t15*t19*t20*t41-
2586     9         4.063357973836682d-5*t15*t17*t18*t41-2.7089053158911214d-
2587     :         5*t14*t15*t16*t41)*wght+Cmat3(iq,D3_RB_GBB_GBB)
2588            Cmat3(iq,D3_GBB_GBB_GBB) = Cmat3(iq,D3_GBB_GBB_GBB)-9.305257
2589     1         363491002d-1*t1*(t26*t3*(-1.9939108471767675d-6*t12*t34*t
2590     2         41*t9-5.981732541530301d-7*t34*t41*t6*t8+2.99086627076515
2591     3         07d-7*t32*t41*t8-1.1963465083060604d-6*t10*t34*t41*t7+2.9
2592     4         90866270765151d-6*t12*t32*t41*t7+1.1963465083060604d-6*t1
2593     5         0*t32*t41*t6-1.9939108471767675d-7*t2*t34*t41)+t3*(2.9908
2594     6         66270765151d-6*t17*t22*t41*t49+1.1963465083060604d-6*t16*
2595     7         t20*t41*t49+2.9908662707651507d-7*t18*t41*t49-2.990866270
2596     8         765151d-6*t19*t22*t36*t41-1.7945197624590903d-6*t17*t20*t
2597     9         36*t41-8.972598812295453d-7*t16*t18*t36*t41-2.99086627076
2598     :         5151d-7*t14*t36*t41+4.9847771179419187d-8*t15*t23*t41+2.4
2599     ;         923885589709593d-7*t15*t21*t22*t41+1.9939108471767675d-7*
2600     <         t15*t19*t20*t41+1.4954331353825753d-7*t15*t17*t18*t41+9.9
2601     =         69554235883837d-8*t14*t15*t16*t41)-1.4721155704335503d-2*
2602     >         t34*t41*t42+8.126715947673364d-5*t29*t32/rhob**8)*wght
2603          endif ! rhoa.gt.tol_rho.and.rhob.gt.tol_rho
2604        endif ! ipol.eq.1
2605      enddo ! iq
2606      end
2607C> @}
2608