1C> \ingroup nwxc 2C> @{ 3C> 4C> \file nwxcm_c_lyp.F 5C> The nwxcm_c_lyp functional 6C> 7C> @} 8C> 9C> \ingroup nwxc_priv 10C> @{ 11C> 12C> \brief Evaluate the nwxcm_c_lyp functional [1] 13C> 14C> \f{eqnarray*}{ 15C> {\it t_1} &=& \rho_\beta+\rho_\alpha\\\\ 16C> {\it t_2} &=& {{1}\over{{\it t_1}}}\\\\ 17C> {\it t_3} &=& {{1}\over{{\it t_1}^{{{1}\over{3}}}}}\\\\ 18C> {\it t_4} &=& {{1}\over{0.349\,{\it t_3}+1.0}}\\\\ 19C> {\it t_5} &=& {{1}\over{{\it t_1}^{{{11}\over{3}}}}}\\\\ 20C> {\it t_6} &=& 0.2533\,{\it t_3}\\\\ 21C> {\it t_7} &=& e^ {- {\it t_6} }\\\\ 22C> {\it t_8} &=& 0.349\,{\it t_3}\,{\it t_4}\\\\ 23C> {\it t_9} &=& {\it t_8}+{\it t_6}-11.0\\\\ 24C> {\it t_{10}} &=& {\it t_8}+{\it t_6}\\\\ 25C> {\it t_{11}} &=& -3.0\,{\it t_{10}}\\\\ 26C> f &=& 1.0\,\left(-0.006491760000000001\,{\it t_5}\,{\it t_4} 27C> \,\left(0.1111111111111111\,\rho_\alpha\,\rho_\beta\,\left({ 28C> \it t_{11}}-1.0\,\rho_\beta\,{\it t_2}\,{\it t_9}+1.0\right) 29C> -\rho_\alpha^2\right)\,{\it t_7}\,\sigma_{\beta\beta} 30C> -0.006491760000000001\,{\it t_5}\,{\it t_4}\, 31C> \left(0.1111111111111111\,\rho_\alpha\,\rho_\beta\,\left(47.0 32C> -7.0\,{\it t_{10}}\right)-1.333333333333333\,{ 33C> \it t_1}^2\right)\,{\it t_7}\,\sigma_{\alpha\beta} 34C> -0.006491760000000001\,{\it t_5}\,{\it t_4}\, 35C> \left(0.1111111111111111\,\rho_\alpha\,\rho_\beta\,\left({ 36C> \it t_{11}}-1.0\,\rho_\alpha\,{\it t_2}\,{\it t_9}+1.0\right) 37C> -\rho_\beta^2\right)\,{\it t_7}\,\sigma_{\alpha\alpha} 38C> -0.2367051431943866\,\rho_\alpha\,\rho_\beta\,{\it t_5} 39C> \,\left(\rho_\beta^{{{8}\over{3}}}+\rho_\alpha^{{{8} 40C> \over{3}}}\right)\,{\it t_4}\,{\it t_7}-0.19672\,\rho_\alpha 41C> \,\rho_\beta\,{\it t_2}\,{\it t_4}\right)\\\\ 42C> g &=& 0\\\\ 43C> G &=& 0.0\\\\ 44C> \f} 45C> 46C> Code generated with Maxima 5.34.0 [2] 47C> driven by autoxc [3]. 48C> 49C> ### References ### 50C> 51C> [1] C Lee, W Yang, RG Parr, Phys.Rev.B 37, 785 (1988) , DOI: 52C> <a href="https://doi.org/10.1103/PhysRevB.37.785 "> 53C> 10.1103/PhysRevB.37.785 </a> 54C> 55C> [2] Maxima, a computer algebra system, 56C> <a href="http://maxima.sourceforge.net/"> 57C> http://maxima.sourceforge.net/</a> 58C> 59C> [3] autoxc, revision 27097 2015-05-08 60C> 61 subroutine nwxcm_c_lyp(param,tol_rho,ipol,nq,wght, 62 +rho,rgamma,fnc,Amat,Cmat) 63c $Id: $ 64#ifdef NWXC_QUAD_PREC 65 implicit real(kind=selected_real_kind(30))(a-h,o-z),integer(i-n) 66 integer, parameter :: rk=selected_real_kind(30) 67#else 68 implicit real(kind=selected_real_kind(15))(a-h,o-z),integer(i-n) 69 integer, parameter :: rk=selected_real_kind(15) 70#endif 71#include "nwxc_param.fh" 72 double precision param(*) !< [Input] Parameters of functional 73 double precision tol_rho !< [Input] The lower limit on the density 74 integer ipol !< [Input] The number of spin channels 75 integer nq !< [Input] The number of points 76 double precision wght !< [Input] The weight of the functional 77 double precision rho(nq,*) !< [Input] The density 78 double precision rgamma(nq,*) !< [Input] The norm of the density 79 !< gradients 80 double precision fnc(nq) !< [Output] The value of the functional 81c 82c Sampling Matrices for the XC Kernel 83c 84 double precision Amat(nq,*) !< [Output] The derivative wrt rho 85 double precision Cmat(nq,*) !< [Output] The derivative wrt rgamma 86 integer iq 87 double precision tmp 88 double precision rhoa,rhob 89 double precision gammaaa,gammaab,gammabb 90 double precision taua,taub 91 double precision nwxcm_heaviside 92 external nwxcm_heaviside 93CDIR$ NOVECTOR 94 do iq = 1, nq 95 if (ipol.eq.1) then 96 rhoa = 0.5d0*rho(iq,R_T) 97 gammaaa = 0.25d0*rgamma(iq,G_TT) 98 if (rhoa.gt.tol_rho) then 99 t1 = 1/rhoa**3.333333333333333d-1 100 t2 = 1/(2.7700148356845083d-1*t1+1.0d+0) 101 t3 = 2.010443432317725d-1*t1 102 t4 = exp(-t3) 103 t5 = rhoa**3.6666666666666664d+0 104 t6 = 1/t5 105 t7 = rhoa**2 106 t8 = -5.333333333333332d+0*t7 107 t9 = 2.7700148356845083d-1*t1*t2 108 t10 = t9+t3 109 t11 = -t7 110 t12 = 7.937005259840998d-1 111 t13 = 3.49d-1*t1*t12+1.0d+0 112 t14 = 1/t13 113 t15 = 1/t13**2 114 t16 = 7.874506561842959d-2 115 t17 = 2.533d-1*t1*t12 116 t18 = 3.49d-1*t1*t12*t14 117 t19 = t18+t17 118 t20 = 4.7d+1-7.0d+0*t19 119 t21 = 3.968502629920499d-1 120 t22 = 1/rhoa**1.3333333333333333d+0 121 t23 = -1.1633333333333332d-1*t14*t21*t22-8.443333333333334d- 122 1 2*t21*t22+1.2788303649853786d-2*t15/rhoa**1.6666666666666 123 2 669d+0 124 t24 = exp(-t17) 125 t25 = t18+t17-1.1d+1 126 t26 = -5.0d-1*t25-3.0d+0*t19+1.0d+0 127 t27 = 1.111111111111111d-1*rhoa*t26 128 t28 = -3.5d+0*t23 129 t29 = 1/rhoa 130 t30 = 1/rhoa**5 131 t31 = t8+1.111111111111111d-1*t20*t7 132 t32 = 1.111111111111111d-1*t26*t7+t11 133 t33 = rhoa**4.666666666666667d+0 134 fnc(iq) = 1.0d+0*(-1.0223881343581931d-3*gammaaa*t2*t4*t6*(1 135 1 .111111111111111d-1*t7*(-5.0d-1*(t9+t3-1.1d+1)-3.0d+0*t10 136 2 +1.0d+0)+t11)-5.111940671790965d-4*gammaaa*t2*t4*t6*(t8+1 137 3 .111111111111111d-1*(4.7d+1-7.0d+0*t10)*t7)-3.72787240661 138 4 23486d-2*rhoa*t2*t4-9.836d-2*rhoa*t2)*wght+fnc(iq) 139 Amat(iq,D1_RA) = 1.0d+0*(t14*t16*t6*(-6.491760000000001d-3*g 140 1 ammaaa*t24*(1.111111111111111d-1*(2.5d-1*t25*t29+t28)*t7+ 141 2 t27-2*rhoa)-6.491760000000001d-3*gammaaa*t24*(1.111111111 142 3 111111d-1*(t28-2.5d-1*t25*t29)*t7+t27)-6.491760000000001d 143 4 -3*gammaaa*t24*(-7.777777777777777d-1*t23*t7+1.1111111111 144 5 11111d-1*rhoa*t20-5.333333333333332d+0*rhoa)-1.1046240015 145 6 73804d+0*t24*t5)+3.937253280921478d-2*t14*(1.735837716758 146 7 835d+0*t24*t33+4.760624000000001d-2*gammaaa*t24*t32+2.380 147 8 3120000000005d-2*gammaaa*t24*t31)/t33+3.125d-2*t14*t30*(- 148 9 3.9971608514092083d-2*t24*t33-1.0962418720000001d-3*gamma 149 : aa*t24*t32-5.481209360000001d-4*gammaaa*t24*t31)+3.125d-2 150 ; *t15*t30*(-5.507339664989395d-2*t24*t33-1.51041616d-3*gam 151 < maaa*t24*t32-7.5520808d-4*gammaaa*t24*t31)-4.540977653965 152 = 4695d-3*t1*t15-4.918d-2*t14)*wght+Amat(iq,D1_RA) 153 Cmat(iq,D1_GAA) = Cmat(iq,D1_GAA)-6.491760000000001d-3*t14*t 154 1 16*t24*t32*t6*wght 155 Cmat(iq,D1_GAB) = Cmat(iq,D1_GAB)-6.491760000000001d-3*t14*t 156 1 16*t24*t31*t6*wght 157 endif ! rhoa.gt.tol_rho 158 else ! ipol.eq.1 159 rhoa = rho(iq,R_A) 160 rhob = rho(iq,R_B) 161 gammaaa = rgamma(iq,G_AA) 162 gammaab = rgamma(iq,G_AB) 163 gammabb = rgamma(iq,G_BB) 164 if (rhoa.gt.tol_rho.and.rhob.gt.tol_rho) then 165 t1 = rhob+rhoa 166 t2 = 1/t1 167 t3 = 1/t1**3.333333333333333d-1 168 t4 = 3.49d-1*t3+1.0d+0 169 t5 = 1/t4 170 t6 = 1/t1**3.6666666666666664d+0 171 t7 = rhoa**2.6666666666666666d+0 172 t8 = rhob**2.6666666666666666d+0 173 t9 = t8+t7 174 t10 = 2.533d-1*t3 175 t11 = exp(-t10) 176 t12 = t1**2 177 t13 = 3.49d-1*t3*t5 178 t14 = t13+t10 179 t15 = 4.7d+1-7.0d+0*t14 180 t16 = 1.111111111111111d-1*rhoa*rhob*t15-1.333333333333333d+ 181 1 0*t12 182 t17 = t13+t10-1.1d+1 183 t18 = -3.0d+0*t14 184 t19 = -1.0d+0*rhoa*t17*t2+t18+1.0d+0 185 t20 = 1.111111111111111d-1*rhoa*rhob*t19-rhob**2 186 t21 = -1.0d+0*rhob*t17*t2+t18+1.0d+0 187 t22 = 1.111111111111111d-1*rhoa*rhob*t21-rhoa**2 188 t23 = 1/t4**2 189 t24 = -2.288509333333333d-2*rhoa*rhob*t23/t1**2.333333333333 190 1 3334d+0 191 t25 = 1/t12 192 t26 = 1.9672d-1*rhoa*rhob*t25*t5 193 t27 = -2.666666666666666d+0*t1 194 t28 = 1/t1**1.3333333333333333d+0 195 t29 = -1.1633333333333332d-1*t28*t5-8.443333333333334d-2*t28 196 1 +4.060033333333332d-2*t23/t1**1.6666666666666669d+0 197 t30 = -7.777777777777777d-1*rhoa*rhob*t29 198 t31 = -3.0d+0*t29 199 t32 = -1.0d+0*rhoa*t2*t29 200 t33 = 1.0d+0*rhoa*t17*t25 201 t34 = -1.0d+0*t17*t2 202 t35 = -1.0d+0*rhob*t2*t29 203 t36 = 1.0d+0*rhob*t17*t25 204 t37 = 1/t1**5 205 t38 = t23*t37*(-2.7536698324946973d-2*rhoa*rhob*t11*t9-7.552 206 1 0808d-4*gammabb*t11*t22-7.5520808d-4*gammaaa*t11*t20-7.55 207 2 20808d-4*gammaab*t11*t16) 208 t39 = t37*t5*(-1.9985804257046041d-2*rhoa*rhob*t11*t9-5.4812 209 1 09360000001d-4*gammabb*t11*t22-5.481209360000001d-4*gamma 210 2 aa*t11*t20-5.481209360000001d-4*gammaab*t11*t16) 211 t40 = t5*(8.679188583794175d-1*rhoa*rhob*t11*t9+2.3803120000 212 1 000005d-2*gammabb*t11*t22+2.3803120000000005d-2*gammaaa*t 213 2 11*t20+2.3803120000000005d-2*gammaab*t11*t16)/t1**4.66666 214 3 6666666667d+0 215 fnc(iq) = 1.0d+0*(-2.367051431943866d-1*rhoa*rhob*t11*t5*t6* 216 1 t9-6.491760000000001d-3*gammabb*t11*t22*t5*t6-6.491760000 217 2 000001d-3*gammaaa*t11*t20*t5*t6-6.491760000000001d-3*gamm 218 3 aab*t11*t16*t5*t6-1.9672d-1*rhoa*rhob*t2*t5)*wght+fnc(iq) 219 Amat(iq,D1_RA) = 1.0d+0*(t5*t6*(-2.367051431943866d-1*rhob*t 220 1 11*t9-6.312137151850309d-1*rhob*t11*t7-6.491760000000001d 221 2 -3*gammabb*t11*(1.111111111111111d-1*rhoa*rhob*(t36+t35+t 222 3 31)+1.111111111111111d-1*rhob*t21-2*rhoa)-6.4917600000000 223 4 01d-3*gammaaa*t11*(1.111111111111111d-1*rhoa*rhob*(t34+t3 224 5 3+t32+t31)+1.111111111111111d-1*rhob*t19)-6.4917600000000 225 6 01d-3*gammaab*t11*(t30+t27+1.111111111111111d-1*rhob*t15) 226 7 )-1.9672d-1*rhob*t2*t5+t40+t39+t38+t26+t24)*wght+Amat(iq, 227 8 D1_RA) 228 Amat(iq,D1_RB) = 1.0d+0*(t5*t6*(-2.367051431943866d-1*rhoa*t 229 1 11*t9-6.312137151850309d-1*rhoa*t11*t8-6.491760000000001d 230 2 -3*gammabb*t11*(1.111111111111111d-1*rhoa*rhob*(t36+t35+t 231 3 34+t31)+1.111111111111111d-1*rhoa*t21)-6.491760000000001d 232 4 -3*gammaaa*t11*(1.111111111111111d-1*rhoa*rhob*(t33+t32+t 233 5 31)+1.111111111111111d-1*rhoa*t19-2*rhob)-6.4917600000000 234 6 01d-3*gammaab*t11*(t30+t27+1.111111111111111d-1*rhoa*t15) 235 7 )-1.9672d-1*rhoa*t2*t5+t40+t39+t38+t26+t24)*wght+Amat(iq, 236 8 D1_RB) 237 Cmat(iq,D1_GAA) = Cmat(iq,D1_GAA)-6.491760000000001d-3*t11*t 238 1 20*t5*t6*wght 239 Cmat(iq,D1_GAB) = Cmat(iq,D1_GAB)-6.491760000000001d-3*t11*t 240 1 16*t5*t6*wght 241 Cmat(iq,D1_GBB) = Cmat(iq,D1_GBB)-6.491760000000001d-3*t11*t 242 1 22*t5*t6*wght 243 elseif (rhoa.gt.tol_rho.and.rhob.le.tol_rho) then 244 fnc(iq) = fnc(iq) 245 Amat(iq,D1_RA) = Amat(iq,D1_RA) 246 Cmat(iq,D1_GAA) = Cmat(iq,D1_GAA) 247 elseif (rhoa.le.tol_rho.and.rhob.gt.tol_rho) then 248 fnc(iq) = fnc(iq) 249 Amat(iq,D1_RB) = Amat(iq,D1_RB) 250 Cmat(iq,D1_GBB) = Cmat(iq,D1_GBB) 251 endif ! rhoa.gt.tol_rho.and.rhob.gt.tol_rho 252 endif ! ipol.eq.1 253 enddo ! iq 254 end 255C> 256C> \brief Evaluate the nwxcm_c_lyp functional [1] 257C> 258C> \f{eqnarray*}{ 259C> {\it t_1} &=& \rho_\beta+\rho_\alpha\\\\ 260C> {\it t_2} &=& {{1}\over{{\it t_1}}}\\\\ 261C> {\it t_3} &=& {{1}\over{{\it t_1}^{{{1}\over{3}}}}}\\\\ 262C> {\it t_4} &=& {{1}\over{0.349\,{\it t_3}+1.0}}\\\\ 263C> {\it t_5} &=& {{1}\over{{\it t_1}^{{{11}\over{3}}}}}\\\\ 264C> {\it t_6} &=& 0.2533\,{\it t_3}\\\\ 265C> {\it t_7} &=& e^ {- {\it t_6} }\\\\ 266C> {\it t_8} &=& 0.349\,{\it t_3}\,{\it t_4}\\\\ 267C> {\it t_9} &=& {\it t_8}+{\it t_6}-11.0\\\\ 268C> {\it t_{10}} &=& {\it t_8}+{\it t_6}\\\\ 269C> {\it t_{11}} &=& -3.0\,{\it t_{10}}\\\\ 270C> f &=& 1.0\,\left(-0.006491760000000001\,{\it t_5}\,{\it t_4} 271C> \,\left(0.1111111111111111\,\rho_\alpha\,\rho_\beta\,\left({ 272C> \it t_{11}}-1.0\,\rho_\beta\,{\it t_2}\,{\it t_9}+1.0\right) 273C> -\rho_\alpha^2\right)\,{\it t_7}\,\sigma_{\beta\beta} 274C> -0.006491760000000001\,{\it t_5}\,{\it t_4}\, 275C> \left(0.1111111111111111\,\rho_\alpha\,\rho_\beta\,\left(47.0 276C> -7.0\,{\it t_{10}}\right)-1.333333333333333\,{ 277C> \it t_1}^2\right)\,{\it t_7}\,\sigma_{\alpha\beta} 278C> -0.006491760000000001\,{\it t_5}\,{\it t_4}\, 279C> \left(0.1111111111111111\,\rho_\alpha\,\rho_\beta\,\left({ 280C> \it t_{11}}-1.0\,\rho_\alpha\,{\it t_2}\,{\it t_9}+1.0\right) 281C> -\rho_\beta^2\right)\,{\it t_7}\,\sigma_{\alpha\alpha} 282C> -0.2367051431943866\,\rho_\alpha\,\rho_\beta\,{\it t_5} 283C> \,\left(\rho_\beta^{{{8}\over{3}}}+\rho_\alpha^{{{8} 284C> \over{3}}}\right)\,{\it t_4}\,{\it t_7}-0.19672\,\rho_\alpha 285C> \,\rho_\beta\,{\it t_2}\,{\it t_4}\right)\\\\ 286C> g &=& 0\\\\ 287C> G &=& 0.0\\\\ 288C> \f} 289C> 290C> Code generated with Maxima 5.34.0 [2] 291C> driven by autoxc [3]. 292C> 293C> ### References ### 294C> 295C> [1] C Lee, W Yang, RG Parr, Phys.Rev.B 37, 785 (1988) , DOI: 296C> <a href="https://doi.org/10.1103/PhysRevB.37.785 "> 297C> 10.1103/PhysRevB.37.785 </a> 298C> 299C> [2] Maxima, a computer algebra system, 300C> <a href="http://maxima.sourceforge.net/"> 301C> http://maxima.sourceforge.net/</a> 302C> 303C> [3] autoxc, revision 27097 2015-05-08 304C> 305 subroutine nwxcm_c_lyp_d2(param,tol_rho,ipol,nq,wght, 306 +rho,rgamma,fnc,Amat,Amat2,Cmat,Cmat2) 307c $Id: $ 308#ifdef NWXC_QUAD_PREC 309 implicit real(kind=selected_real_kind(30))(a-h,o-z),integer(i-n) 310 integer, parameter :: rk=selected_real_kind(30) 311#else 312 implicit real(kind=selected_real_kind(15))(a-h,o-z),integer(i-n) 313 integer, parameter :: rk=selected_real_kind(15) 314#endif 315#include "nwxc_param.fh" 316 double precision param(*) !< [Input] Parameters of functional 317 double precision tol_rho !< [Input] The lower limit on the density 318 integer ipol !< [Input] The number of spin channels 319 integer nq !< [Input] The number of points 320 double precision wght !< [Input] The weight of the functional 321 double precision rho(nq,*) !< [Input] The density 322 double precision rgamma(nq,*) !< [Input] The norm of the density 323 !< gradients 324 double precision fnc(nq) !< [Output] The value of the functional 325c 326c Sampling Matrices for the XC Kernel 327c 328 double precision Amat(nq,*) !< [Output] The derivative wrt rho 329 double precision Cmat(nq,*) !< [Output] The derivative wrt rgamma 330c 331c Sampling Matrices for the XC Kernel 332c 333 double precision Amat2(nq,*) !< [Output] The 2nd derivative wrt rho 334 double precision Cmat2(nq,*) !< [Output] The 2nd derivative wrt rgamma 335 !< and possibly rho 336 integer iq 337 double precision tmp 338 double precision rhoa,rhob 339 double precision gammaaa,gammaab,gammabb 340 double precision taua,taub 341 double precision nwxcm_heaviside 342 external nwxcm_heaviside 343CDIR$ NOVECTOR 344 do iq = 1, nq 345 if (ipol.eq.1) then 346 rhoa = 0.5d0*rho(iq,R_T) 347 gammaaa = 0.25d0*rgamma(iq,G_TT) 348 if (rhoa.gt.tol_rho) then 349 t1 = 1/rhoa**3.333333333333333d-1 350 t2 = 1/(2.7700148356845083d-1*t1+1.0d+0) 351 t3 = 2.010443432317725d-1*t1 352 t4 = exp(-t3) 353 t5 = rhoa**3.6666666666666664d+0 354 t6 = 1/t5 355 t7 = rhoa**2 356 t8 = -5.333333333333332d+0*t7 357 t9 = 2.7700148356845083d-1*t1*t2 358 t10 = t9+t3 359 t11 = -t7 360 t12 = 7.937005259840998d-1 361 t13 = 3.49d-1*t1*t12+1.0d+0 362 t14 = 1/t13 363 t15 = 1.9842513149602492d-1 364 t16 = 1/t13**2 365 t17 = 7.874506561842959d-2 366 t18 = 2.533d-1*t1*t12 367 t19 = 3.49d-1*t1*t12*t14 368 t20 = t19+t18 369 t21 = 4.7d+1-7.0d+0*t20 370 t22 = 1/rhoa**1.6666666666666669d+0 371 t23 = 3.968502629920499d-1 372 t24 = 1/rhoa**1.3333333333333333d+0 373 t25 = -1.1633333333333332d-1*t14*t23*t24-8.443333333333334d- 374 1 2*t23*t24+1.2788303649853786d-2*t16*t22 375 t26 = -7.777777777777777d-1*t25*t7+1.111111111111111d-1*rhoa 376 1 *t21-5.333333333333332d+0*rhoa 377 t27 = exp(-t18) 378 t28 = t19+t18-1.1d+1 379 t29 = -5.0d-1*t28-3.0d+0*t20+1.0d+0 380 t30 = 1.111111111111111d-1*rhoa*t29 381 t31 = -3.5d+0*t25 382 t32 = 1/rhoa 383 t33 = t31-2.5d-1*t28*t32 384 t34 = 1.111111111111111d-1*t33*t7+t30 385 t35 = 2.5d-1*t28*t32+t31 386 t36 = 1.111111111111111d-1*t35*t7+t30-2*rhoa 387 t37 = -1.104624001573804d+0*t27*t5-6.491760000000001d-3*gamm 388 1 aaa*t27*t36-6.491760000000001d-3*gammaaa*t27*t34-6.491760 389 2 000000001d-3*gammaaa*t26*t27 390 t38 = 1/rhoa**5 391 t39 = t8+1.111111111111111d-1*t21*t7 392 t40 = 1.111111111111111d-1*t29*t7+t11 393 t41 = rhoa**4.666666666666667d+0 394 t42 = -5.507339664989395d-2*t27*t41-1.51041616d-3*gammaaa*t2 395 1 7*t40-7.5520808d-4*gammaaa*t27*t39 396 t43 = -3.9971608514092083d-2*t27*t41-1.0962418720000001d-3*g 397 1 ammaaa*t27*t40-5.481209360000001d-4*gammaaa*t27*t39 398 t44 = 3.937253280921477d-2 399 t45 = 1/t41 400 t46 = 1.735837716758835d+0*t27*t41+4.760624000000001d-2*gamm 401 1 aaa*t27*t40+2.3803120000000005d-2*gammaaa*t27*t39 402 t47 = 1/t13**3 403 t48 = -5.324598382222221d-3*t17*t22*t47 404 t49 = 7.568296089942449d-3*t16*t24 405 t50 = -4.577018666666666d-2*t15*t16*t24 406 t51 = -1.5555555555555553d+0*rhoa*t25 407 t52 = rhoa**2.6666666666666666d+0 408 t53 = rhoa**2.3333333333333334d+0 409 t54 = 1/t53 410 t55 = 1.551111111111111d-1*t14*t15*t54+1.1257777777777778d-1 411 1 *t15*t54-1.2788303649853786d-2*t16/t52+1.1807930277777773 412 2 d-3*t47/rhoa**3 413 t56 = -7.777777777777777d-1*t55*t7 414 t57 = -5.481209360000001d-4*gammaaa*t23*t24*t26*t27 415 t58 = -3.5d+0*t55 416 t59 = 1/t7 417 t60 = -5.481209360000001d-4*gammaaa*t23*t24*t27*t34 418 t61 = -5.481209360000001d-4*gammaaa*t23*t24*t27*t36 419 t62 = -5.329547801878944d-2*t23*t27*t53 420 t63 = -1.9985804257046041d-2*t12*t27*t53 421 t64 = -3.6666666666666664d+0*t14*t37*t44*t45 422 t65 = rhoa**3.3333333333333337d+0 423 t66 = 3.125d-2*t14*t38*(-1.6874680727699207d-3*t12*t27*t65-9 424 1 .326708653288152d-2*t27*t5-9.255935539253335d-5*gammaaa*t 425 2 23*t24*t27*t40-4.6279677696266674d-5*gammaaa*t23*t24*t27* 426 3 t39-5.481209360000001d-4*gammaaa*t27*t36-5.48120936000000 427 4 1d-4*gammaaa*t27*t34-5.481209360000001d-4*gammaaa*t26*t27 428 5 ) 429 t67 = t14*t44*t45*(7.328128227583548d-2*t12*t27*t65+4.050288 430 1 005770614d+0*t27*t5+4.0195535306666674d-3*gammaaa*t23*t24 431 2 *t27*t40+2.0097767653333337d-3*gammaaa*t23*t24*t27*t39+2. 432 3 3803120000000005d-2*gammaaa*t27*t36+2.3803120000000005d-2 433 4 *gammaaa*t27*t34+2.3803120000000005d-2*gammaaa*t26*t27) 434 t68 = 1.2401570718501566d-2 435 t69 = 1/rhoa**6.333333333333333d+0 436 t70 = 2.3266666666666663d-1*t42*t47*t68*t69 437 t71 = 1.1633333333333332d-1*t16*t43*t68*t69 438 t72 = 1/rhoa**6 439 t73 = -7.8125d-2*t14*t43*t72 440 t74 = -2.3333333333333334d+0*t14*t44*t46/rhoa**5.66666666666 441 1 6667d+0 442 t75 = 3.125d-2*t16*t38*(-2.3250152285696893d-3*t12*t27*t65-1 443 1 .2850459218308585d-1*t27*t5-1.2752947110933333d-4*gammaaa 444 2 *t23*t24*t27*t40-6.376473555466666d-5*gammaaa*t23*t24*t27 445 3 *t39+1.1633333333333332d-1*t37-7.5520808d-4*gammaaa*t27*t 446 4 36-7.5520808d-4*gammaaa*t27*t34-7.5520808d-4*gammaaa*t26* 447 5 t27) 448 t76 = 1.5625d-2*t16*(1.1633333333333332d-1*t46-5*t42)*t72 449 t77 = -2.36002525d-5*t16*t27*t38*t40 450 t78 = -1.7128779250000004d-5*t14*t27*t38*t40 451 t79 = 2.3803120000000005d-2*t14*t27*t40*t44*t45 452 fnc(iq) = 1.0d+0*(-1.0223881343581931d-3*gammaaa*t2*t4*t6*(1 453 1 .111111111111111d-1*t7*(-5.0d-1*(t9+t3-1.1d+1)-3.0d+0*t10 454 2 +1.0d+0)+t11)-5.111940671790965d-4*gammaaa*t2*t4*t6*(t8+1 455 3 .111111111111111d-1*(4.7d+1-7.0d+0*t10)*t7)-3.72787240661 456 4 23486d-2*rhoa*t2*t4-9.836d-2*rhoa*t2)*wght+fnc(iq) 457 Amat(iq,D1_RA) = 1.0d+0*(t14*t17*t37*t6+t14*t44*t45*t46+3.12 458 1 5d-2*t14*t38*t43+3.125d-2*t16*t38*t42-2.288509333333333d- 459 2 2*t1*t15*t16-4.918d-2*t14)*wght+Amat(iq,D1_RA) 460 Cmat(iq,D1_GAA) = Cmat(iq,D1_GAA)-6.491760000000001d-3*t14*t 461 1 17*t27*t40*t6*wght 462 Cmat(iq,D1_GAB) = Cmat(iq,D1_GAB)-6.491760000000001d-3*t14*t 463 1 17*t27*t39*t6*wght 464 Amat2(iq,D2_RA_RA) = 1.0d+0*(t76+t75+t74+t73+t71+t70+t14*t17 465 1 *t6*(-6.491760000000001d-3*gammaaa*t27*(1.111111111111111 466 2 d-1*(2.5d-1*t28*t59+t58-5.0d-1*t25*t32)*t7+2.222222222222 467 3 222d-1*rhoa*t33)-6.491760000000001d-3*gammaaa*t27*(1.1111 468 4 11111111111d-1*(-2.5d-1*t28*t59+t58+5.0d-1*t25*t32)*t7+2. 469 5 222222222222222d-1*rhoa*t35-2)+t63+t62+t61+t60+t57-6.4917 470 6 60000000001d-3*gammaaa*t27*(t56+t51-2.666666666666666d+0) 471 7 -2.3144502890117796d+0*t27*t52)+t67+t66+t64+t50+t49+t48+4 472 8 .918d-2*t14*t32)*wght+Amat2(iq,D2_RA_RA) 473 Amat2(iq,D2_RA_RB) = 1.0d+0*(t76+t75+t74+t73+t71+t70+t14*t17 474 1 *t6*(-1.298352d-2*gammaaa*t27*(-3.8888888888888884d-1*t55 475 2 *t7+1.111111111111111d-1*rhoa*t35+1.111111111111111d-1*rh 476 3 oa*t33+1.111111111111111d-1*t29)+t63+t62+t61+t60+t57-6.49 477 4 1760000000001d-3*gammaaa*t27*(t56+t51+1.111111111111111d- 478 5 1*t21-2.666666666666666d+0)-1.735837716758835d+0*t27*t52) 479 6 +t67+t66+t64+t50+t49+t48-4.918d-2*t14*t32)*wght+Amat2(iq, 480 7 D2_RA_RB) 481 Cmat2(iq,D2_RA_GAA) = 1.0d+0*(t79+t78+t77-6.491760000000001d 482 1 -3*t14*t17*t27*t34*t6)*wght+Cmat2(iq,D2_RA_GAA) 483 Cmat2(iq,D2_RA_GAB) = 1.0d+0*(-6.491760000000001d-3*t14*t17* 484 1 t26*t27*t6+2.3803120000000005d-2*t14*t27*t39*t44*t45-2.36 485 2 002525d-5*t16*t27*t38*t39-1.7128779250000004d-5*t14*t27*t 486 3 38*t39)*wght+Cmat2(iq,D2_RA_GAB) 487 Cmat2(iq,D2_RA_GBB) = 1.0d+0*(t79+t78+t77-6.491760000000001d 488 1 -3*t14*t17*t27*t36*t6)*wght+Cmat2(iq,D2_RA_GBB) 489 Cmat2(iq,D2_GAA_GAA) = Cmat2(iq,D2_GAA_GAA) 490 Cmat2(iq,D2_GAA_GAB) = Cmat2(iq,D2_GAA_GAB) 491 Cmat2(iq,D2_GAA_GBB) = Cmat2(iq,D2_GAA_GBB) 492 Cmat2(iq,D2_GAB_GAB) = Cmat2(iq,D2_GAB_GAB) 493 endif ! rhoa.gt.tol_rho 494 else ! ipol.eq.1 495 rhoa = rho(iq,R_A) 496 rhob = rho(iq,R_B) 497 gammaaa = rgamma(iq,G_AA) 498 gammaab = rgamma(iq,G_AB) 499 gammabb = rgamma(iq,G_BB) 500 if (rhoa.gt.tol_rho.and.rhob.gt.tol_rho) then 501 t1 = rhob+rhoa 502 t2 = 1/t1 503 t3 = 1/t1**3.333333333333333d-1 504 t4 = 3.49d-1*t3+1.0d+0 505 t5 = 1/t4 506 t6 = 1/t1**3.6666666666666664d+0 507 t7 = rhoa**2.6666666666666666d+0 508 t8 = rhob**2.6666666666666666d+0 509 t9 = t8+t7 510 t10 = 2.533d-1*t3 511 t11 = exp(-t10) 512 t12 = t1**2 513 t13 = 3.49d-1*t3*t5 514 t14 = t13+t10 515 t15 = 4.7d+1-7.0d+0*t14 516 t16 = 1.111111111111111d-1*rhoa*rhob*t15-1.333333333333333d+ 517 1 0*t12 518 t17 = t13+t10-1.1d+1 519 t18 = -3.0d+0*t14 520 t19 = -1.0d+0*rhoa*t17*t2+t18+1.0d+0 521 t20 = 1.111111111111111d-1*rhoa*rhob*t19-rhob**2 522 t21 = -1.0d+0*rhob*t17*t2+t18+1.0d+0 523 t22 = 1.111111111111111d-1*rhoa*rhob*t21-rhoa**2 524 t23 = 1/t1**2.3333333333333334d+0 525 t24 = 1/t4**2 526 t25 = -2.288509333333333d-2*rhoa*rhob*t23*t24 527 t26 = 1/t12 528 t27 = 1.9672d-1*rhoa*rhob*t26*t5 529 t28 = -2.666666666666666d+0*t1 530 t29 = 1/t1**1.3333333333333333d+0 531 t30 = -1.1633333333333332d-1*t29*t5-8.443333333333334d-2*t29 532 1 +4.060033333333332d-2*t24/t1**1.6666666666666669d+0 533 t31 = -7.777777777777777d-1*rhoa*rhob*t30 534 t32 = t31+t28+1.111111111111111d-1*rhob*t15 535 t33 = -3.0d+0*t30 536 t34 = -1.0d+0*rhoa*t2*t30 537 t35 = 1.0d+0*rhoa*t17*t26 538 t36 = -1.0d+0*t17*t2 539 t37 = t36+t35+t34+t33 540 t38 = 1.111111111111111d-1*rhoa*rhob*t37+1.111111111111111d- 541 1 1*rhob*t19 542 t39 = -1.0d+0*rhob*t2*t30 543 t40 = 1.0d+0*rhob*t17*t26 544 t41 = t40+t39+t33 545 t42 = 1.111111111111111d-1*rhoa*rhob*t41+1.111111111111111d- 546 1 1*rhob*t21-2*rhoa 547 t43 = -2.367051431943866d-1*rhob*t11*t9-6.312137151850309d-1 548 1 *rhob*t11*t7-6.491760000000001d-3*gammabb*t11*t42-6.49176 549 2 0000000001d-3*gammaaa*t11*t38-6.491760000000001d-3*gammaa 550 3 b*t11*t32 551 t44 = 1/t1**5 552 t45 = -2.7536698324946973d-2*rhoa*rhob*t11*t9-7.5520808d-4*g 553 1 ammabb*t11*t22-7.5520808d-4*gammaaa*t11*t20-7.5520808d-4* 554 2 gammaab*t11*t16 555 t46 = t24*t44*t45 556 t47 = -1.9985804257046041d-2*rhoa*rhob*t11*t9-5.481209360000 557 1 001d-4*gammabb*t11*t22-5.481209360000001d-4*gammaaa*t11*t 558 2 20-5.481209360000001d-4*gammaab*t11*t16 559 t48 = t44*t47*t5 560 t49 = 1/t1**4.666666666666667d+0 561 t50 = 8.679188583794175d-1*rhoa*rhob*t11*t9+2.38031200000000 562 1 05d-2*gammabb*t11*t22+2.3803120000000005d-2*gammaaa*t11*t 563 2 20+2.3803120000000005d-2*gammaab*t11*t16 564 t51 = t49*t5*t50 565 t52 = t31+t28+1.111111111111111d-1*rhoa*t15 566 t53 = t35+t34+t33 567 t54 = 1.111111111111111d-1*rhoa*rhob*t53+1.111111111111111d- 568 1 1*rhoa*t19-2*rhob 569 t55 = t40+t39+t36+t33 570 t56 = 1.111111111111111d-1*rhoa*rhob*t55+1.111111111111111d- 571 1 1*rhoa*t21 572 t57 = -2.367051431943866d-1*rhoa*t11*t9-6.312137151850309d-1 573 1 *rhoa*t11*t8-6.491760000000001d-3*gammabb*t11*t56-6.49176 574 2 0000000001d-3*gammaaa*t11*t54-6.491760000000001d-3*gammaa 575 3 b*t11*t52 576 t58 = 1/t4**3 577 t59 = -5.324598382222221d-3*rhoa*rhob*t58*t6 578 t60 = 7.628364444444444d-2*rhoa*rhob*t24/t1**3.3333333333333 579 1 337d+0 580 t61 = 1/t1**3 581 t62 = -3.9344d-1*rhoa*rhob*t5*t61 582 t63 = -3.6666666666666664d+0*t43*t49*t5 583 t64 = -1.9985804257046041d-2*rhob*t11*t29*t9 584 t65 = -5.329547801878944d-2*rhob*t11*t29*t7 585 t66 = 9.446344222222218d-3*t58*t61+1.551111111111111d-1*t23* 586 1 t5-8.120066666666665d-2*t24/t1**2.6666666666666666d+0+1.1 587 2 257777777777778d-1*t23 588 t67 = -7.777777777777777d-1*rhoa*rhob*t66 589 t68 = -3.0d+0*t66 590 t69 = -1.0d+0*rhob*t2*t66 591 t70 = 2.0d+0*rhob*t26*t30 592 t71 = -2.0d+0*rhob*t17*t61 593 t72 = -1.0d+0*rhoa*t2*t66 594 t73 = 2.0d+0*rhoa*t26*t30 595 t74 = -2.0d+0*t2*t30 596 t75 = -2.0d+0*rhoa*t17*t61 597 t76 = 2.0d+0*t17*t26 598 t77 = -5.481209360000001d-4*gammaab*t11*t29*t32 599 t78 = -5.481209360000001d-4*gammaaa*t11*t29*t38 600 t79 = -5.481209360000001d-4*gammabb*t11*t29*t42 601 t80 = 1/t1**6.333333333333333d+0 602 t81 = 2.3266666666666663d-1*t45*t58*t80 603 t82 = 1.1633333333333332d-1*t24*t47*t80 604 t83 = 1/t1**6 605 t84 = -5*t47*t5*t83 606 t85 = -4.666666666666667d+0*t5*t50/t1**5.666666666666667d+0 607 t86 = -1.6874680727699207d-3*rhoa*rhob*t11*t29*t9 608 t87 = -4.6279677696266674d-5*gammaab*t11*t16*t29 609 t88 = -4.6279677696266674d-5*gammaaa*t11*t20*t29 610 t89 = -4.6279677696266674d-5*gammabb*t11*t22*t29 611 t90 = 7.328128227583548d-2*rhoa*rhob*t11*t29*t9 612 t91 = 2.0097767653333337d-3*gammaab*t11*t16*t29 613 t92 = 2.0097767653333337d-3*gammaaa*t11*t20*t29 614 t93 = 2.0097767653333337d-3*gammabb*t11*t22*t29 615 t94 = -2.3250152285696893d-3*rhoa*rhob*t11*t29*t9 616 t95 = -6.376473555466666d-5*gammaab*t11*t16*t29 617 t96 = -6.376473555466666d-5*gammaaa*t11*t20*t29 618 t97 = -6.376473555466666d-5*gammabb*t11*t22*t29 619 t98 = 1.1633333333333332d-1*t43 620 t99 = t24*(1.1633333333333332d-1*t50-5*t45)*t83 621 t100 = -1.0d+0*t2*t30 622 t101 = 1.0d+0*t17*t26 623 t102 = t44*t5*(-1.9985804257046041d-2*rhoa*t11*t9+t89+t88+t8 624 1 7+t86-5.329547801878943d-2*rhoa*t11*t8-5.481209360000001d 625 2 -4*gammabb*t11*t56-5.481209360000001d-4*gammaaa*t11*t54-5 626 3 .481209360000001d-4*gammaab*t11*t52) 627 t103 = t49*t5*(t93+t92+t91+t90+8.679188583794175d-1*rhoa*t11 628 1 *t9+2.3144502890117796d+0*rhoa*t11*t8+2.3803120000000005d 629 2 -2*gammabb*t11*t56+2.3803120000000005d-2*gammaaa*t11*t54+ 630 3 2.3803120000000005d-2*gammaab*t11*t52) 631 t104 = -7.343119553319191d-2*rhoa*t11*t8 632 t105 = -2.7536698324946973d-2*rhoa*t11*t9 633 t106 = -7.5520808d-4*gammaab*t11*t52 634 t107 = -7.5520808d-4*gammaaa*t11*t54 635 t108 = -7.5520808d-4*gammabb*t11*t56 636 t109 = -7.5520808d-4*t11*t20*t24*t44 637 t110 = -5.481209360000001d-4*t11*t20*t44*t5 638 t111 = 2.3803120000000005d-2*t11*t20*t49*t5 639 t112 = -7.5520808d-4*t11*t16*t24*t44 640 t113 = -5.481209360000001d-4*t11*t16*t44*t5 641 t114 = 2.3803120000000005d-2*t11*t16*t49*t5 642 t115 = -7.5520808d-4*t11*t22*t24*t44 643 t116 = -5.481209360000001d-4*t11*t22*t44*t5 644 t117 = 2.3803120000000005d-2*t11*t22*t49*t5 645 fnc(iq) = 1.0d+0*(-2.367051431943866d-1*rhoa*rhob*t11*t5*t6* 646 1 t9-6.491760000000001d-3*gammabb*t11*t22*t5*t6-6.491760000 647 2 000001d-3*gammaaa*t11*t20*t5*t6-6.491760000000001d-3*gamm 648 3 aab*t11*t16*t5*t6-1.9672d-1*rhoa*rhob*t2*t5)*wght+fnc(iq) 649 Amat(iq,D1_RA) = 1.0d+0*(t43*t5*t6+t51-1.9672d-1*rhob*t2*t5+ 650 1 t48+t46+t27+t25)*wght+Amat(iq,D1_RA) 651 Amat(iq,D1_RB) = 1.0d+0*(t5*t57*t6+t51-1.9672d-1*rhoa*t2*t5+ 652 1 t48+t46+t27+t25)*wght+Amat(iq,D1_RB) 653 Cmat(iq,D1_GAA) = Cmat(iq,D1_GAA)-6.491760000000001d-3*t11*t 654 1 20*t5*t6*wght 655 Cmat(iq,D1_GAB) = Cmat(iq,D1_GAB)-6.491760000000001d-3*t11*t 656 1 16*t5*t6*wght 657 Cmat(iq,D1_GBB) = Cmat(iq,D1_GBB)-6.491760000000001d-3*t11*t 658 1 22*t5*t6*wght 659 Amat2(iq,D2_RA_RA) = 1.0d+0*(t99+t24*t44*(t98+t97+t96+t95+t9 660 1 4-2.7536698324946973d-2*rhob*t11*t9-7.343119553319191d-2* 661 2 rhob*t11*t7-7.5520808d-4*gammabb*t11*t42-7.5520808d-4*gam 662 3 maaa*t11*t38-7.5520808d-4*gammaab*t11*t32)+t49*t5*(t93+t9 663 4 2+t91+t90+8.679188583794175d-1*rhob*t11*t9+2.314450289011 664 5 7796d+0*rhob*t11*t7+2.3803120000000005d-2*gammabb*t11*t42 665 6 +2.3803120000000005d-2*gammaaa*t11*t38+2.3803120000000005 666 7 d-2*gammaab*t11*t32)+t44*t5*(-1.9985804257046041d-2*rhob* 667 8 t11*t9+t89+t88+t87+t86-5.329547801878943d-2*rhob*t11*t7-5 668 9 .481209360000001d-4*gammabb*t11*t42-5.481209360000001d-4* 669 : gammaaa*t11*t38-5.481209360000001d-4*gammaab*t11*t32)+t85 670 ; +t84+t82+t81+t5*t6*(t79+t78+t77-6.491760000000001d-3*gamm 671 < aaa*t11*(1.111111111111111d-1*rhoa*rhob*(t76+t75+t74+t73+ 672 = t72+t68)+2.222222222222222d-1*rhob*t37)-6.491760000000001 673 > d-3*gammabb*t11*(1.111111111111111d-1*rhoa*rhob*(t71+t70+ 674 ? t69+t68)+2.222222222222222d-1*rhob*t41-2)-6.4917600000000 675 @ 01d-3*gammaab*t11*(t67-1.5555555555555553d+0*rhob*t30-2.6 676 1 66666666666666d+0)+t65+t64-2.3144502890117796d+0*rhoa**1. 677 2 6666666666666669d+0*rhob*t11)+t63+t62+t60+t59+3.9344d-1*r 678 3 hob*t26*t5-4.577018666666666d-2*rhob*t23*t24)*wght+Amat2( 679 4 iq,D2_RA_RA) 680 Amat2(iq,D2_RA_RB) = 1.0d+0*(t99+t24*t44*(t98+t97+t96+t95+t9 681 1 4+t108+t107+t106+t105+t104)+t5*t6*(-2.367051431943866d-1* 682 2 t11*t9-6.312137151850309d-1*t11*t8+t79+t78+t77-6.49176000 683 3 0000001d-3*gammaaa*t11*(1.111111111111111d-1*rhoa*rhob*(t 684 4 75+t73+t72+t68+t101+t100)+1.111111111111111d-1*rhob*t53+1 685 5 .111111111111111d-1*rhoa*t37+1.111111111111111d-1*t19)-6. 686 6 491760000000001d-3*gammabb*t11*(1.111111111111111d-1*rhoa 687 7 *rhob*(t71+t70+t69+t68+t101+t100)+1.111111111111111d-1*rh 688 8 ob*t55+1.111111111111111d-1*rhoa*t41+1.111111111111111d-1 689 9 *t21)-6.312137151850309d-1*t11*t7-6.491760000000001d-3*ga 690 : mmaab*t11*(t67-7.777777777777777d-1*rhob*t30-7.7777777777 691 ; 77777d-1*rhoa*t30+1.111111111111111d-1*t15-2.666666666666 692 < 666d+0)+t65+t64)+t85+t84+t82+t81+t63+t62+t60+t59+(1.9672d 693 = -1*rhob+1.9672d-1*rhoa)*t26*t5-1.9672d-1*t2*t5+(-2.288509 694 > 333333333d-2*rhob-2.288509333333333d-2*rhoa)*t23*t24+t103 695 ? +t102)*wght+Amat2(iq,D2_RA_RB) 696 Amat2(iq,D2_RB_RB) = 1.0d+0*(t99+t24*t44*(t97+t96+t95+t94+1. 697 1 1633333333333332d-1*t57+t108+t107+t106+t105+t104)+t5*t6*( 698 2 -1.9985804257046041d-2*rhoa*t11*t29*t9-5.329547801878944d 699 3 -2*rhoa*t11*t29*t8-6.491760000000001d-3*gammabb*t11*(1.11 700 4 1111111111111d-1*rhoa*rhob*(t76+t74+t71+t70+t69+t68)+2.22 701 5 2222222222222d-1*rhoa*t55)-6.491760000000001d-3*gammaaa*t 702 6 11*(1.111111111111111d-1*rhoa*rhob*(t75+t73+t72+t68)+2.22 703 7 2222222222222d-1*rhoa*t53-2)-6.491760000000001d-3*gammaab 704 8 *t11*(t67-1.5555555555555553d+0*rhoa*t30-2.66666666666666 705 9 6d+0)-5.481209360000001d-4*gammabb*t11*t29*t56-5.48120936 706 : 0000001d-4*gammaaa*t11*t29*t54-5.481209360000001d-4*gamma 707 ; ab*t11*t29*t52-2.3144502890117796d+0*rhoa*rhob**1.6666666 708 < 666666669d+0*t11)+t85+t84+t82+t81+t62+t60+t59-3.666666666 709 = 6666664d+0*t49*t5*t57+3.9344d-1*rhoa*t26*t5-4.57701866666 710 > 6666d-2*rhoa*t23*t24+t103+t102)*wght+Amat2(iq,D2_RB_RB) 711 Cmat2(iq,D2_RA_GAA) = 1.0d+0*(-6.491760000000001d-3*t11*t38* 712 1 t5*t6+t111+t110+t109)*wght+Cmat2(iq,D2_RA_GAA) 713 Cmat2(iq,D2_RA_GAB) = 1.0d+0*(-6.491760000000001d-3*t11*t32* 714 1 t5*t6+t114+t113+t112)*wght+Cmat2(iq,D2_RA_GAB) 715 Cmat2(iq,D2_RA_GBB) = 1.0d+0*(-6.491760000000001d-3*t11*t42* 716 1 t5*t6+t117+t116+t115)*wght+Cmat2(iq,D2_RA_GBB) 717 Cmat2(iq,D2_RB_GAA) = 1.0d+0*(-6.491760000000001d-3*t11*t5*t 718 1 54*t6+t111+t110+t109)*wght+Cmat2(iq,D2_RB_GAA) 719 Cmat2(iq,D2_RB_GAB) = 1.0d+0*(-6.491760000000001d-3*t11*t5*t 720 1 52*t6+t114+t113+t112)*wght+Cmat2(iq,D2_RB_GAB) 721 Cmat2(iq,D2_RB_GBB) = 1.0d+0*(-6.491760000000001d-3*t11*t5*t 722 1 56*t6+t117+t116+t115)*wght+Cmat2(iq,D2_RB_GBB) 723 Cmat2(iq,D2_GAA_GAA) = Cmat2(iq,D2_GAA_GAA) 724 Cmat2(iq,D2_GAA_GAB) = Cmat2(iq,D2_GAA_GAB) 725 Cmat2(iq,D2_GAA_GBB) = Cmat2(iq,D2_GAA_GBB) 726 Cmat2(iq,D2_GAB_GAB) = Cmat2(iq,D2_GAB_GAB) 727 Cmat2(iq,D2_GAB_GBB) = Cmat2(iq,D2_GAB_GBB) 728 Cmat2(iq,D2_GBB_GBB) = Cmat2(iq,D2_GBB_GBB) 729 elseif (rhoa.gt.tol_rho.and.rhob.le.tol_rho) then 730 fnc(iq) = fnc(iq) 731 Amat(iq,D1_RA) = Amat(iq,D1_RA) 732 Cmat(iq,D1_GAA) = Cmat(iq,D1_GAA) 733 Amat2(iq,D2_RA_RA) = Amat2(iq,D2_RA_RA) 734 Cmat2(iq,D2_RA_GAA) = Cmat2(iq,D2_RA_GAA) 735 Cmat2(iq,D2_GAA_GAA) = Cmat2(iq,D2_GAA_GAA) 736 elseif (rhoa.le.tol_rho.and.rhob.gt.tol_rho) then 737 fnc(iq) = fnc(iq) 738 Amat(iq,D1_RB) = Amat(iq,D1_RB) 739 Cmat(iq,D1_GBB) = Cmat(iq,D1_GBB) 740 Amat2(iq,D2_RB_RB) = Amat2(iq,D2_RB_RB) 741 Cmat2(iq,D2_RB_GBB) = Cmat2(iq,D2_RB_GBB) 742 Cmat2(iq,D2_GBB_GBB) = Cmat2(iq,D2_GBB_GBB) 743 endif ! rhoa.gt.tol_rho.and.rhob.gt.tol_rho 744 endif ! ipol.eq.1 745 enddo ! iq 746 end 747C> 748C> \brief Evaluate the nwxcm_c_lyp functional [1] 749C> 750C> \f{eqnarray*}{ 751C> {\it t_1} &=& \rho_\beta+\rho_\alpha\\\\ 752C> {\it t_2} &=& {{1}\over{{\it t_1}}}\\\\ 753C> {\it t_3} &=& {{1}\over{{\it t_1}^{{{1}\over{3}}}}}\\\\ 754C> {\it t_4} &=& {{1}\over{0.349\,{\it t_3}+1.0}}\\\\ 755C> {\it t_5} &=& {{1}\over{{\it t_1}^{{{11}\over{3}}}}}\\\\ 756C> {\it t_6} &=& 0.2533\,{\it t_3}\\\\ 757C> {\it t_7} &=& e^ {- {\it t_6} }\\\\ 758C> {\it t_8} &=& 0.349\,{\it t_3}\,{\it t_4}\\\\ 759C> {\it t_9} &=& {\it t_8}+{\it t_6}-11.0\\\\ 760C> {\it t_{10}} &=& {\it t_8}+{\it t_6}\\\\ 761C> {\it t_{11}} &=& -3.0\,{\it t_{10}}\\\\ 762C> f &=& 1.0\,\left(-0.006491760000000001\,{\it t_5}\,{\it t_4} 763C> \,\left(0.1111111111111111\,\rho_\alpha\,\rho_\beta\,\left({ 764C> \it t_{11}}-1.0\,\rho_\beta\,{\it t_2}\,{\it t_9}+1.0\right) 765C> -\rho_\alpha^2\right)\,{\it t_7}\,\sigma_{\beta\beta} 766C> -0.006491760000000001\,{\it t_5}\,{\it t_4}\, 767C> \left(0.1111111111111111\,\rho_\alpha\,\rho_\beta\,\left(47.0 768C> -7.0\,{\it t_{10}}\right)-1.333333333333333\,{ 769C> \it t_1}^2\right)\,{\it t_7}\,\sigma_{\alpha\beta} 770C> -0.006491760000000001\,{\it t_5}\,{\it t_4}\, 771C> \left(0.1111111111111111\,\rho_\alpha\,\rho_\beta\,\left({ 772C> \it t_{11}}-1.0\,\rho_\alpha\,{\it t_2}\,{\it t_9}+1.0\right) 773C> -\rho_\beta^2\right)\,{\it t_7}\,\sigma_{\alpha\alpha} 774C> -0.2367051431943866\,\rho_\alpha\,\rho_\beta\,{\it t_5} 775C> \,\left(\rho_\beta^{{{8}\over{3}}}+\rho_\alpha^{{{8} 776C> \over{3}}}\right)\,{\it t_4}\,{\it t_7}-0.19672\,\rho_\alpha 777C> \,\rho_\beta\,{\it t_2}\,{\it t_4}\right)\\\\ 778C> g &=& 0\\\\ 779C> G &=& 0.0\\\\ 780C> \f} 781C> 782C> Code generated with Maxima 5.34.0 [2] 783C> driven by autoxc [3]. 784C> 785C> ### References ### 786C> 787C> [1] C Lee, W Yang, RG Parr, Phys.Rev.B 37, 785 (1988) , DOI: 788C> <a href="https://doi.org/10.1103/PhysRevB.37.785 "> 789C> 10.1103/PhysRevB.37.785 </a> 790C> 791C> [2] Maxima, a computer algebra system, 792C> <a href="http://maxima.sourceforge.net/"> 793C> http://maxima.sourceforge.net/</a> 794C> 795C> [3] autoxc, revision 27097 2015-05-08 796C> 797 subroutine nwxcm_c_lyp_d3(param,tol_rho,ipol,nq,wght, 798 +rho,rgamma,fnc,Amat,Amat2,Amat3, 799 +Cmat,Cmat2,Cmat3) 800c $Id: $ 801#ifdef NWXC_QUAD_PREC 802 implicit real(kind=selected_real_kind(30))(a-h,o-z),integer(i-n) 803 integer, parameter :: rk=selected_real_kind(30) 804#else 805 implicit real(kind=selected_real_kind(15))(a-h,o-z),integer(i-n) 806 integer, parameter :: rk=selected_real_kind(15) 807#endif 808#include "nwxc_param.fh" 809 double precision param(*) !< [Input] Parameters of functional 810 double precision tol_rho !< [Input] The lower limit on the density 811 integer ipol !< [Input] The number of spin channels 812 integer nq !< [Input] The number of points 813 double precision wght !< [Input] The weight of the functional 814 double precision rho(nq,*) !< [Input] The density 815 double precision rgamma(nq,*) !< [Input] The norm of the density 816 !< gradients 817 double precision fnc(nq) !< [Output] The value of the functional 818c 819c Sampling Matrices for the XC Kernel 820c 821 double precision Amat(nq,*) !< [Output] The derivative wrt rho 822 double precision Cmat(nq,*) !< [Output] The derivative wrt rgamma 823c 824c Sampling Matrices for the XC Kernel 825c 826 double precision Amat2(nq,*) !< [Output] The 2nd derivative wrt rho 827 double precision Cmat2(nq,*) !< [Output] The 2nd derivative wrt rgamma 828 !< and possibly rho 829c 830c Sampling Matrices for the XC Kernel 831c 832 double precision Amat3(nq,*) !< [Output] The 3rd derivative wrt rho 833 double precision Cmat3(nq,*) !< [Output] The 3rd derivative wrt rgamma 834 !< and possibly rho 835 integer iq 836 double precision tmp 837 double precision rhoa,rhob 838 double precision gammaaa,gammaab,gammabb 839 double precision taua,taub 840 double precision nwxcm_heaviside 841 external nwxcm_heaviside 842CDIR$ NOVECTOR 843 do iq = 1, nq 844 if (ipol.eq.1) then 845 rhoa = 0.5d0*rho(iq,R_T) 846 gammaaa = 0.25d0*rgamma(iq,G_TT) 847 if (rhoa.gt.tol_rho) then 848 t1 = 1/rhoa**3.333333333333333d-1 849 t2 = 1/(2.7700148356845083d-1*t1+1.0d+0) 850 t3 = 2.010443432317725d-1*t1 851 t4 = exp(-t3) 852 t5 = rhoa**3.6666666666666664d+0 853 t6 = 1/t5 854 t7 = rhoa**2 855 t8 = -5.333333333333332d+0*t7 856 t9 = 2.7700148356845083d-1*t1*t2 857 t10 = t9+t3 858 t11 = -t7 859 t12 = 7.937005259840998d-1 860 t13 = 3.49d-1*t1*t12+1.0d+0 861 t14 = 1/t13 862 t15 = 1.9842513149602492d-1 863 t16 = 1/t13**2 864 t17 = 7.874506561842959d-2 865 t18 = 2.533d-1*t1*t12 866 t19 = 3.49d-1*t1*t12*t14 867 t20 = t19+t18 868 t21 = 4.7d+1-7.0d+0*t20 869 t22 = 3.149802624737183d-1 870 t23 = rhoa**1.6666666666666669d+0 871 t24 = 1/t23 872 t25 = 3.968502629920499d-1 873 t26 = rhoa**1.3333333333333333d+0 874 t27 = 1/t26 875 t28 = -1.1633333333333332d-1*t14*t25*t27-8.443333333333334d- 876 1 2*t25*t27+4.060033333333332d-2*t16*t22*t24 877 t29 = -7.777777777777777d-1*t28*t7+1.111111111111111d-1*rhoa 878 1 *t21-5.333333333333332d+0*rhoa 879 t30 = exp(-t18) 880 t31 = t19+t18-1.1d+1 881 t32 = -5.0d-1*t31-3.0d+0*t20+1.0d+0 882 t33 = 1.111111111111111d-1*rhoa*t32 883 t34 = -3.5d+0*t28 884 t35 = 1/rhoa 885 t36 = t34-2.5d-1*t31*t35 886 t37 = 1.111111111111111d-1*t36*t7+t33 887 t38 = 2.5d-1*t31*t35+t34 888 t39 = 1.111111111111111d-1*t38*t7+t33-2*rhoa 889 t40 = -1.104624001573804d+0*t30*t5-6.491760000000001d-3*gamm 890 1 aaa*t30*t39-6.491760000000001d-3*gammaaa*t30*t37-6.491760 891 2 000000001d-3*gammaaa*t29*t30 892 t41 = 1/rhoa**5 893 t42 = t8+1.111111111111111d-1*t21*t7 894 t43 = 1.111111111111111d-1*t32*t7+t11 895 t44 = rhoa**4.666666666666667d+0 896 t45 = -5.507339664989395d-2*t30*t44-1.51041616d-3*gammaaa*t3 897 1 0*t43-7.5520808d-4*gammaaa*t30*t42 898 t46 = -3.9971608514092083d-2*t30*t44-1.0962418720000001d-3*g 899 1 ammaaa*t30*t43-5.481209360000001d-4*gammaaa*t30*t42 900 t47 = 3.937253280921477d-2 901 t48 = 1/t44 902 t49 = 1.735837716758835d+0*t30*t44+4.760624000000001d-2*gamm 903 1 aaa*t30*t43+2.3803120000000005d-2*gammaaa*t30*t42 904 t50 = 1/t13**3 905 t51 = -5.324598382222221d-3*t17*t24*t50 906 t52 = 9.921256574801247d-2 907 t53 = 7.628364444444444d-2*t16*t27*t52 908 t54 = -4.577018666666666d-2*t15*t16*t27 909 t55 = -1.5555555555555553d+0*rhoa*t28 910 t56 = 1/rhoa**3 911 t57 = 1.5749013123685918d-1 912 t58 = rhoa**2.6666666666666666d+0 913 t59 = 1/t58 914 t60 = rhoa**2.3333333333333334d+0 915 t61 = 1/t60 916 t62 = 1.551111111111111d-1*t14*t15*t61+1.1257777777777778d-1 917 1 *t15*t61-8.120066666666665d-2*t16*t57*t59+1.1807930277777 918 2 773d-3*t50*t56 919 t63 = -7.777777777777777d-1*t62*t7 920 t64 = t63+t55-2.666666666666666d+0 921 t65 = -5.481209360000001d-4*gammaaa*t25*t27*t29*t30 922 t66 = -3.5d+0*t62 923 t67 = 1/t7 924 t68 = 2.5d-1*t31*t67+t66-5.0d-1*t28*t35 925 t69 = 1.111111111111111d-1*t68*t7+2.222222222222222d-1*rhoa* 926 1 t36 927 t70 = -2.5d-1*t31*t67+t66+5.0d-1*t28*t35 928 t71 = 1.111111111111111d-1*t7*t70+2.222222222222222d-1*rhoa* 929 1 t38-2 930 t72 = -5.481209360000001d-4*gammaaa*t25*t27*t30*t37 931 t73 = -5.481209360000001d-4*gammaaa*t25*t27*t30*t39 932 t74 = -5.329547801878944d-2*t25*t30*t60 933 t75 = -1.9985804257046041d-2*t12*t30*t60 934 t76 = t75+t74+t73+t72-6.491760000000001d-3*gammaaa*t30*t71-6 935 1 .491760000000001d-3*gammaaa*t30*t69+t65-6.491760000000001 936 2 d-3*gammaaa*t30*t64-2.3144502890117796d+0*t30*t58 937 t77 = -3.6666666666666664d+0*t14*t40*t47*t48 938 t78 = rhoa**3.3333333333333337d+0 939 t79 = -1.6874680727699207d-3*t12*t30*t78-9.326708653288152d- 940 1 2*t30*t5-9.255935539253335d-5*gammaaa*t25*t27*t30*t43-4.6 941 2 279677696266674d-5*gammaaa*t25*t27*t30*t42-5.481209360000 942 3 001d-4*gammaaa*t30*t39-5.481209360000001d-4*gammaaa*t30*t 943 4 37-5.481209360000001d-4*gammaaa*t29*t30 944 t80 = 3.125d-2*t14*t41*t79 945 t81 = 7.328128227583548d-2*t12*t30*t78+4.050288005770614d+0* 946 1 t30*t5+4.0195535306666674d-3*gammaaa*t25*t27*t30*t43+2.00 947 2 97767653333337d-3*gammaaa*t25*t27*t30*t42+2.3803120000000 948 3 005d-2*gammaaa*t30*t39+2.3803120000000005d-2*gammaaa*t30* 949 4 t37+2.3803120000000005d-2*gammaaa*t29*t30 950 t82 = t14*t47*t48*t81 951 t83 = 1.2401570718501566d-2 952 t84 = 1/rhoa**6.333333333333333d+0 953 t85 = 2.3266666666666663d-1*t45*t50*t83*t84 954 t86 = 1.1633333333333332d-1*t16*t46*t83*t84 955 t87 = 1/rhoa**6 956 t88 = -7.8125d-2*t14*t46*t87 957 t89 = 1/rhoa**5.666666666666667d+0 958 t90 = -2.3333333333333334d+0*t14*t47*t49*t89 959 t91 = -7.5520808d-4*gammaaa*t29*t30 960 t92 = -7.5520808d-4*gammaaa*t30*t37 961 t93 = -7.5520808d-4*gammaaa*t30*t39 962 t94 = -6.376473555466666d-5*gammaaa*t25*t27*t30*t42 963 t95 = -1.2752947110933333d-4*gammaaa*t25*t27*t30*t43 964 t96 = -2.3250152285696893d-3*t12*t30*t78 965 t97 = -1.2850459218308585d-1*t30*t5 966 t98 = t97+t96+t95+t94+t93+t92+t91+1.1633333333333332d-1*t40 967 t99 = 3.125d-2*t16*t41*t98 968 t100 = 1.1633333333333332d-1*t49-5*t45 969 t101 = 1.5625d-2*t100*t16*t87 970 t102 = t63+t55+1.111111111111111d-1*t21-2.666666666666666d+0 971 t103 = -3.8888888888888884d-1*t62*t7+1.111111111111111d-1*rh 972 1 oa*t38+1.111111111111111d-1*rhoa*t36+1.111111111111111d-1 973 2 *t32 974 t104 = t75+t74+t73+t72+t65-1.735837716758835d+0*t30*t58-1.29 975 1 8352d-2*gammaaa*t103*t30-6.491760000000001d-3*gammaaa*t10 976 2 2*t30 977 t105 = -2.36002525d-5*t16*t30*t41*t43 978 t106 = -1.7128779250000004d-5*t14*t30*t41*t43 979 t107 = 2.3803120000000005d-2*t14*t30*t43*t47*t48 980 t108 = 1/t13**4 981 t109 = -5.80714011061111d-5*t108*t56 982 t110 = 3.7272188675555545d-2*t47*t50*t59 983 t111 = -1.5973795146666664d-2*t17*t50*t59 984 t112 = 4.960628287400625d-2 985 t113 = -3.0004900148148145d-1*t112*t16*t61 986 t114 = 2.288509333333333d-1*t16*t52*t61 987 t115 = -4.4999148607197886d-3*rhoa*t30*t57 988 t116 = -1.6874680727699207d-3*rhoa*t22*t30 989 t117 = 7.106063735838591d-2*t15*t26*t30 990 t118 = -2.333333333333333d+0*rhoa*t62 991 t119 = 1/t78 992 t120 = 2.3457970370370365d-1*t16*t17*t6-3.619259259259259d-1 993 1 *t119*t14*t52-2.626814814814815d-1*t119*t52-2.95198256944 994 2 44435d-3*t50/rhoa**4+3.296774133555554d-3*t108*t112/rhoa* 995 3 *4.333333333333333d+0 996 t121 = -7.777777777777777d-1*t120*t7 997 t122 = -4.6279677696266674d-5*gammaaa*t29*t30*t57*t59 998 t123 = 7.308279146666667d-4*gammaaa*t15*t29*t30*t61 999 t124 = -3.5d+0*t120 1000 t125 = -4.6279677696266674d-5*gammaaa*t30*t37*t57*t59 1001 t126 = 7.308279146666667d-4*gammaaa*t15*t30*t37*t61 1002 t127 = -4.6279677696266674d-5*gammaaa*t30*t39*t57*t59 1003 t128 = 7.308279146666667d-4*gammaaa*t15*t30*t39*t61 1004 t129 = -1.4247855427754028d-4*t22*t30*t7 1005 t130 = -9.255935539253335d-5*gammaaa*t25*t27*t29*t30 1006 t131 = -9.255935539253335d-5*gammaaa*t25*t27*t30*t37 1007 t132 = -9.255935539253335d-5*gammaaa*t25*t27*t30*t39 1008 t133 = -3.907547453488116d-6*gammaaa*t30*t42*t57*t59 1009 t134 = 6.170623692835556d-5*gammaaa*t15*t30*t42*t61 1010 t135 = -7.815094906976232d-6*gammaaa*t30*t43*t57*t59 1011 t136 = 1.2341247385671113d-4*gammaaa*t15*t30*t43*t61 1012 t137 = -6.749872291079682d-3*t25*t30*t60 1013 t138 = -3.3749361455398413d-3*t12*t30*t60 1014 t139 = 6.187382933489709d-3*t22*t30*t7 1015 t140 = 4.0195535306666674d-3*gammaaa*t25*t27*t29*t30 1016 t141 = 4.0195535306666674d-3*gammaaa*t25*t27*t30*t37 1017 t142 = 4.0195535306666674d-3*gammaaa*t25*t27*t30*t39 1018 t143 = 1.6969215155297782d-4*gammaaa*t30*t42*t57*t59 1019 t144 = -2.679702353777778d-3*gammaaa*t15*t30*t42*t61 1020 t145 = 3.3938430310595563d-4*gammaaa*t30*t43*t57*t59 1021 t146 = -5.359404707555556d-3*gammaaa*t15*t30*t43*t61 1022 t147 = 2.9312512910334193d-1*t25*t30*t60 1023 t148 = 1.4656256455167097d-1*t12*t30*t60 1024 t149 = 8.555555555555555d+0*t14*t40*t47*t89 1025 t150 = 2.3266666666666663d-1*t16*t79*t83*t84 1026 t151 = -1.5625d-1*t14*t79*t87 1027 t152 = -2.3333333333333334d+0*t14*t17*t81*t89 1028 t153 = 4.921566601151848d-3 1029 t154 = 1/rhoa**7.666666666666667d+0 1030 t155 = 8.120066666666664d-2*t108*t153*t154*t45 1031 t156 = 2.706688888888888d-2*t153*t154*t46*t50 1032 t157 = 6.200785359250782d-3 1033 t158 = 1/rhoa**7.333333333333333d+0 1034 t159 = -1.3184444444444443d+0*t157*t158*t16*t46 1035 t160 = 1/rhoa**7 1036 t161 = 2.34375d-1*t14*t160*t46 1037 t162 = 1.9686266404607394d-2 1038 t163 = 1.322222222222222d+1*t14*t162*t49/rhoa**6.66666666666 1039 1 6667d+0 1040 t164 = -1.9630878579890076d-4*t22*t30*t7 1041 t165 = -1.2752947110933333d-4*gammaaa*t25*t27*t29*t30 1042 t166 = -1.2752947110933333d-4*gammaaa*t25*t27*t30*t37 1043 t167 = -1.2752947110933333d-4*gammaaa*t25*t27*t30*t39 1044 t168 = -5.383869171999022d-6*gammaaa*t30*t42*t57*t59 1045 t169 = 8.501964740622221d-5*gammaaa*t15*t30*t42*t61 1046 t170 = -1.0767738343998043d-5*gammaaa*t30*t43*t57*t59 1047 t171 = 1.7003929481244442d-4*gammaaa*t15*t30*t43*t61 1048 t172 = -4.6500304571393786d-3*t12*t30*t60 1049 t173 = t97+t96+t95+t94+t93+t92+t91 1050 t174 = 1.5625d-2*t16*t87*(-5*t98+2.3266666666666663d-1*t81-4 1051 1 .2655555555555547d-1*t40-5*t173) 1052 t175 = t50*t83*t84*(2.3266666666666663d-1*t98+2.326666666666 1053 1 6663d-1*t173) 1054 t176 = 7.8125d-3*t16*t160*(-5.428888888888889d-1*t49-6*t100) 1055 t177 = t157*t158*(2.3266666666666663d-1*t100-1.4735555555555 1056 1 552d+0*t45)*t50 1057 t178 = -7.777777777777777d-1*rhoa*t62 1058 t179 = -1.757117466133333d-4*t30*t43*t50*t83*t84 1059 t180 = -6.376473555466666d-5*t16*t30*t43*t83*t84 1060 t181 = 1.0226776083333333d-4*t16*t30*t43*t87 1061 t182 = 4.282194812500001d-5*t14*t30*t43*t87 1062 t183 = -1.110812266666667d-1*t14*t162*t30*t43*t89 1063 t184 = -5.481209360000001d-4*t25*t27*t30*t37 1064 t185 = -6.376473555466666d-5*t25*t27*t30*t43 1065 t186 = -4.6279677696266674d-5*t25*t27*t30*t43 1066 t187 = 2.0097767653333337d-3*t25*t27*t30*t43 1067 t188 = -1.757117466133333d-4*t30*t42*t50*t83*t84 1068 t189 = -6.376473555466666d-5*t16*t30*t42*t83*t84 1069 t190 = 1.0226776083333333d-4*t16*t30*t42*t87 1070 t191 = 4.282194812500001d-5*t14*t30*t42*t87 1071 t192 = -1.110812266666667d-1*t14*t162*t30*t42*t89 1072 t193 = -5.481209360000001d-4*t25*t27*t29*t30 1073 t194 = 3.125d-2*t16*t41*(-6.376473555466666d-5*t25*t27*t30*t 1074 1 42-1.51041616d-3*t29*t30) 1075 t195 = 3.125d-2*t14*t41*(-4.6279677696266674d-5*t25*t27*t30* 1076 1 t42-5.481209360000001d-4*t29*t30) 1077 t196 = t14*(2.0097767653333337d-3*t25*t27*t30*t42+4.76062400 1078 1 0000001d-2*t29*t30)*t47*t48 1079 t197 = 3.125d-2*t14*(t186-5.481209360000001d-4*t30*t39)*t41 1080 fnc(iq) = 1.0d+0*(-1.0223881343581931d-3*gammaaa*t2*t4*t6*(1 1081 1 .111111111111111d-1*t7*(-5.0d-1*(t9+t3-1.1d+1)-3.0d+0*t10 1082 2 +1.0d+0)+t11)-5.111940671790965d-4*gammaaa*t2*t4*t6*(t8+1 1083 3 .111111111111111d-1*(4.7d+1-7.0d+0*t10)*t7)-3.72787240661 1084 4 23486d-2*rhoa*t2*t4-9.836d-2*rhoa*t2)*wght+fnc(iq) 1085 Amat(iq,D1_RA) = 1.0d+0*(t14*t17*t40*t6+t14*t47*t48*t49+3.12 1086 1 5d-2*t14*t41*t46+3.125d-2*t16*t41*t45-2.288509333333333d- 1087 2 2*t1*t15*t16-4.918d-2*t14)*wght+Amat(iq,D1_RA) 1088 Cmat(iq,D1_GAA) = Cmat(iq,D1_GAA)-6.491760000000001d-3*t14*t 1089 1 17*t30*t43*t6*wght 1090 Cmat(iq,D1_GAB) = Cmat(iq,D1_GAB)-6.491760000000001d-3*t14*t 1091 1 17*t30*t42*t6*wght 1092 Amat2(iq,D2_RA_RA) = 1.0d+0*(t99+t90+t88+t86+t85+t82+t80+t77 1093 1 +t14*t17*t6*t76+t54+t53+t51+4.918d-2*t14*t35+t101)*wght+A 1094 2 mat2(iq,D2_RA_RA) 1095 Amat2(iq,D2_RA_RB) = 1.0d+0*(t99+t90+t88+t86+t85+t82+t80+t77 1096 1 +t104*t14*t17*t6+t54+t53+t51-4.918d-2*t14*t35+t101)*wght+ 1097 2 Amat2(iq,D2_RA_RB) 1098 Cmat2(iq,D2_RA_GAA) = 1.0d+0*(-6.491760000000001d-3*t14*t17* 1099 1 t30*t37*t6+t107+t106+t105)*wght+Cmat2(iq,D2_RA_GAA) 1100 Cmat2(iq,D2_RA_GAB) = 1.0d+0*(-6.491760000000001d-3*t14*t17* 1101 1 t29*t30*t6+2.3803120000000005d-2*t14*t30*t42*t47*t48-2.36 1102 2 002525d-5*t16*t30*t41*t42-1.7128779250000004d-5*t14*t30*t 1103 3 41*t42)*wght+Cmat2(iq,D2_RA_GAB) 1104 Cmat2(iq,D2_RA_GBB) = 1.0d+0*(-6.491760000000001d-3*t14*t17* 1105 1 t30*t39*t6+t107+t106+t105)*wght+Cmat2(iq,D2_RA_GBB) 1106 Cmat2(iq,D2_GAA_GAA) = Cmat2(iq,D2_GAA_GAA) 1107 Cmat2(iq,D2_GAA_GAB) = Cmat2(iq,D2_GAA_GAB) 1108 Cmat2(iq,D2_GAA_GBB) = Cmat2(iq,D2_GAA_GBB) 1109 Cmat2(iq,D2_GAB_GAB) = Cmat2(iq,D2_GAB_GAB) 1110 Amat3(iq,D3_RA_RA_RA) = 1.0d+0*(-3.6666666666666664d+0*t14*t 1111 1 17*t48*t76+3.125d-2*t16*t41*(2.3266666666666663d-1*t76-7. 1112 2 5520808d-4*gammaaa*t30*t71-7.5520808d-4*gammaaa*t30*t69-7 1113 3 .5520808d-4*gammaaa*t30*t64-9.300060914278757d-3*t25*t30* 1114 4 t60-2.69247716955037d-1*t30*t58+t172+t171+t170+t169+t168+ 1115 5 t167+t166+t165+t164)+t14*t17*t6*(-1.0962418720000001d-3*g 1116 6 ammaaa*t25*t27*t30*t71-6.491760000000001d-3*gammaaa*t30*( 1117 7 3.333333333333333d-1*rhoa*t70+1.111111111111111d-1*(-7.5d 1118 8 -1*t28*t67+7.5d-1*t35*t62+3.75d-1*t31*t56+t124)*t7)-6.491 1119 9 760000000001d-3*gammaaa*t30*(1.111111111111111d-1*(7.5d-1 1120 : *t28*t67-7.5d-1*t35*t62-3.75d-1*t31*t56+t124)*t7+3.333333 1121 ; 333333333d-1*rhoa*t68)-1.0962418720000001d-3*gammaaa*t25* 1122 < t27*t30*t69-1.0962418720000001d-3*gammaaa*t25*t27*t30*t64 1123 = -3.6418576646172784d-1*t25*t26*t30-3.857417148352966d+0*t 1124 > 23*t30-6.491760000000001d-3*gammaaa*(t121+t118)*t30+t128+ 1125 ? t127+t126+t125+t123+t122+t117+t116+t115)+t14*t47*t48*(2.3 1126 @ 803120000000005d-2*gammaaa*t30*t71+2.3803120000000005d-2* 1127 1 gammaaa*t30*t69+2.3803120000000005d-2*gammaaa*t30*t64+8.4 1128 2 86317726376527d+0*t30*t58+t148+t147+t146+t145+t144+t143+t 1129 3 142+t141+t140+t139)+3.125d-2*t14*t41*(-5.481209360000001d 1130 4 -4*gammaaa*t30*t71-5.481209360000001d-4*gammaaa*t30*t69-5 1131 5 .481209360000001d-4*gammaaa*t30*t64-1.9541675273556125d-1 1132 6 *t30*t58+t138+t137+t136+t135+t134+t133+t132+t131+t130+t12 1133 7 9)-7.377d-2*t14*t67+t177+t176+t175+t174+t163+t161+t159+t1 1134 8 56+t155+t152+t151+t150+t149+t114+t113+t111+t110+t109)*wgh 1135 9 t+Amat3(iq,D3_RA_RA_RA) 1136 Amat3(iq,D3_RA_RA_RB) = 1.0d+0*(3.125d-2*t16*t41*(1.16333333 1137 1 33333332d-1*t76-9.300060914278756d-3*t25*t30*t60-2.019357 1138 2 8771627776d-1*t30*t58-1.51041616d-3*gammaaa*t103*t30-7.55 1139 3 20808d-4*gammaaa*t102*t30+t172+t171+t170+t169+t168+t167+t 1140 4 166+t165+t164+1.1633333333333332d-1*t104)+3.3333333333333 1141 5 33d-1*t14*t47*t48*(-11*t76-11*t104)+t14*t17*t6*(-5.481209 1142 6 360000001d-4*gammaaa*t25*t27*t30*t71-6.491760000000001d-3 1143 7 *gammaaa*t30*(1.111111111111111d-1*rhoa*t70+1.11111111111 1144 8 1111d-1*(-2.5d-1*t28*t67+2.5d-1*t35*t62+1.25d-1*t31*t56+t 1145 9 124)*t7+2.222222222222222d-1*t38+t178)-6.491760000000001d 1146 : -3*gammaaa*t30*(1.111111111111111d-1*(2.5d-1*t28*t67-2.5d 1147 ; -1*t35*t62-1.25d-1*t31*t56+t124)*t7+1.111111111111111d-1* 1148 < rhoa*t68+2.222222222222222d-1*t36+t178)-5.481209360000001 1149 = d-4*gammaaa*t25*t27*t30*t69-5.481209360000001d-4*gammaaa* 1150 > t25*t27*t30*t64-6.491760000000001d-3*gammaaa*(-1.55555555 1151 ? 55555553d+0*t28+t121+t118)*t30-1.0962418720000001d-3*gamm 1152 @ aaa*t103*t25*t27*t30-5.481209360000001d-4*gammaaa*t102*t2 1153 1 5*t27*t30-2.7535996976374544d-1*t25*t26*t30-1.99858042570 1154 2 46041d-2*t12*t26*t30-2.3144502890117796d+0*t23*t30+t128+t 1155 3 127+t126+t125+t123+t122+t117+t116+t115)+2.459d-2*t14*t67- 1156 4 4.577018666666666d-2*t15*t16*t61+t14*t47*t48*(6.364738294 1157 5 782395d+0*t30*t58+4.760624000000001d-2*gammaaa*t103*t30+2 1158 6 .3803120000000005d-2*gammaaa*t102*t30+t148+t147+t146+t145 1159 7 +t144+t143+t142+t141+t140+t139)+3.125d-2*t14*t41*(-1.4656 1160 8 256455167094d-1*t30*t58-1.0962418720000001d-3*gammaaa*t10 1161 9 3*t30-5.481209360000001d-4*gammaaa*t102*t30+t138+t137+t13 1162 : 6+t135+t134+t133+t132+t131+t130+t129)+t177+t176+t175+t174 1163 ; +t163+t161+t159+t156+t155+t152+t151+t150+t149+t114+t113+t 1164 < 111+t110+t109)*wght+Amat3(iq,D3_RA_RA_RB) 1165 Cmat3(iq,D3_RA_RA_GAA) = 1.0d+0*(t14*t17*t6*(t184-6.49176000 1166 1 0000001d-3*t30*t69)+t14*(4.760624000000001d-2*t30*t37+t18 1167 2 7)*t47*t48+3.125d-2*t14*(t186-5.481209360000001d-4*t30*t3 1168 3 7)*t41+3.125d-2*t16*(t185-1.51041616d-3*t30*t37)*t41+t183 1169 4 +t182+t181+t180+t179)*wght+Cmat3(iq,D3_RA_RA_GAA) 1170 Cmat3(iq,D3_RA_RA_GAB) = 1.0d+0*(t14*t17*t6*(t193-6.49176000 1171 1 0000001d-3*t30*t64)+t196+t195+t194+t192+t191+t190+t189+t1 1172 2 88)*wght+Cmat3(iq,D3_RA_RA_GAB) 1173 Cmat3(iq,D3_RA_RA_GBB) = 1.0d+0*(t14*t17*t6*(-6.491760000000 1174 1 001d-3*t30*t71-5.481209360000001d-4*t25*t27*t30*t39)+t14* 1175 2 (4.760624000000001d-2*t30*t39+t187)*t47*t48+3.125d-2*t16* 1176 3 (t185-1.51041616d-3*t30*t39)*t41+t197+t183+t182+t181+t180 1177 4 +t179)*wght+Cmat3(iq,D3_RA_RA_GBB) 1178 Cmat3(iq,D3_RA_RB_GAA) = 1.0d+0*(t14*t17*(t184-6.49176000000 1179 1 0001d-3*t103*t30)*t6+t14*(2.3803120000000005d-2*t30*t39+2 1180 2 .3803120000000005d-2*t30*t37+t187)*t47*t48+3.125d-2*t16*( 1181 3 -7.5520808d-4*t30*t39-7.5520808d-4*t30*t37+t185)*t41+t197 1182 4 +t183+t182+t181+t180+t179)*wght+Cmat3(iq,D3_RA_RB_GAA) 1183 Cmat3(iq,D3_RA_RB_GAB) = 1.0d+0*(t14*t17*(t193-6.49176000000 1184 1 0001d-3*t102*t30)*t6+t196+t195+t194+t192+t191+t190+t189+t 1185 2 188)*wght+Cmat3(iq,D3_RA_RB_GAB) 1186 Cmat3(iq,D3_RA_GAA_GAA) = Cmat3(iq,D3_RA_GAA_GAA) 1187 Cmat3(iq,D3_RA_GAA_GAB) = Cmat3(iq,D3_RA_GAA_GAB) 1188 Cmat3(iq,D3_RA_GAA_GBB) = Cmat3(iq,D3_RA_GAA_GBB) 1189 Cmat3(iq,D3_RA_GAB_GAB) = Cmat3(iq,D3_RA_GAB_GAB) 1190 Cmat3(iq,D3_RA_GAB_GBB) = Cmat3(iq,D3_RA_GAB_GBB) 1191 Cmat3(iq,D3_RA_GBB_GBB) = Cmat3(iq,D3_RA_GBB_GBB) 1192 Cmat3(iq,D3_GAA_GAA_GAA) = Cmat3(iq,D3_GAA_GAA_GAA) 1193 Cmat3(iq,D3_GAA_GAA_GAB) = Cmat3(iq,D3_GAA_GAA_GAB) 1194 Cmat3(iq,D3_GAA_GAA_GBB) = Cmat3(iq,D3_GAA_GAA_GBB) 1195 Cmat3(iq,D3_GAA_GAB_GAB) = Cmat3(iq,D3_GAA_GAB_GAB) 1196 Cmat3(iq,D3_GAA_GAB_GBB) = Cmat3(iq,D3_GAA_GAB_GBB) 1197 Cmat3(iq,D3_GAB_GAB_GAB) = Cmat3(iq,D3_GAB_GAB_GAB) 1198 endif ! rhoa.gt.tol_rho 1199 else ! ipol.eq.1 1200 rhoa = rho(iq,R_A) 1201 rhob = rho(iq,R_B) 1202 gammaaa = rgamma(iq,G_AA) 1203 gammaab = rgamma(iq,G_AB) 1204 gammabb = rgamma(iq,G_BB) 1205 if (rhoa.gt.tol_rho.and.rhob.gt.tol_rho) then 1206 t1 = rhob+rhoa 1207 t2 = 1/t1 1208 t3 = 1/t1**3.333333333333333d-1 1209 t4 = 3.49d-1*t3+1.0d+0 1210 t5 = 1/t4 1211 t6 = 1/t1**3.6666666666666664d+0 1212 t7 = rhoa**2.6666666666666666d+0 1213 t8 = rhob**2.6666666666666666d+0 1214 t9 = t8+t7 1215 t10 = 2.533d-1*t3 1216 t11 = exp(-t10) 1217 t12 = t1**2 1218 t13 = 3.49d-1*t3*t5 1219 t14 = t13+t10 1220 t15 = 4.7d+1-7.0d+0*t14 1221 t16 = 1.111111111111111d-1*rhoa*rhob*t15-1.333333333333333d+ 1222 1 0*t12 1223 t17 = t13+t10-1.1d+1 1224 t18 = -3.0d+0*t14 1225 t19 = -1.0d+0*rhoa*t17*t2+t18+1.0d+0 1226 t20 = 1.111111111111111d-1*rhoa*rhob*t19-rhob**2 1227 t21 = -1.0d+0*rhob*t17*t2+t18+1.0d+0 1228 t22 = 1.111111111111111d-1*rhoa*rhob*t21-rhoa**2 1229 t23 = 1/t1**2.3333333333333334d+0 1230 t24 = 1/t4**2 1231 t25 = -2.288509333333333d-2*rhoa*rhob*t23*t24 1232 t26 = 1/t12 1233 t27 = 1.9672d-1*rhoa*rhob*t26*t5 1234 t28 = -2.666666666666666d+0*t1 1235 t29 = 1/t1**1.3333333333333333d+0 1236 t30 = -1.1633333333333332d-1*t29*t5-8.443333333333334d-2*t29 1237 1 +4.060033333333332d-2*t24/t1**1.6666666666666669d+0 1238 t31 = -7.777777777777777d-1*rhoa*rhob*t30 1239 t32 = t31+t28+1.111111111111111d-1*rhob*t15 1240 t33 = -3.0d+0*t30 1241 t34 = -1.0d+0*rhoa*t2*t30 1242 t35 = 1.0d+0*rhoa*t17*t26 1243 t36 = -1.0d+0*t17*t2 1244 t37 = t36+t35+t34+t33 1245 t38 = 1.111111111111111d-1*rhoa*rhob*t37+1.111111111111111d- 1246 1 1*rhob*t19 1247 t39 = -1.0d+0*rhob*t2*t30 1248 t40 = 1.0d+0*rhob*t17*t26 1249 t41 = t40+t39+t33 1250 t42 = 1.111111111111111d-1*rhoa*rhob*t41+1.111111111111111d- 1251 1 1*rhob*t21-2*rhoa 1252 t43 = -2.367051431943866d-1*rhob*t11*t9-6.312137151850309d-1 1253 1 *rhob*t11*t7-6.491760000000001d-3*gammabb*t11*t42-6.49176 1254 2 0000000001d-3*gammaaa*t11*t38-6.491760000000001d-3*gammaa 1255 3 b*t11*t32 1256 t44 = 1/t1**5 1257 t45 = -2.7536698324946973d-2*rhoa*rhob*t11*t9-7.5520808d-4*g 1258 1 ammabb*t11*t22-7.5520808d-4*gammaaa*t11*t20-7.5520808d-4* 1259 2 gammaab*t11*t16 1260 t46 = t24*t44*t45 1261 t47 = -1.9985804257046041d-2*rhoa*rhob*t11*t9-5.481209360000 1262 1 001d-4*gammabb*t11*t22-5.481209360000001d-4*gammaaa*t11*t 1263 2 20-5.481209360000001d-4*gammaab*t11*t16 1264 t48 = t44*t47*t5 1265 t49 = 1/t1**4.666666666666667d+0 1266 t50 = 8.679188583794175d-1*rhoa*rhob*t11*t9+2.38031200000000 1267 1 05d-2*gammabb*t11*t22+2.3803120000000005d-2*gammaaa*t11*t 1268 2 20+2.3803120000000005d-2*gammaab*t11*t16 1269 t51 = t49*t5*t50 1270 t52 = t31+t28+1.111111111111111d-1*rhoa*t15 1271 t53 = t35+t34+t33 1272 t54 = 1.111111111111111d-1*rhoa*rhob*t53+1.111111111111111d- 1273 1 1*rhoa*t19-2*rhob 1274 t55 = t40+t39+t36+t33 1275 t56 = 1.111111111111111d-1*rhoa*rhob*t55+1.111111111111111d- 1276 1 1*rhoa*t21 1277 t57 = -2.367051431943866d-1*rhoa*t11*t9-6.312137151850309d-1 1278 1 *rhoa*t11*t8-6.491760000000001d-3*gammabb*t11*t56-6.49176 1279 2 0000000001d-3*gammaaa*t11*t54-6.491760000000001d-3*gammaa 1280 3 b*t11*t52 1281 t58 = 1/t4**3 1282 t59 = -5.324598382222221d-3*rhoa*rhob*t58*t6 1283 t60 = 1/t1**3.3333333333333337d+0 1284 t61 = 7.628364444444444d-2*rhoa*rhob*t24*t60 1285 t62 = 1/t1**3 1286 t63 = -3.9344d-1*rhoa*rhob*t5*t62 1287 t64 = -3.6666666666666664d+0*t43*t49*t5 1288 t65 = rhoa**1.6666666666666669d+0 1289 t66 = -1.9985804257046041d-2*rhob*t11*t29*t9 1290 t67 = -5.329547801878944d-2*rhob*t11*t29*t7 1291 t68 = 1/t1**2.6666666666666666d+0 1292 t69 = -8.120066666666665d-2*t24*t68+9.446344222222218d-3*t58 1293 1 *t62+1.551111111111111d-1*t23*t5+1.1257777777777778d-1*t2 1294 2 3 1295 t70 = -7.777777777777777d-1*rhoa*rhob*t69 1296 t71 = t70-1.5555555555555553d+0*rhob*t30-2.666666666666666d+ 1297 1 0 1298 t72 = -3.0d+0*t69 1299 t73 = -1.0d+0*rhob*t2*t69 1300 t74 = 2.0d+0*rhob*t26*t30 1301 t75 = -2.0d+0*rhob*t17*t62 1302 t76 = t75+t74+t73+t72 1303 t77 = 1.111111111111111d-1*rhoa*rhob*t76+2.222222222222222d- 1304 1 1*rhob*t41-2 1305 t78 = -1.0d+0*rhoa*t2*t69 1306 t79 = 2.0d+0*rhoa*t26*t30 1307 t80 = -2.0d+0*t2*t30 1308 t81 = -2.0d+0*rhoa*t17*t62 1309 t82 = 2.0d+0*t17*t26 1310 t83 = t82+t81+t80+t79+t78+t72 1311 t84 = 1.111111111111111d-1*rhoa*rhob*t83+2.222222222222222d- 1312 1 1*rhob*t37 1313 t85 = -5.481209360000001d-4*gammaab*t11*t29*t32 1314 t86 = -5.481209360000001d-4*gammaaa*t11*t29*t38 1315 t87 = -5.481209360000001d-4*gammabb*t11*t29*t42 1316 t88 = t87+t86+t85-6.491760000000001d-3*gammaaa*t11*t84-6.491 1317 1 760000000001d-3*gammabb*t11*t77-6.491760000000001d-3*gamm 1318 2 aab*t11*t71+t67+t66-2.3144502890117796d+0*rhob*t11*t65 1319 t89 = 1/t1**6.333333333333333d+0 1320 t90 = 2.3266666666666663d-1*t45*t58*t89 1321 t91 = 1.1633333333333332d-1*t24*t47*t89 1322 t92 = 1/t1**6 1323 t93 = -5*t47*t5*t92 1324 t94 = 1/t1**5.666666666666667d+0 1325 t95 = -4.666666666666667d+0*t5*t50*t94 1326 t96 = -1.6874680727699207d-3*rhoa*rhob*t11*t29*t9 1327 t97 = -4.6279677696266674d-5*gammaab*t11*t16*t29 1328 t98 = -4.6279677696266674d-5*gammaaa*t11*t20*t29 1329 t99 = -4.6279677696266674d-5*gammabb*t11*t22*t29 1330 t100 = t99+t98+t97+t96-1.9985804257046041d-2*rhob*t11*t9-5.3 1331 1 29547801878943d-2*rhob*t11*t7-5.481209360000001d-4*gammab 1332 2 b*t11*t42-5.481209360000001d-4*gammaaa*t11*t38-5.48120936 1333 3 0000001d-4*gammaab*t11*t32 1334 t101 = 7.328128227583548d-2*rhoa*rhob*t11*t29*t9 1335 t102 = 2.0097767653333337d-3*gammaab*t11*t16*t29 1336 t103 = 2.0097767653333337d-3*gammaaa*t11*t20*t29 1337 t104 = 2.0097767653333337d-3*gammabb*t11*t22*t29 1338 t105 = 8.679188583794175d-1*rhob*t11*t9+2.3144502890117796d+ 1339 1 0*rhob*t11*t7+2.3803120000000005d-2*gammabb*t11*t42+2.380 1340 2 3120000000005d-2*gammaaa*t11*t38+2.3803120000000005d-2*ga 1341 3 mmaab*t11*t32+t104+t103+t102+t101 1342 t106 = -7.343119553319191d-2*rhob*t11*t7 1343 t107 = -2.7536698324946973d-2*rhob*t11*t9 1344 t108 = -2.3250152285696893d-3*rhoa*rhob*t11*t29*t9 1345 t109 = -7.5520808d-4*gammaab*t11*t32 1346 t110 = -6.376473555466666d-5*gammaab*t11*t16*t29 1347 t111 = -7.5520808d-4*gammaaa*t11*t38 1348 t112 = -6.376473555466666d-5*gammaaa*t11*t20*t29 1349 t113 = -7.5520808d-4*gammabb*t11*t42 1350 t114 = -6.376473555466666d-5*gammabb*t11*t22*t29 1351 t115 = 1.1633333333333332d-1*t43 1352 t116 = t115+t114+t113+t112+t111+t110+t109+t108+t107+t106 1353 t117 = 1.1633333333333332d-1*t50-5*t45 1354 t118 = t117*t24*t92 1355 t119 = -2.288509333333333d-2*rhob-2.288509333333333d-2*rhoa 1356 t120 = 1.9672d-1*rhob+1.9672d-1*rhoa 1357 t121 = t70-7.777777777777777d-1*rhob*t30-7.777777777777777d- 1358 1 1*rhoa*t30+1.111111111111111d-1*t15-2.666666666666666d+0 1359 t122 = -1.0d+0*t2*t30 1360 t123 = 1.0d+0*t17*t26 1361 t124 = t81+t79+t78+t72+t123+t122 1362 t125 = 1.111111111111111d-1*rhob*t53+1.111111111111111d-1*rh 1363 1 oa*t37+1.111111111111111d-1*t19+1.111111111111111d-1*rhoa 1364 2 *rhob*t124 1365 t126 = t75+t74+t73+t72+t123+t122 1366 t127 = 1.111111111111111d-1*rhob*t55+1.111111111111111d-1*rh 1367 1 oa*t41+1.111111111111111d-1*t21+1.111111111111111d-1*rhoa 1368 2 *rhob*t126 1369 t128 = -2.367051431943866d-1*t11*t9+t87+t86+t85-6.3121371518 1370 1 50309d-1*t11*t8-6.312137151850309d-1*t11*t7+t67+t66-6.491 1371 2 760000000001d-3*gammabb*t11*t127-6.491760000000001d-3*gam 1372 3 maaa*t11*t125-6.491760000000001d-3*gammaab*t11*t121 1373 t129 = t99+t98+t97+t96-1.9985804257046041d-2*rhoa*t11*t9-5.3 1374 1 29547801878943d-2*rhoa*t11*t8-5.481209360000001d-4*gammab 1375 2 b*t11*t56-5.481209360000001d-4*gammaaa*t11*t54-5.48120936 1376 3 0000001d-4*gammaab*t11*t52 1377 t130 = t129*t44*t5 1378 t131 = 8.679188583794175d-1*rhoa*t11*t9+2.3144502890117796d+ 1379 1 0*rhoa*t11*t8+2.3803120000000005d-2*gammabb*t11*t56+2.380 1380 2 3120000000005d-2*gammaaa*t11*t54+2.3803120000000005d-2*ga 1381 3 mmaab*t11*t52+t104+t103+t102+t101 1382 t132 = t131*t49*t5 1383 t133 = -7.343119553319191d-2*rhoa*t11*t8 1384 t134 = -2.7536698324946973d-2*rhoa*t11*t9 1385 t135 = -7.5520808d-4*gammaab*t11*t52 1386 t136 = -7.5520808d-4*gammaaa*t11*t54 1387 t137 = -7.5520808d-4*gammabb*t11*t56 1388 t138 = t137+t136+t135+t134+t133+t115+t114+t112+t110+t108 1389 t139 = rhob**1.6666666666666669d+0 1390 t140 = t70-1.5555555555555553d+0*rhoa*t30-2.666666666666666d 1391 1 +0 1392 t141 = t81+t79+t78+t72 1393 t142 = 2.222222222222222d-1*rhoa*t53+1.111111111111111d-1*rh 1394 1 oa*rhob*t141-2 1395 t143 = t82+t80+t75+t74+t73+t72 1396 t144 = 2.222222222222222d-1*rhoa*t55+1.111111111111111d-1*rh 1397 1 oa*rhob*t143 1398 t145 = -1.9985804257046041d-2*rhoa*t11*t29*t9-5.329547801878 1399 1 944d-2*rhoa*t11*t29*t8-5.481209360000001d-4*gammabb*t11*t 1400 2 29*t56-5.481209360000001d-4*gammaaa*t11*t29*t54-5.4812093 1401 3 60000001d-4*gammaab*t11*t29*t52-6.491760000000001d-3*gamm 1402 4 abb*t11*t144-6.491760000000001d-3*gammaaa*t11*t142-6.4917 1403 5 60000000001d-3*gammaab*t11*t140-2.3144502890117796d+0*rho 1404 6 a*t11*t139 1405 t146 = 1.1633333333333332d-1*t57+t137+t136+t135+t134+t133+t1 1406 1 14+t112+t110+t108 1407 t147 = -7.5520808d-4*t11*t20*t24*t44 1408 t148 = -5.481209360000001d-4*t11*t20*t44*t5 1409 t149 = 2.3803120000000005d-2*t11*t20*t49*t5 1410 t150 = -7.5520808d-4*t11*t16*t24*t44 1411 t151 = -5.481209360000001d-4*t11*t16*t44*t5 1412 t152 = 2.3803120000000005d-2*t11*t16*t49*t5 1413 t153 = -7.5520808d-4*t11*t22*t24*t44 1414 t154 = -5.481209360000001d-4*t11*t22*t44*t5 1415 t155 = 2.3803120000000005d-2*t11*t22*t49*t5 1416 t156 = 1/t4**4 1417 t157 = -1.858284835395555d-3*rhoa*rhob*t156*t44 1418 t158 = 3.7272188675555545d-2*rhoa*rhob*t49*t58 1419 t159 = 1/t1**4.333333333333333d+0 1420 t160 = -3.0004900148148145d-1*rhoa*rhob*t159*t24 1421 t161 = 1/t1**4 1422 t162 = 1.18032d+0*rhoa*rhob*t161*t5 1423 t163 = 1.711111111111111d+1*t43*t5*t94 1424 t164 = -1.6874680727699207d-3*rhob*t11*t68*t9 1425 t165 = 2.6647739009394716d-2*rhob*t11*t23*t9 1426 t166 = -4.4999148607197886d-3*rhob*t11*t68*t7 1427 t167 = 7.106063735838591d-2*rhob*t11*t23*t7 1428 t168 = -3.619259259259259d-1*t5*t60-2.626814814814815d-1*t60 1429 1 +2.3457970370370365d-1*t24*t6-4.7231721111111097d-2*t161* 1430 2 t58+3.296774133555554d-3*t156*t159 1431 t169 = -7.777777777777777d-1*rhoa*rhob*t168 1432 t170 = -3.0d+0*t168 1433 t171 = -1.0d+0*rhob*t168*t2 1434 t172 = 3.0d+0*rhob*t26*t69 1435 t173 = -6.0d+0*rhob*t30*t62 1436 t174 = 6.0d+0*rhob*t161*t17 1437 t175 = -1.0d+0*rhoa*t168*t2 1438 t176 = 3.0d+0*rhoa*t26*t69 1439 t177 = -3.0d+0*t2*t69 1440 t178 = -6.0d+0*rhoa*t30*t62 1441 t179 = 6.0d+0*t26*t30 1442 t180 = 6.0d+0*rhoa*t161*t17 1443 t181 = -6.0d+0*t17*t62 1444 t182 = -4.6279677696266674d-5*gammaab*t11*t32*t68 1445 t183 = 7.308279146666667d-4*gammaab*t11*t23*t32 1446 t184 = -4.6279677696266674d-5*gammaaa*t11*t38*t68 1447 t185 = 7.308279146666667d-4*gammaaa*t11*t23*t38 1448 t186 = -4.6279677696266674d-5*gammabb*t11*t42*t68 1449 t187 = 7.308279146666667d-4*gammabb*t11*t23*t42 1450 t188 = 1/t1**7.666666666666667d+0 1451 t189 = 8.120066666666664d-2*t156*t188*t45 1452 t190 = 2.706688888888888d-2*t188*t47*t58 1453 t191 = 1/t1**7.333333333333333d+0 1454 t192 = -1.3184444444444443d+0*t191*t24*t47 1455 t193 = 1/t1**7 1456 t194 = 30*t193*t47*t5 1457 t195 = 2.644444444444444d+1*t5*t50/t1**6.666666666666667d+0 1458 t196 = 6.187382933489709d-3*rhoa*rhob*t11*t68*t9 1459 t197 = -9.770837636778064d-2*rhoa*rhob*t11*t23*t9 1460 t198 = 1.6969215155297782d-4*gammaab*t11*t16*t68 1461 t199 = -2.679702353777778d-3*gammaab*t11*t16*t23 1462 t200 = 1.6969215155297782d-4*gammaaa*t11*t20*t68 1463 t201 = -2.679702353777778d-3*gammaaa*t11*t20*t23 1464 t202 = 1.6969215155297782d-4*gammabb*t11*t22*t68 1465 t203 = -2.679702353777778d-3*gammabb*t11*t22*t23 1466 t204 = -1.4247855427754028d-4*rhoa*rhob*t11*t68*t9 1467 t205 = 2.249957430359894d-3*rhoa*rhob*t11*t23*t9 1468 t206 = -3.907547453488116d-6*gammaab*t11*t16*t68 1469 t207 = 6.170623692835556d-5*gammaab*t11*t16*t23 1470 t208 = -3.907547453488116d-6*gammaaa*t11*t20*t68 1471 t209 = 6.170623692835556d-5*gammaaa*t11*t20*t23 1472 t210 = -3.907547453488116d-6*gammabb*t11*t22*t68 1473 t211 = 6.170623692835556d-5*gammabb*t11*t22*t23 1474 t212 = -1.9630878579890076d-4*rhoa*rhob*t11*t68*t9 1475 t213 = 3.100020304759586d-3*rhoa*rhob*t11*t23*t9 1476 t214 = -5.383869171999022d-6*gammaab*t11*t16*t68 1477 t215 = 8.501964740622221d-5*gammaab*t11*t16*t23 1478 t216 = -5.383869171999022d-6*gammaaa*t11*t20*t68 1479 t217 = 8.501964740622221d-5*gammaaa*t11*t20*t23 1480 t218 = -5.383869171999022d-6*gammabb*t11*t22*t68 1481 t219 = 8.501964740622221d-5*gammabb*t11*t22*t23 1482 t220 = -4.2655555555555547d-1*t43 1483 t221 = t114+t113+t112+t111+t110+t109+t108+t107+t106 1484 t222 = -5*t116 1485 t223 = 2.3266666666666663d-1*t116 1486 t224 = t193*t24*(-5.428888888888889d-1*t50-6*t117) 1487 t225 = t191*(2.3266666666666663d-1*t117-1.4735555555555552d+ 1488 1 0*t45)*t58 1489 t226 = -5.324598382222221d-3*rhoa 1490 t227 = 7.628364444444444d-2*rhoa 1491 t228 = -4.577018666666666d-2*t23*t24 1492 t229 = -3.9344d-1*rhoa 1493 t230 = 3.9344d-1*t26*t5 1494 t231 = -1.5555555555555553d+0*t30 1495 t232 = -1.0d+0*t2*t69 1496 t233 = 2.0d+0*t26*t30 1497 t234 = -2.0d+0*t17*t62 1498 t235 = -2.0d+0*t2*t69 1499 t236 = 4.0d+0*t26*t30 1500 t237 = -4.0d+0*t17*t62 1501 t238 = t137+t136+t135+t134+t133+t114+t112+t110+t108 1502 t239 = -5*t238 1503 t240 = 2.3266666666666663d-1*t238 1504 t241 = t49*t5*(1.4656256455167097d-1*rhoa*t11*t29*t9+3.90833 1505 1 50547112256d-1*rhoa*t11*t29*t8+4.0195535306666674d-3*gamm 1506 2 abb*t11*t29*t56+4.0195535306666674d-3*gammaaa*t11*t29*t54 1507 3 +4.0195535306666674d-3*gammaab*t11*t29*t52+t203+t202+t201 1508 4 +t200+t199+t198+t197+t196+2.3803120000000005d-2*gammabb*t 1509 5 11*t144+2.3803120000000005d-2*gammaaa*t11*t142+2.38031200 1510 6 00000005d-2*gammaab*t11*t140+8.486317726376527d+0*rhoa*t1 1511 7 1*t139) 1512 t242 = t44*t5*(-3.3749361455398413d-3*rhoa*t11*t29*t9-8.9998 1513 1 29721439576d-3*rhoa*t11*t29*t8-9.255935539253335d-5*gamma 1514 2 bb*t11*t29*t56-9.255935539253335d-5*gammaaa*t11*t29*t54-9 1515 3 .255935539253335d-5*gammaab*t11*t29*t52+t211+t210+t209+t2 1516 4 08+t207+t206+t205+t204-5.481209360000001d-4*gammabb*t11*t 1517 5 144-5.481209360000001d-4*gammaaa*t11*t142-5.4812093600000 1518 6 01d-4*gammaab*t11*t140-1.9541675273556125d-1*rhoa*t11*t13 1519 7 9) 1520 t243 = 2.3266666666666663d-1*t129*t24*t89 1521 t244 = -10*t129*t5*t92 1522 t245 = -9.333333333333333d+0*t131*t5*t94 1523 t246 = -2.69247716955037d-1*rhoa*t11*t139 1524 t247 = -4.6500304571393786d-3*rhoa*t11*t29*t9 1525 t248 = -1.2400081219038343d-2*rhoa*t11*t29*t8 1526 t249 = -7.5520808d-4*gammaab*t11*t140 1527 t250 = -7.5520808d-4*gammaaa*t11*t142 1528 t251 = -7.5520808d-4*gammabb*t11*t144 1529 t252 = -1.2752947110933333d-4*gammaab*t11*t29*t52 1530 t253 = -1.2752947110933333d-4*gammaaa*t11*t29*t54 1531 t254 = -1.2752947110933333d-4*gammabb*t11*t29*t56 1532 t255 = 2.3266666666666663d-1*t131 1533 t256 = -1.757117466133333d-4*t11*t20*t58*t89 1534 t257 = -6.376473555466666d-5*t11*t20*t24*t89 1535 t258 = 6.545136693333333d-3*t11*t20*t24*t92 1536 t259 = 2.7406046800000006d-3*t11*t20*t5*t92 1537 t260 = -1.110812266666667d-1*t11*t20*t5*t94 1538 t261 = -5.481209360000001d-4*t11*t29*t38 1539 t262 = -6.376473555466666d-5*t11*t20*t29 1540 t263 = -4.6279677696266674d-5*t11*t20*t29 1541 t264 = 2.0097767653333337d-3*t11*t20*t29 1542 t265 = -1.757117466133333d-4*t11*t16*t58*t89 1543 t266 = -6.376473555466666d-5*t11*t16*t24*t89 1544 t267 = 6.545136693333333d-3*t11*t16*t24*t92 1545 t268 = 2.7406046800000006d-3*t11*t16*t5*t92 1546 t269 = -1.110812266666667d-1*t11*t16*t5*t94 1547 t270 = -5.481209360000001d-4*t11*t29*t32 1548 t271 = -6.376473555466666d-5*t11*t16*t29 1549 t272 = -4.6279677696266674d-5*t11*t16*t29 1550 t273 = 2.0097767653333337d-3*t11*t16*t29 1551 t274 = -1.757117466133333d-4*t11*t22*t58*t89 1552 t275 = -6.376473555466666d-5*t11*t22*t24*t89 1553 t276 = 6.545136693333333d-3*t11*t22*t24*t92 1554 t277 = 2.7406046800000006d-3*t11*t22*t5*t92 1555 t278 = -1.110812266666667d-1*t11*t22*t5*t94 1556 t279 = -5.481209360000001d-4*t11*t29*t42 1557 t280 = -6.376473555466666d-5*t11*t22*t29 1558 t281 = -4.6279677696266674d-5*t11*t22*t29 1559 t282 = 2.0097767653333337d-3*t11*t22*t29 1560 t283 = t44*t5*(t263-5.481209360000001d-4*t11*t54) 1561 t284 = t44*t5*(t272-5.481209360000001d-4*t11*t52) 1562 t285 = t44*t5*(t281-5.481209360000001d-4*t11*t56) 1563 fnc(iq) = 1.0d+0*(-2.367051431943866d-1*rhoa*rhob*t11*t5*t6* 1564 1 t9-6.491760000000001d-3*gammabb*t11*t22*t5*t6-6.491760000 1565 2 000001d-3*gammaaa*t11*t20*t5*t6-6.491760000000001d-3*gamm 1566 3 aab*t11*t16*t5*t6-1.9672d-1*rhoa*rhob*t2*t5)*wght+fnc(iq) 1567 Amat(iq,D1_RA) = 1.0d+0*(t43*t5*t6+t51-1.9672d-1*rhob*t2*t5+ 1568 1 t48+t46+t27+t25)*wght+Amat(iq,D1_RA) 1569 Amat(iq,D1_RB) = 1.0d+0*(t5*t57*t6+t51-1.9672d-1*rhoa*t2*t5+ 1570 1 t48+t46+t27+t25)*wght+Amat(iq,D1_RB) 1571 Cmat(iq,D1_GAA) = Cmat(iq,D1_GAA)-6.491760000000001d-3*t11*t 1572 1 20*t5*t6*wght 1573 Cmat(iq,D1_GAB) = Cmat(iq,D1_GAB)-6.491760000000001d-3*t11*t 1574 1 16*t5*t6*wght 1575 Cmat(iq,D1_GBB) = Cmat(iq,D1_GBB)-6.491760000000001d-3*t11*t 1576 1 22*t5*t6*wght 1577 Amat2(iq,D2_RA_RA) = 1.0d+0*(t95+t93+t91+t90+t5*t6*t88+t64+t 1578 1 63+t61+t59+t105*t49*t5+t100*t44*t5+3.9344d-1*rhob*t26*t5+ 1579 2 t116*t24*t44-4.577018666666666d-2*rhob*t23*t24+t118)*wght 1580 3 +Amat2(iq,D2_RA_RA) 1581 Amat2(iq,D2_RA_RB) = 1.0d+0*(t95+t93+t91+t90+t64+t63+t61+t12 1582 1 8*t5*t6+t59+t120*t26*t5-1.9672d-1*t2*t5+t138*t24*t44+t119 1583 2 *t23*t24+t132+t130+t118)*wght+Amat2(iq,D2_RA_RB) 1584 Amat2(iq,D2_RB_RB) = 1.0d+0*(t95+t93+t91+t90+t63+t61+t145*t5 1585 1 *t6+t59-3.6666666666666664d+0*t49*t5*t57+3.9344d-1*rhoa*t 1586 2 26*t5+t146*t24*t44-4.577018666666666d-2*rhoa*t23*t24+t132 1587 3 +t130+t118)*wght+Amat2(iq,D2_RB_RB) 1588 Cmat2(iq,D2_RA_GAA) = 1.0d+0*(-6.491760000000001d-3*t11*t38* 1589 1 t5*t6+t149+t148+t147)*wght+Cmat2(iq,D2_RA_GAA) 1590 Cmat2(iq,D2_RA_GAB) = 1.0d+0*(-6.491760000000001d-3*t11*t32* 1591 1 t5*t6+t152+t151+t150)*wght+Cmat2(iq,D2_RA_GAB) 1592 Cmat2(iq,D2_RA_GBB) = 1.0d+0*(-6.491760000000001d-3*t11*t42* 1593 1 t5*t6+t155+t154+t153)*wght+Cmat2(iq,D2_RA_GBB) 1594 Cmat2(iq,D2_RB_GAA) = 1.0d+0*(-6.491760000000001d-3*t11*t5*t 1595 1 54*t6+t149+t148+t147)*wght+Cmat2(iq,D2_RB_GAA) 1596 Cmat2(iq,D2_RB_GAB) = 1.0d+0*(-6.491760000000001d-3*t11*t5*t 1597 1 52*t6+t152+t151+t150)*wght+Cmat2(iq,D2_RB_GAB) 1598 Cmat2(iq,D2_RB_GBB) = 1.0d+0*(-6.491760000000001d-3*t11*t5*t 1599 1 56*t6+t155+t154+t153)*wght+Cmat2(iq,D2_RB_GBB) 1600 Cmat2(iq,D2_GAA_GAA) = Cmat2(iq,D2_GAA_GAA) 1601 Cmat2(iq,D2_GAA_GAB) = Cmat2(iq,D2_GAA_GAB) 1602 Cmat2(iq,D2_GAA_GBB) = Cmat2(iq,D2_GAA_GBB) 1603 Cmat2(iq,D2_GAB_GAB) = Cmat2(iq,D2_GAB_GAB) 1604 Cmat2(iq,D2_GAB_GBB) = Cmat2(iq,D2_GAB_GBB) 1605 Cmat2(iq,D2_GBB_GBB) = Cmat2(iq,D2_GBB_GBB) 1606 Amat3(iq,D3_RA_RA_RA) = 1.0d+0*(-9.333333333333333d+0*t105*t 1607 1 5*t94-10*t100*t5*t92+(t222-5*t221+t220+2.3266666666666663 1608 2 d-1*t105)*t24*t92+t49*t5*(1.4656256455167097d-1*rhob*t11* 1609 3 t29*t9+2.3803120000000005d-2*gammaaa*t11*t84+2.3803120000 1610 4 000005d-2*gammabb*t11*t77+2.3803120000000005d-2*gammaab*t 1611 5 11*t71+3.9083350547112256d-1*rhob*t11*t29*t7+8.4863177263 1612 6 76527d+0*rhob*t11*t65+4.0195535306666674d-3*gammabb*t11*t 1613 7 29*t42+4.0195535306666674d-3*gammaaa*t11*t29*t38+4.019553 1614 8 5306666674d-3*gammaab*t11*t29*t32+t203+t202+t201+t200+t19 1615 9 9+t198+t197+t196)+t44*t5*(-3.3749361455398413d-3*rhob*t11 1616 : *t29*t9-5.481209360000001d-4*gammaaa*t11*t84-5.4812093600 1617 ; 00001d-4*gammabb*t11*t77-5.481209360000001d-4*gammaab*t11 1618 < *t71-8.999829721439576d-3*rhob*t11*t29*t7-1.9541675273556 1619 = 125d-1*rhob*t11*t65-9.255935539253335d-5*gammabb*t11*t29* 1620 > t42-9.255935539253335d-5*gammaaa*t11*t29*t38-9.2559355392 1621 ? 53335d-5*gammaab*t11*t29*t32+t211+t210+t209+t208+t207+t20 1622 @ 6+t205+t204)+t24*t44*(-4.6500304571393786d-3*rhob*t11*t29 1623 1 *t9+2.3266666666666663d-1*t88-7.5520808d-4*gammaaa*t11*t8 1624 2 4-7.5520808d-4*gammabb*t11*t77-7.5520808d-4*gammaab*t11*t 1625 3 71-1.2400081219038343d-2*rhob*t11*t29*t7-2.69247716955037 1626 4 d-1*rhob*t11*t65-1.2752947110933333d-4*gammabb*t11*t29*t4 1627 5 2-1.2752947110933333d-4*gammaaa*t11*t29*t38-1.27529471109 1628 6 33333d-4*gammaab*t11*t29*t32+t219+t218+t217+t216+t215+t21 1629 7 4+t213+t212)+(t223+2.3266666666666663d-1*t221)*t58*t89+2. 1630 8 3266666666666663d-1*t100*t24*t89-7.333333333333333d+0*t49 1631 9 *t5*t88+t5*t6*(-1.0962418720000001d-3*gammaaa*t11*t29*t84 1632 : -6.491760000000001d-3*gammaaa*t11*(3.333333333333333d-1*r 1633 ; hob*t83+1.111111111111111d-1*rhoa*rhob*(t181+t180+t179+t1 1634 < 78+t177+t176+t175+t170))-1.0962418720000001d-3*gammabb*t1 1635 = 1*t29*t77-6.491760000000001d-3*gammabb*t11*(3.33333333333 1636 > 3333d-1*rhob*t76+1.111111111111111d-1*rhoa*rhob*(t174+t17 1637 ? 3+t172+t171+t170))-1.0962418720000001d-3*gammaab*t11*t29* 1638 @ t71-6.491760000000001d-3*gammaab*t11*(t169-2.333333333333 1639 1 333d+0*rhob*t69)-3.9083350547112256d-1*rhob*t11*t29*t65+t 1640 2 187+t186+t185+t184+t183+t182+t167+t166+t165+t164-3.857417 1641 3 148352966d+0*rhoa**6.666666666666666d-1*rhob*t11)-1.18032 1642 4 d+0*rhob*t5*t62+2.288509333333333d-1*rhob*t24*t60-1.59737 1643 5 95146666664d-2*rhob*t58*t6+t225+t224+t195+t194+t192+t190+ 1644 6 t189+t163+t162+t160+t158+t157)*wght+Amat3(iq,D3_RA_RA_RA) 1645 Amat3(iq,D3_RA_RA_RB) = 1.0d+0*(3.333333333333333d-1*(-14*t1 1646 1 31-14*t105)*t5*t94+(-5*t129-5*t100)*t5*t92+(t239+t222+t22 1647 2 0+1.1633333333333332d-1*t131+1.1633333333333332d-1*t105)* 1648 3 t24*t92+t49*t5*(7.328128227583548d-2*rhob*t11*t29*t9+7.32 1649 4 8128227583548d-2*rhoa*t11*t29*t9+8.679188583794175d-1*t11 1650 5 *t9+1.9541675273556128d-1*rhoa*t11*t29*t8+2.3144502890117 1651 6 796d+0*t11*t8+1.9541675273556128d-1*rhob*t11*t29*t7+2.314 1652 7 4502890117796d+0*t11*t7+2.0097767653333337d-3*gammabb*t11 1653 8 *t29*t56+2.0097767653333337d-3*gammaaa*t11*t29*t54+2.0097 1654 9 767653333337d-3*gammaab*t11*t29*t52+2.0097767653333337d-3 1655 : *gammabb*t11*t29*t42+2.0097767653333337d-3*gammaaa*t11*t2 1656 ; 9*t38+2.0097767653333337d-3*gammaab*t11*t29*t32+t203+t202 1657 < +t201+t200+t199+t198+t197+t196+2.3803120000000005d-2*gamm 1658 = abb*t11*t127+2.3803120000000005d-2*gammaaa*t11*t125+2.380 1659 > 3120000000005d-2*gammaab*t11*t121)+t44*t5*(-1.68746807276 1660 ? 99207d-3*rhob*t11*t29*t9-1.6874680727699207d-3*rhoa*t11*t 1661 @ 29*t9-1.9985804257046041d-2*t11*t9-4.499914860719788d-3*r 1662 1 hoa*t11*t29*t8-5.329547801878943d-2*t11*t8-4.499914860719 1663 2 788d-3*rhob*t11*t29*t7-5.329547801878943d-2*t11*t7-4.6279 1664 3 677696266674d-5*gammabb*t11*t29*t56-4.6279677696266674d-5 1665 4 *gammaaa*t11*t29*t54-4.6279677696266674d-5*gammaab*t11*t2 1666 5 9*t52-4.6279677696266674d-5*gammabb*t11*t29*t42-4.6279677 1667 6 696266674d-5*gammaaa*t11*t29*t38-4.6279677696266674d-5*ga 1668 7 mmaab*t11*t29*t32+t211+t210+t209+t208+t207+t206+t205+t204 1669 8 -5.481209360000001d-4*gammabb*t11*t127-5.481209360000001d 1670 9 -4*gammaaa*t11*t125-5.481209360000001d-4*gammaab*t11*t121 1671 : )+t24*t44*(-2.3250152285696893d-3*rhob*t11*t29*t9-2.32501 1672 ; 52285696893d-3*rhoa*t11*t29*t9-2.7536698324946973d-2*t11* 1673 < t9+1.1633333333333332d-1*t88-6.200040609519172d-3*rhoa*t1 1674 = 1*t29*t8-7.343119553319191d-2*t11*t8-6.200040609519171d-3 1675 > *rhob*t11*t29*t7-7.343119553319191d-2*t11*t7-6.3764735554 1676 ? 66666d-5*gammabb*t11*t29*t56-6.376473555466666d-5*gammaaa 1677 @ *t11*t29*t54-6.376473555466666d-5*gammaab*t11*t29*t52-6.3 1678 1 76473555466666d-5*gammabb*t11*t29*t42-6.376473555466666d- 1679 2 5*gammaaa*t11*t29*t38-6.376473555466666d-5*gammaab*t11*t2 1680 3 9*t32+t219+t218+t217+t216+t215+t214+t213+t212+1.163333333 1681 4 3333332d-1*t128-7.5520808d-4*gammabb*t11*t127-7.5520808d- 1682 5 4*gammaaa*t11*t125-7.5520808d-4*gammaab*t11*t121)+t5*t6*( 1683 6 -1.9985804257046041d-2*t11*t29*t9-5.481209360000001d-4*ga 1684 7 mmaaa*t11*t29*t84-6.491760000000001d-3*gammaaa*t11*(1.111 1685 8 111111111111d-1*rhoa*t83+2.222222222222222d-1*t37+1.11111 1686 9 1111111111d-1*rhoa*rhob*(t237+t236+t235+t180+t178+t176+t1 1687 : 75+t170)+2.222222222222222d-1*rhob*t124)-5.32954780187894 1688 ; 3d-2*t11*t29*t8-5.481209360000001d-4*gammabb*t11*t29*t77- 1689 < 6.491760000000001d-3*gammabb*t11*(1.111111111111111d-1*rh 1690 = oa*t76+2.222222222222222d-1*t41+1.111111111111111d-1*rhoa 1691 > *rhob*(t234+t233+t232+t174+t173+t172+t171+t170)+2.2222222 1692 ? 22222222d-1*rhob*t126)-5.481209360000001d-4*gammaab*t11*t 1693 @ 29*t71-5.329547801878944d-2*t11*t29*t7-6.491760000000001d 1694 1 -3*gammaab*t11*(-1.5555555555555553d+0*rhob*t69-7.7777777 1695 2 77777777d-1*rhoa*t69+t231+t169)-1.9541675273556128d-1*rho 1696 3 b*t11*t29*t65-2.3144502890117796d+0*t11*t65-5.48120936000 1697 4 0001d-4*gammabb*t11*t127*t29-5.481209360000001d-4*gammaaa 1698 5 *t11*t125*t29-5.481209360000001d-4*gammaab*t11*t121*t29+t 1699 6 187+t186+t185+t184+t183+t182+t167+t166+t165+t164)+(t240+t 1700 7 223)*t58*t89+(1.1633333333333332d-1*t129+1.16333333333333 1701 8 32d-1*t100)*t24*t89+3.333333333333333d-1*t49*t5*(-11*t88- 1702 9 11*t128)+(t229-7.8688d-1*rhob)*t5*t62+(t227+1.52567288888 1703 : 88887d-1*rhob)*t24*t60+(t226-1.0649196764444441d-2*rhob)* 1704 ; t58*t6+t230+t228+t225+t224+t195+t194+t192+t190+t189+t163+ 1705 < t162+t160+t158+t157)*wght+Amat3(iq,D3_RA_RA_RB) 1706 Amat3(iq,D3_RA_RB_RB) = 1.0d+0*(t24*(t255+t239+t220-5*t138)* 1707 1 t92+t5*t6*(-3.9971608514092083d-2*t11*t29*t9-1.0659095603 1708 2 757887d-1*t11*t29*t8-1.0659095603757889d-1*t11*t29*t7-6.4 1709 3 91760000000001d-3*gammaab*t11*(-7.777777777777777d-1*rhob 1710 4 *t69-1.5555555555555553d+0*rhoa*t69+t231+t169)-6.49176000 1711 5 0000001d-3*gammabb*t11*(2.222222222222222d-1*t55+1.111111 1712 6 111111111d-1*rhoa*rhob*(t237+t236+t235+t174+t173+t172+t17 1713 7 1+t170)+1.111111111111111d-1*rhob*t143+2.222222222222222d 1714 8 -1*rhoa*t126)-6.491760000000001d-3*gammaaa*t11*(2.2222222 1715 9 22222222d-1*t53+1.111111111111111d-1*rhoa*rhob*(t234+t233 1716 : +t232+t180+t178+t176+t175+t170)+1.111111111111111d-1*rhob 1717 ; *t141+2.222222222222222d-1*rhoa*t124)-1.0962418720000001d 1718 < -3*gammabb*t11*t127*t29-1.0962418720000001d-3*gammaaa*t11 1719 = *t125*t29-1.0962418720000001d-3*gammaab*t11*t121*t29+t187 1720 > +t186+t185+t184+t183+t182+t167+t166+t165+t164-2.314450289 1721 ? 0117796d+0*t11*t139)+(t240+2.3266666666666663d-1*t138)*t5 1722 @ 8*t89+(t229-2*t120)*t5*t62+(t227+1.1633333333333332d-1*t1 1723 1 20)*t24*t60-2.3333333333333334d+0*t119*t24*t60+(t226+2.32 1724 2 66666666666663d-1*t119)*t58*t6-7.333333333333333d+0*t128* 1725 3 t49*t5+t24*(t254+t253+t252+t251+t250+t249+t248+t247+t246+ 1726 4 t219+t218+t217+t216+t215+t214+t213+t212+2.326666666666666 1727 5 3d-1*t128)*t44+t245+t244+t243+t242+t241+t230+t228+t225+t2 1728 6 24+t195+t194+t192+t190+t189+t163+t162+t160+t158+t157)*wgh 1729 7 t+Amat3(iq,D3_RA_RB_RB) 1730 Amat3(iq,D3_RB_RB_RB) = 1.0d+0*(1.711111111111111d+1*t5*t57* 1731 1 t94+t24*(-4.2655555555555547d-1*t57+t255+t239-5*t146)*t92 1732 2 +t5*t6*(-1.6874680727699207d-3*rhoa*t11*t68*t9+2.66477390 1733 3 09394716d-2*rhoa*t11*t23*t9-4.4999148607197886d-3*rhoa*t1 1734 4 1*t68*t8+7.106063735838591d-2*rhoa*t11*t23*t8-6.491760000 1735 5 000001d-3*gammaab*t11*(t169-2.333333333333333d+0*rhoa*t69 1736 6 )-4.6279677696266674d-5*gammabb*t11*t56*t68-4.62796776962 1737 7 66674d-5*gammaaa*t11*t54*t68-4.6279677696266674d-5*gammaa 1738 8 b*t11*t52*t68+7.308279146666667d-4*gammabb*t11*t23*t56+7. 1739 9 308279146666667d-4*gammaaa*t11*t23*t54+7.308279146666667d 1740 : -4*gammaab*t11*t23*t52-1.0962418720000001d-3*gammabb*t11* 1741 ; t144*t29-1.0962418720000001d-3*gammaaa*t11*t142*t29-1.096 1742 < 2418720000001d-3*gammaab*t11*t140*t29-3.9083350547112256d 1743 = -1*rhoa*t11*t139*t29-6.491760000000001d-3*gammabb*t11*(1. 1744 > 111111111111111d-1*rhoa*rhob*(t181+t179+t177+t174+t173+t1 1745 ? 72+t171+t170)+3.333333333333333d-1*rhoa*t143)-6.491760000 1746 @ 000001d-3*gammaaa*t11*(1.111111111111111d-1*rhoa*rhob*(t1 1747 1 80+t178+t176+t175+t170)+3.333333333333333d-1*rhoa*t141)-3 1748 2 .857417148352966d+0*rhoa*rhob**6.666666666666666d-1*t11)+ 1749 3 (t240+2.3266666666666663d-1*t146)*t58*t89-1.18032d+0*rhoa 1750 4 *t5*t62+2.288509333333333d-1*rhoa*t24*t60-1.5973795146666 1751 5 664d-2*rhoa*t58*t6-7.333333333333333d+0*t145*t49*t5+t24*( 1752 6 t254+t253+t252+t251+t250+t249+t248+t247+t246+t219+t218+t2 1753 7 17+t216+t215+t214+t213+t212+2.3266666666666663d-1*t145)*t 1754 8 44+t245+t244+t243+t242+t241+t225+t224+t195+t194+t192+t190 1755 9 +t189+t162+t160+t158+t157)*wght+Amat3(iq,D3_RB_RB_RB) 1756 Cmat3(iq,D3_RA_RA_GAA) = 1.0d+0*(t5*t6*(t261-6.4917600000000 1757 1 01d-3*t11*t84)+(4.760624000000001d-2*t11*t38+t264)*t49*t5 1758 2 +(t263-5.481209360000001d-4*t11*t38)*t44*t5+t24*(t262-1.5 1759 3 1041616d-3*t11*t38)*t44+t260+t259+t258+t257+t256)*wght+Cm 1760 4 at3(iq,D3_RA_RA_GAA) 1761 Cmat3(iq,D3_RA_RA_GAB) = 1.0d+0*(t5*t6*(t270-6.4917600000000 1762 1 01d-3*t11*t71)+(4.760624000000001d-2*t11*t32+t273)*t49*t5 1763 2 +(t272-5.481209360000001d-4*t11*t32)*t44*t5+t24*(t271-1.5 1764 3 1041616d-3*t11*t32)*t44+t269+t268+t267+t266+t265)*wght+Cm 1765 4 at3(iq,D3_RA_RA_GAB) 1766 Cmat3(iq,D3_RA_RA_GBB) = 1.0d+0*(t5*t6*(t279-6.4917600000000 1767 1 01d-3*t11*t77)+(4.760624000000001d-2*t11*t42+t282)*t49*t5 1768 2 +(t281-5.481209360000001d-4*t11*t42)*t44*t5+t24*(t280-1.5 1769 3 1041616d-3*t11*t42)*t44+t278+t277+t276+t275+t274)*wght+Cm 1770 4 at3(iq,D3_RA_RA_GBB) 1771 Cmat3(iq,D3_RA_RB_GAA) = 1.0d+0*((t261-6.491760000000001d-3* 1772 1 t11*t125)*t5*t6+t49*t5*(2.3803120000000005d-2*t11*t54+2.3 1773 2 803120000000005d-2*t11*t38+t264)+t24*t44*(-7.5520808d-4*t 1774 3 11*t54-7.5520808d-4*t11*t38+t262)+t283+t260+t259+t258+t25 1775 4 7+t256)*wght+Cmat3(iq,D3_RA_RB_GAA) 1776 Cmat3(iq,D3_RA_RB_GAB) = 1.0d+0*((t270-6.491760000000001d-3* 1777 1 t11*t121)*t5*t6+t49*t5*(2.3803120000000005d-2*t11*t52+2.3 1778 2 803120000000005d-2*t11*t32+t273)+t24*t44*(-7.5520808d-4*t 1779 3 11*t52-7.5520808d-4*t11*t32+t271)+t284+t269+t268+t267+t26 1780 4 6+t265)*wght+Cmat3(iq,D3_RA_RB_GAB) 1781 Cmat3(iq,D3_RA_RB_GBB) = 1.0d+0*((t279-6.491760000000001d-3* 1782 1 t11*t127)*t5*t6+t49*t5*(2.3803120000000005d-2*t11*t56+2.3 1783 2 803120000000005d-2*t11*t42+t282)+t24*t44*(-7.5520808d-4*t 1784 3 11*t56-7.5520808d-4*t11*t42+t280)+t285+t278+t277+t276+t27 1785 4 5+t274)*wght+Cmat3(iq,D3_RA_RB_GBB) 1786 Cmat3(iq,D3_RB_RB_GAA) = 1.0d+0*(t5*(-5.481209360000001d-4*t 1787 1 11*t29*t54-6.491760000000001d-3*t11*t142)*t6+t49*t5*(4.76 1788 2 0624000000001d-2*t11*t54+t264)+t24*t44*(t262-1.51041616d- 1789 3 3*t11*t54)+t283+t260+t259+t258+t257+t256)*wght+Cmat3(iq,D 1790 4 3_RB_RB_GAA) 1791 Cmat3(iq,D3_RB_RB_GAB) = 1.0d+0*(t5*(-5.481209360000001d-4*t 1792 1 11*t29*t52-6.491760000000001d-3*t11*t140)*t6+t49*t5*(4.76 1793 2 0624000000001d-2*t11*t52+t273)+t24*t44*(t271-1.51041616d- 1794 3 3*t11*t52)+t284+t269+t268+t267+t266+t265)*wght+Cmat3(iq,D 1795 4 3_RB_RB_GAB) 1796 Cmat3(iq,D3_RB_RB_GBB) = 1.0d+0*(t5*(-5.481209360000001d-4*t 1797 1 11*t29*t56-6.491760000000001d-3*t11*t144)*t6+t49*t5*(4.76 1798 2 0624000000001d-2*t11*t56+t282)+t24*t44*(t280-1.51041616d- 1799 3 3*t11*t56)+t285+t278+t277+t276+t275+t274)*wght+Cmat3(iq,D 1800 4 3_RB_RB_GBB) 1801 Cmat3(iq,D3_RA_GAA_GAA) = Cmat3(iq,D3_RA_GAA_GAA) 1802 Cmat3(iq,D3_RA_GAA_GAB) = Cmat3(iq,D3_RA_GAA_GAB) 1803 Cmat3(iq,D3_RA_GAA_GBB) = Cmat3(iq,D3_RA_GAA_GBB) 1804 Cmat3(iq,D3_RA_GAB_GAB) = Cmat3(iq,D3_RA_GAB_GAB) 1805 Cmat3(iq,D3_RA_GAB_GBB) = Cmat3(iq,D3_RA_GAB_GBB) 1806 Cmat3(iq,D3_RA_GBB_GBB) = Cmat3(iq,D3_RA_GBB_GBB) 1807 Cmat3(iq,D3_RB_GAA_GAA) = Cmat3(iq,D3_RB_GAA_GAA) 1808 Cmat3(iq,D3_RB_GAA_GAB) = Cmat3(iq,D3_RB_GAA_GAB) 1809 Cmat3(iq,D3_RB_GAA_GBB) = Cmat3(iq,D3_RB_GAA_GBB) 1810 Cmat3(iq,D3_RB_GAB_GAB) = Cmat3(iq,D3_RB_GAB_GAB) 1811 Cmat3(iq,D3_RB_GAB_GBB) = Cmat3(iq,D3_RB_GAB_GBB) 1812 Cmat3(iq,D3_RB_GBB_GBB) = Cmat3(iq,D3_RB_GBB_GBB) 1813 Cmat3(iq,D3_GAA_GAA_GAA) = Cmat3(iq,D3_GAA_GAA_GAA) 1814 Cmat3(iq,D3_GAA_GAA_GAB) = Cmat3(iq,D3_GAA_GAA_GAB) 1815 Cmat3(iq,D3_GAA_GAA_GBB) = Cmat3(iq,D3_GAA_GAA_GBB) 1816 Cmat3(iq,D3_GAA_GAB_GAB) = Cmat3(iq,D3_GAA_GAB_GAB) 1817 Cmat3(iq,D3_GAA_GAB_GBB) = Cmat3(iq,D3_GAA_GAB_GBB) 1818 Cmat3(iq,D3_GAA_GBB_GBB) = Cmat3(iq,D3_GAA_GBB_GBB) 1819 Cmat3(iq,D3_GAB_GAB_GAB) = Cmat3(iq,D3_GAB_GAB_GAB) 1820 Cmat3(iq,D3_GAB_GAB_GBB) = Cmat3(iq,D3_GAB_GAB_GBB) 1821 Cmat3(iq,D3_GAB_GBB_GBB) = Cmat3(iq,D3_GAB_GBB_GBB) 1822 Cmat3(iq,D3_GBB_GBB_GBB) = Cmat3(iq,D3_GBB_GBB_GBB) 1823 elseif (rhoa.gt.tol_rho.and.rhob.le.tol_rho) then 1824 fnc(iq) = fnc(iq) 1825 Amat(iq,D1_RA) = Amat(iq,D1_RA) 1826 Cmat(iq,D1_GAA) = Cmat(iq,D1_GAA) 1827 Amat2(iq,D2_RA_RA) = Amat2(iq,D2_RA_RA) 1828 Cmat2(iq,D2_RA_GAA) = Cmat2(iq,D2_RA_GAA) 1829 Cmat2(iq,D2_GAA_GAA) = Cmat2(iq,D2_GAA_GAA) 1830 Amat3(iq,D3_RA_RA_RA) = Amat3(iq,D3_RA_RA_RA) 1831 Cmat3(iq,D3_RA_RA_GAA) = Cmat3(iq,D3_RA_RA_GAA) 1832 Cmat3(iq,D3_RA_GAA_GAA) = Cmat3(iq,D3_RA_GAA_GAA) 1833 Cmat3(iq,D3_GAA_GAA_GAA) = Cmat3(iq,D3_GAA_GAA_GAA) 1834 elseif (rhoa.le.tol_rho.and.rhob.gt.tol_rho) then 1835 fnc(iq) = fnc(iq) 1836 Amat(iq,D1_RB) = Amat(iq,D1_RB) 1837 Cmat(iq,D1_GBB) = Cmat(iq,D1_GBB) 1838 Amat2(iq,D2_RB_RB) = Amat2(iq,D2_RB_RB) 1839 Cmat2(iq,D2_RB_GBB) = Cmat2(iq,D2_RB_GBB) 1840 Cmat2(iq,D2_GBB_GBB) = Cmat2(iq,D2_GBB_GBB) 1841 Amat3(iq,D3_RB_RB_RB) = Amat3(iq,D3_RB_RB_RB) 1842 Cmat3(iq,D3_RB_RB_GBB) = Cmat3(iq,D3_RB_RB_GBB) 1843 Cmat3(iq,D3_RB_GBB_GBB) = Cmat3(iq,D3_RB_GBB_GBB) 1844 Cmat3(iq,D3_GBB_GBB_GBB) = Cmat3(iq,D3_GBB_GBB_GBB) 1845 endif ! rhoa.gt.tol_rho.and.rhob.gt.tol_rho 1846 endif ! ipol.eq.1 1847 enddo ! iq 1848 end 1849C> @} 1850