1C> \ingroup nwxc 2C> @{ 3C> 4C> \file nwxcm_x_pw91.F 5C> The nwxcm_x_pw91 functional 6C> 7C> @} 8C> 9C> \ingroup nwxc_priv 10C> @{ 11C> 12C> \brief Evaluate the nwxcm_x_pw91 functional [1] 13C> 14C> \f{eqnarray*}{ 15C> {\it t_1} &=& 1.0\,\rho_\alpha^{{{1}\over{3}}}\\\\ 16C> {\it t_2} &=& {\it t_1}^{4.0}\\\\ 17C> {\it t_3} &=& 1.0\,\rho_\beta^{{{1}\over{3}}}\\\\ 18C> {\it t_4} &=& {\it t_3}^{4.0}\\\\ 19C> {\it t_5} &=& {\it param}\left(2\right)\\\\ 20C> {\it t_6} &=& {{1}\over{{\it t_2}}}\\\\ 21C> {\it t_7} &=& \sqrt{\sigma_{\alpha\alpha}}\\\\ 22C> {\it t_8} &=& {\it param}\left(1\right)\\\\ 23C> {\it t_9} &=& {{1}\over{{\it t_2}^{{\it t_8}}}}\\\\ 24C> {\it t_{10}} &=& {{{\it t_8}}\over{2}}\\\\ 25C> {\it t_{11}} &=& \sigma_{\alpha\alpha}^{{\it t_{10}}}\\\\ 26C> {\it t_{12}} &=& {{1}\over{{\it t_1}^{8.0}}}\\\\ 27C> {\it t_{13}} &=& {\it t_5}-0.001890381166699926\\\\ 28C> {\it t_{14}} &=& {{1}\over{{\it t_4}}}\\\\ 29C> {\it t_{15}} &=& \sqrt{\sigma_{\beta\beta}}\\\\ 30C> {\it t_{16}} &=& {{1}\over{{\it t_4}^{{\it t_8}}}}\\\\ 31C> {\it t_{17}} &=& \sigma_{\beta\beta}^{{\it t_{10}}}\\\\ 32C> {\it t_{18}} &=& {{1}\over{{\it t_3}^{8.0}}}\\\\ 33C> {\it t_{19}} &=& 1.0\,\rho_s^{{{1}\over{3}}}\\\\ 34C> {\it t_{20}} &=& {\it t_{19}}^{4.0}\\\\ 35C> {\it t_{21}} &=& {{1}\over{{\it t_{20}}}}\\\\ 36C> {\it t_{22}} &=& \sqrt{\sigma_{ss}}\\\\ 37C> {\it t_{23}} &=& {{1}\over{{\it t_{20}}^{{\it t_8}}}}\\\\ 38C> {\it t_{24}} &=& \sigma_{ss}^{{\it t_{10}}}\\\\ 39C> {\it t_{25}} &=& {{1}\over{{\it t_{19}}^{8.0}}}\\\\ 40C> f &=& {{1.0\,{\it t_4}\,\left({{{\it t_{13}}\,{\it t_{18}} 41C> \,\sigma_{\beta\beta}}\over{e^{1.6455\,{\it t_{18}}\, 42C> \sigma_{\beta\beta}}}}+1.0 \times 10^{-6}\,{\it t_{16}}\,{ 43C> \it t_{17}}-{\it t_5}\,{\it t_{18}}\, 44C> \sigma_{\beta\beta}\right)}\over{1.074661302677647 \times 10^{ 45C> -6}\,{\it t_{16}}\,{\it t_{17}}+6.0\,{\it t_5}\,{\it t_{14}} 46C> \,{\rm asinh}\; \left({\it t_{14}}\,{\it t_{15}}\right)\,{ 47C> \it t_{15}}+1.0}}+{{1.0\,{\it t_2}\,\left({{{\it t_{13}}\,{ 48C> \it t_{12}}\,\sigma_{\alpha\alpha}}\over{e^{1.6455\,{ 49C> \it t_{12}}\,\sigma_{\alpha\alpha}}}}+1.0 \times 10^{-6}\,{ 50C> \it t_9}\,{\it t_{11}}-{\it t_5}\,{\it t_{12}}\, 51C> \sigma_{\alpha\alpha}\right)}\over{1.074661302677647 \times 10^{ 52C> -6}\,{\it t_9}\,{\it t_{11}}+6.0\,{\it t_5}\,{\it t_6} 53C> \,{\rm asinh}\; \left({\it t_6}\,{\it t_7}\right)\,{\it t_7} 54C> +1.0}}-0.9305257363491\,{\it t_4}-0.9305257363491\,{\it t_2}\\\\ 55C> g &=& 0\\\\ 56C> G &=& {{1.0\,{\it t_{20}}\,\left({{{\it t_{13}}\,{\it t_{25}} 57C> \,\sigma_{ss}}\over{e^{1.6455\,{\it t_{25}}\,\sigma_{ss}}}} 58C> +1.0 \times 10^{-6}\,{\it t_{23}}\,{\it t_{24}}-{\it t_5}\,{ 59C> \it t_{25}}\,\sigma_{ss}\right)} 60C> \over{1.074661302677647 \times 10^{-6}\,{\it t_{23}}\,{ 61C> \it t_{24}}+6.0\,{\it t_5}\,{\it t_{21}}\,{\rm asinh}\; 62C> \left({\it t_{21}}\,{\it t_{22}}\right)\,{\it t_{22}}+1.0}} 63C> -0.9305257363491\,{\it t_{20}}\\\\ 64C> \f} 65C> 66C> Code generated with Maxima 5.34.0 [2] 67C> driven by autoxc [3]. 68C> 69C> ### References ### 70C> 71C> [1] JP Perdew, JA Chevary, SH Vosko, KA Jackson, MR Pederson 72C> , DJ Singh, C. Fiolhais, Phys.Rev. B 46, 6671 (1992) , DOI: 73C> <a href="https://doi.org/10.1021/jp050536c "> 74C> 10.1021/jp050536c </a> 75C> 76C> [2] Maxima, a computer algebra system, 77C> <a href="http://maxima.sourceforge.net/"> 78C> http://maxima.sourceforge.net/</a> 79C> 80C> [3] autoxc, revision 27097 2015-05-08 81C> 82 subroutine nwxcm_x_pw91(param,tol_rho,ipol,nq,wght, 83 +rho,rgamma,fnc,Amat,Cmat) 84c $Id: $ 85#ifdef NWXC_QUAD_PREC 86 implicit real(kind=selected_real_kind(30))(a-h,o-z),integer(i-n) 87 integer, parameter :: rk=selected_real_kind(30) 88#else 89 implicit real(kind=selected_real_kind(15))(a-h,o-z),integer(i-n) 90 integer, parameter :: rk=selected_real_kind(15) 91#endif 92#include "nwxc_param.fh" 93 double precision param(*) !< [Input] Parameters of functional 94 double precision tol_rho !< [Input] The lower limit on the density 95 integer ipol !< [Input] The number of spin channels 96 integer nq !< [Input] The number of points 97 double precision wght !< [Input] The weight of the functional 98 double precision rho(nq,*) !< [Input] The density 99 double precision rgamma(nq,*) !< [Input] The norm of the density 100 !< gradients 101 double precision fnc(nq) !< [Output] The value of the functional 102c 103c Sampling Matrices for the XC Kernel 104c 105 double precision Amat(nq,*) !< [Output] The derivative wrt rho 106 double precision Cmat(nq,*) !< [Output] The derivative wrt rgamma 107 integer iq 108 double precision tmp 109 double precision rhoa,rhob 110 double precision gammaaa,gammaab,gammabb 111 double precision taua,taub 112 double precision nwxcm_heaviside 113 external nwxcm_heaviside 114CDIR$ NOVECTOR 115 do iq = 1, nq 116 if (ipol.eq.1) then 117 rhoa = 0.5d0*rho(iq,R_T) 118 gammaaa = 0.25d0*rgamma(iq,G_TT) 119 if (rhoa.gt.tol_rho) then 120 t1 = 1.0d+0*rhoa**3.333333333333333d-1 121 t2 = t1**4.0d+0 122 t3 = param(1) 123 t4 = 5.0d-1*t3 124 t5 = gammaaa**t4 125 t6 = 1/t2**t3 126 t7 = param(2) 127 t8 = gammaaa**5.0d-1 128 t9 = 1/t2 129 t10 = asinh(t8*t9) 130 t11 = 6.0d+0*t10*t7*t8*t9+1.0746613026776465d-6*t5*t6+1.0d+0 131 t12 = 1/t11 132 t13 = 1/t1**8.0d+0 133 t14 = t7-1.890381166699926d-3 134 t15 = exp(-1.6455d+0*gammaaa*t13) 135 t16 = -gammaaa*t13*t7+1.0d-6*t5*t6+gammaaa*t13*t14*t15 136 t17 = t1**3.0d+0 137 t18 = 1/rhoa**6.666666666666666d-1 138 t19 = 1/t11**2 139 t20 = 1/rhoa 140 t21 = 1/(gammaaa*t13+1)**5.0d-1 141 t22 = 1/t1**9.0d+0 142 t23 = gammaaa**(t4-1) 143 fnc(iq) = (2.0d+0*t12*t16*t2-1.8610514726982d+0*t2)*wght+fnc 144 1 (iq) 145 Amat(iq,D1_RA) = (-1.0d+0*t16*t19*t2*(-8.0d+0*t10*t18*t7*t8/ 146 1 t1**5.0d+0-8.0d+0*gammaaa*t18*t21*t22*t7-1.43288173690352 147 2 88d-6*t20*t3*t5*t6)+1.0d+0*t12*t2*(2.6666666666666666d+0* 148 3 gammaaa*t18*t22*t7-1.3333333333333333d-6*t20*t3*t5*t6-2.6 149 4 666666666666666d+0*gammaaa*t14*t15*t18*t22+4.387999999999 150 5 9997d+0*gammaaa**2*t14*t15*t18/t1**1.7d+1)+1.333333333333 151 6 3333d+0*t12*t16*t17*t18-1.2407009817987999d+0*t17*t18)*wg 152 7 ht+Amat(iq,D1_RA) 153 Cmat(iq,D1_GAA) = (1.0d+0*t12*t2*(-t13*t7+5.0d-7*t23*t3*t6+t 154 1 13*t14*t15-1.6455d+0*gammaaa*t14*t15/t1**1.6d+1)-1.0d+0*t 155 2 16*t19*t2*(3.0d+0*t10*t7*t9/t8+3.0d+0*t13*t21*t7+5.373306 156 3 513388233d-7*t23*t3*t6))*wght+Cmat(iq,D1_GAA) 157 Cmat(iq,D1_GAB) = Cmat(iq,D1_GAB) 158 endif ! rhoa.gt.tol_rho 159 else ! ipol.eq.1 160 rhoa = rho(iq,R_A) 161 rhob = rho(iq,R_B) 162 gammaaa = rgamma(iq,G_AA) 163 gammaab = rgamma(iq,G_AB) 164 gammabb = rgamma(iq,G_BB) 165 if (rhoa.gt.tol_rho.and.rhob.gt.tol_rho) then 166 t1 = rhoa**3.333333333333333d-1 167 t2 = 1.0d+0*t1 168 t3 = t2**4.0d+0 169 t4 = param(1) 170 t5 = 5.0d-1*t4 171 t6 = gammaaa**t5 172 t7 = 1/t3**t4 173 t8 = param(2) 174 t9 = gammaaa**5.0d-1 175 t10 = 1/t3 176 t11 = asinh(t10*t9) 177 t12 = 6.0d+0*t10*t11*t8*t9+1.0746613026776465d-6*t6*t7+1.0d+ 178 1 0 179 t13 = 1/t12 180 t14 = 1/t2**8.0d+0 181 t15 = t8-1.890381166699926d-3 182 t16 = exp(-1.6455d+0*gammaaa*t14) 183 t17 = -gammaaa*t14*t8+1.0d-6*t6*t7+gammaaa*t14*t15*t16 184 t18 = rhob**3.333333333333333d-1 185 t19 = 1.0d+0*t18 186 t20 = t19**4.0d+0 187 t21 = gammabb**t5 188 t22 = 1/t20**t4 189 t23 = gammabb**5.0d-1 190 t24 = 1/t20 191 t25 = asinh(t23*t24) 192 t26 = 6.0d+0*t23*t24*t25*t8+1.0746613026776465d-6*t21*t22+1. 193 1 0d+0 194 t27 = 1/t26 195 t28 = 1/t19**8.0d+0 196 t29 = exp(-1.6455d+0*gammabb*t28) 197 t30 = -gammabb*t28*t8+gammabb*t15*t28*t29+1.0d-6*t21*t22 198 t31 = t2**3.0d+0 199 t32 = 1/rhoa**6.666666666666666d-1 200 t33 = 1/t12**2 201 t34 = 1/(gammaaa*t14+1)**5.0d-1 202 t35 = 1/t2**9.0d+0 203 t36 = 1.0d+0/t1 204 t37 = t19**3.0d+0 205 t38 = 1/rhob**6.666666666666666d-1 206 t39 = 1/t26**2 207 t40 = 1/(gammabb*t28+1)**5.0d-1 208 t41 = 1/t19**9.0d+0 209 t42 = 1.0d+0/t18 210 t43 = t5-1 211 t44 = gammaaa**t43 212 t45 = gammabb**t43 213 fnc(iq) = (1.0d+0*t20*t27*t30+1.0d+0*t13*t17*t3-9.3052573634 214 1 91d-1*t3-9.305257363491d-1*t20)*wght+fnc(iq) 215 Amat(iq,D1_RA) = (-1.0d+0*t17*t3*t33*(-8.0d+0*t11*t32*t8*t9/ 216 1 t2**5.0d+0-8.0d+0*gammaaa*t32*t34*t35*t8-1.43288173690352 217 2 88d-6*t32*t36*t4*t6*t7)+1.0d+0*t13*t3*(2.6666666666666666 218 3 d+0*gammaaa*t32*t35*t8-1.3333333333333333d-6*t32*t36*t4*t 219 4 6*t7-2.6666666666666666d+0*gammaaa*t15*t16*t32*t35+4.3879 220 5 999999999997d+0*gammaaa**2*t15*t16*t32/t2**1.7d+1)+1.3333 221 6 333333333333d+0*t13*t17*t31*t32-1.2407009817987999d+0*t31 222 7 *t32)*wght+Amat(iq,D1_RA) 223 Amat(iq,D1_RB) = (-1.0d+0*t20*t30*t39*(-8.0d+0*gammabb*t38*t 224 1 40*t41*t8-8.0d+0*t23*t25*t38*t8/t19**5.0d+0-1.43288173690 225 2 35288d-6*t21*t22*t38*t4*t42)+1.0d+0*t20*t27*(2.6666666666 226 3 666666d+0*gammabb*t38*t41*t8-1.3333333333333333d-6*t21*t2 227 4 2*t38*t4*t42-2.6666666666666666d+0*gammabb*t15*t29*t38*t4 228 5 1+4.3879999999999997d+0*gammabb**2*t15*t29*t38/t19**1.7d+ 229 6 1)+1.3333333333333333d+0*t27*t30*t37*t38-1.24070098179879 230 7 99d+0*t37*t38)*wght+Amat(iq,D1_RB) 231 Cmat(iq,D1_GAA) = (1.0d+0*t13*t3*(-t14*t8+5.0d-7*t4*t44*t7-1 232 1 .6455d+0*gammaaa*t15*t16/t2**1.6d+1+t14*t15*t16)-1.0d+0*t 233 2 17*t3*t33*(3.0d+0*t10*t11*t8/t9+3.0d+0*t14*t34*t8+5.37330 234 3 6513388233d-7*t4*t44*t7))*wght+Cmat(iq,D1_GAA) 235 Cmat(iq,D1_GAB) = Cmat(iq,D1_GAB) 236 Cmat(iq,D1_GBB) = (1.0d+0*t20*t27*(-t28*t8+5.0d-7*t22*t4*t45 237 1 +t15*t28*t29-1.6455d+0*gammabb*t15*t29/t19**1.6d+1)-1.0d+ 238 2 0*t20*t30*t39*(3.0d+0*t28*t40*t8+3.0d+0*t24*t25*t8/t23+5. 239 3 373306513388233d-7*t22*t4*t45))*wght+Cmat(iq,D1_GBB) 240 elseif (rhoa.gt.tol_rho.and.rhob.le.tol_rho) then 241 t1 = rhoa**3.333333333333333d-1 242 t2 = 1.0d+0*t1 243 t3 = t2**4.0d+0 244 t4 = param(1) 245 t5 = 5.0d-1*t4 246 t6 = gammaaa**t5 247 t7 = 1/t3**t4 248 t8 = param(2) 249 t9 = gammaaa**5.0d-1 250 t10 = 1/t3 251 t11 = asinh(t10*t9) 252 t12 = 6.0d+0*t10*t11*t8*t9+1.0746613026776465d-6*t6*t7+1.0d+ 253 1 0 254 t13 = 1/t12 255 t14 = 1/t2**8.0d+0 256 t15 = t8-1.890381166699926d-3 257 t16 = exp(-1.6455d+0*gammaaa*t14) 258 t17 = -gammaaa*t14*t8+1.0d-6*t6*t7+gammaaa*t14*t15*t16 259 t18 = t2**3.0d+0 260 t19 = 1/rhoa**6.666666666666666d-1 261 t20 = 1/t12**2 262 t21 = 1/(gammaaa*t14+1)**5.0d-1 263 t22 = 1/t2**9.0d+0 264 t23 = 1.0d+0/t1 265 t24 = gammaaa**(t5-1) 266 fnc(iq) = (1.0d+0*t13*t17*t3-9.305257363491d-1*t3)*wght+fnc( 267 1 iq) 268 Amat(iq,D1_RA) = (-1.0d+0*t17*t20*t3*(-8.0d+0*t11*t19*t8*t9/ 269 1 t2**5.0d+0-8.0d+0*gammaaa*t19*t21*t22*t8-1.43288173690352 270 2 88d-6*t19*t23*t4*t6*t7)+1.0d+0*t13*t3*(2.6666666666666666 271 3 d+0*gammaaa*t19*t22*t8-1.3333333333333333d-6*t19*t23*t4*t 272 4 6*t7-2.6666666666666666d+0*gammaaa*t15*t16*t19*t22+4.3879 273 5 999999999997d+0*gammaaa**2*t15*t16*t19/t2**1.7d+1)+1.3333 274 6 333333333333d+0*t13*t17*t18*t19-1.2407009817987999d+0*t18 275 7 *t19)*wght+Amat(iq,D1_RA) 276 Cmat(iq,D1_GAA) = (1.0d+0*t13*t3*(-t14*t8+5.0d-7*t24*t4*t7-1 277 1 .6455d+0*gammaaa*t15*t16/t2**1.6d+1+t14*t15*t16)-1.0d+0*t 278 2 17*t20*t3*(3.0d+0*t10*t11*t8/t9+3.0d+0*t14*t21*t8+5.37330 279 3 6513388233d-7*t24*t4*t7))*wght+Cmat(iq,D1_GAA) 280 elseif (rhoa.le.tol_rho.and.rhob.gt.tol_rho) then 281 t1 = rhob**3.333333333333333d-1 282 t2 = 1.0d+0*t1 283 t3 = t2**4.0d+0 284 t4 = param(1) 285 t5 = 5.0d-1*t4 286 t6 = gammabb**t5 287 t7 = 1/t3**t4 288 t8 = param(2) 289 t9 = gammabb**5.0d-1 290 t10 = 1/t3 291 t11 = asinh(t10*t9) 292 t12 = 6.0d+0*t10*t11*t8*t9+1.0746613026776465d-6*t6*t7+1.0d+ 293 1 0 294 t13 = 1/t12 295 t14 = 1/t2**8.0d+0 296 t15 = t8-1.890381166699926d-3 297 t16 = exp(-1.6455d+0*gammabb*t14) 298 t17 = -gammabb*t14*t8+1.0d-6*t6*t7+gammabb*t14*t15*t16 299 t18 = t2**3.0d+0 300 t19 = 1/rhob**6.666666666666666d-1 301 t20 = 1/t12**2 302 t21 = 1/(gammabb*t14+1)**5.0d-1 303 t22 = 1/t2**9.0d+0 304 t23 = 1.0d+0/t1 305 t24 = gammabb**(t5-1) 306 fnc(iq) = (1.0d+0*t13*t17*t3-9.305257363491d-1*t3)*wght+fnc( 307 1 iq) 308 Amat(iq,D1_RB) = (-1.0d+0*t17*t20*t3*(-8.0d+0*t11*t19*t8*t9/ 309 1 t2**5.0d+0-8.0d+0*gammabb*t19*t21*t22*t8-1.43288173690352 310 2 88d-6*t19*t23*t4*t6*t7)+1.0d+0*t13*t3*(2.6666666666666666 311 3 d+0*gammabb*t19*t22*t8-1.3333333333333333d-6*t19*t23*t4*t 312 4 6*t7-2.6666666666666666d+0*gammabb*t15*t16*t19*t22+4.3879 313 5 999999999997d+0*gammabb**2*t15*t16*t19/t2**1.7d+1)+1.3333 314 6 333333333333d+0*t13*t17*t18*t19-1.2407009817987999d+0*t18 315 7 *t19)*wght+Amat(iq,D1_RB) 316 Cmat(iq,D1_GBB) = (1.0d+0*t13*t3*(-t14*t8+5.0d-7*t24*t4*t7-1 317 1 .6455d+0*gammabb*t15*t16/t2**1.6d+1+t14*t15*t16)-1.0d+0*t 318 2 17*t20*t3*(3.0d+0*t10*t11*t8/t9+3.0d+0*t14*t21*t8+5.37330 319 3 6513388233d-7*t24*t4*t7))*wght+Cmat(iq,D1_GBB) 320 endif ! rhoa.gt.tol_rho.and.rhob.gt.tol_rho 321 endif ! ipol.eq.1 322 enddo ! iq 323 end 324C> 325C> \brief Evaluate the nwxcm_x_pw91 functional [1] 326C> 327C> \f{eqnarray*}{ 328C> {\it t_1} &=& 1.0\,\rho_\alpha^{{{1}\over{3}}}\\\\ 329C> {\it t_2} &=& {\it t_1}^{4.0}\\\\ 330C> {\it t_3} &=& 1.0\,\rho_\beta^{{{1}\over{3}}}\\\\ 331C> {\it t_4} &=& {\it t_3}^{4.0}\\\\ 332C> {\it t_5} &=& {\it param}\left(2\right)\\\\ 333C> {\it t_6} &=& {{1}\over{{\it t_2}}}\\\\ 334C> {\it t_7} &=& \sqrt{\sigma_{\alpha\alpha}}\\\\ 335C> {\it t_8} &=& {\it param}\left(1\right)\\\\ 336C> {\it t_9} &=& {{1}\over{{\it t_2}^{{\it t_8}}}}\\\\ 337C> {\it t_{10}} &=& {{{\it t_8}}\over{2}}\\\\ 338C> {\it t_{11}} &=& \sigma_{\alpha\alpha}^{{\it t_{10}}}\\\\ 339C> {\it t_{12}} &=& {{1}\over{{\it t_1}^{8.0}}}\\\\ 340C> {\it t_{13}} &=& {\it t_5}-0.001890381166699926\\\\ 341C> {\it t_{14}} &=& {{1}\over{{\it t_4}}}\\\\ 342C> {\it t_{15}} &=& \sqrt{\sigma_{\beta\beta}}\\\\ 343C> {\it t_{16}} &=& {{1}\over{{\it t_4}^{{\it t_8}}}}\\\\ 344C> {\it t_{17}} &=& \sigma_{\beta\beta}^{{\it t_{10}}}\\\\ 345C> {\it t_{18}} &=& {{1}\over{{\it t_3}^{8.0}}}\\\\ 346C> {\it t_{19}} &=& 1.0\,\rho_s^{{{1}\over{3}}}\\\\ 347C> {\it t_{20}} &=& {\it t_{19}}^{4.0}\\\\ 348C> {\it t_{21}} &=& {{1}\over{{\it t_{20}}}}\\\\ 349C> {\it t_{22}} &=& \sqrt{\sigma_{ss}}\\\\ 350C> {\it t_{23}} &=& {{1}\over{{\it t_{20}}^{{\it t_8}}}}\\\\ 351C> {\it t_{24}} &=& \sigma_{ss}^{{\it t_{10}}}\\\\ 352C> {\it t_{25}} &=& {{1}\over{{\it t_{19}}^{8.0}}}\\\\ 353C> f &=& {{1.0\,{\it t_4}\,\left({{{\it t_{13}}\,{\it t_{18}} 354C> \,\sigma_{\beta\beta}}\over{e^{1.6455\,{\it t_{18}}\, 355C> \sigma_{\beta\beta}}}}+1.0 \times 10^{-6}\,{\it t_{16}}\,{ 356C> \it t_{17}}-{\it t_5}\,{\it t_{18}}\, 357C> \sigma_{\beta\beta}\right)}\over{1.074661302677647 \times 10^{ 358C> -6}\,{\it t_{16}}\,{\it t_{17}}+6.0\,{\it t_5}\,{\it t_{14}} 359C> \,{\rm asinh}\; \left({\it t_{14}}\,{\it t_{15}}\right)\,{ 360C> \it t_{15}}+1.0}}+{{1.0\,{\it t_2}\,\left({{{\it t_{13}}\,{ 361C> \it t_{12}}\,\sigma_{\alpha\alpha}}\over{e^{1.6455\,{ 362C> \it t_{12}}\,\sigma_{\alpha\alpha}}}}+1.0 \times 10^{-6}\,{ 363C> \it t_9}\,{\it t_{11}}-{\it t_5}\,{\it t_{12}}\, 364C> \sigma_{\alpha\alpha}\right)}\over{1.074661302677647 \times 10^{ 365C> -6}\,{\it t_9}\,{\it t_{11}}+6.0\,{\it t_5}\,{\it t_6} 366C> \,{\rm asinh}\; \left({\it t_6}\,{\it t_7}\right)\,{\it t_7} 367C> +1.0}}-0.9305257363491\,{\it t_4}-0.9305257363491\,{\it t_2}\\\\ 368C> g &=& 0\\\\ 369C> G &=& {{1.0\,{\it t_{20}}\,\left({{{\it t_{13}}\,{\it t_{25}} 370C> \,\sigma_{ss}}\over{e^{1.6455\,{\it t_{25}}\,\sigma_{ss}}}} 371C> +1.0 \times 10^{-6}\,{\it t_{23}}\,{\it t_{24}}-{\it t_5}\,{ 372C> \it t_{25}}\,\sigma_{ss}\right)} 373C> \over{1.074661302677647 \times 10^{-6}\,{\it t_{23}}\,{ 374C> \it t_{24}}+6.0\,{\it t_5}\,{\it t_{21}}\,{\rm asinh}\; 375C> \left({\it t_{21}}\,{\it t_{22}}\right)\,{\it t_{22}}+1.0}} 376C> -0.9305257363491\,{\it t_{20}}\\\\ 377C> \f} 378C> 379C> Code generated with Maxima 5.34.0 [2] 380C> driven by autoxc [3]. 381C> 382C> ### References ### 383C> 384C> [1] JP Perdew, JA Chevary, SH Vosko, KA Jackson, MR Pederson 385C> , DJ Singh, C. Fiolhais, Phys.Rev. B 46, 6671 (1992) , DOI: 386C> <a href="https://doi.org/10.1021/jp050536c "> 387C> 10.1021/jp050536c </a> 388C> 389C> [2] Maxima, a computer algebra system, 390C> <a href="http://maxima.sourceforge.net/"> 391C> http://maxima.sourceforge.net/</a> 392C> 393C> [3] autoxc, revision 27097 2015-05-08 394C> 395 subroutine nwxcm_x_pw91_d2(param,tol_rho,ipol,nq,wght, 396 +rho,rgamma,fnc,Amat,Amat2,Cmat,Cmat2) 397c $Id: $ 398#ifdef NWXC_QUAD_PREC 399 implicit real(kind=selected_real_kind(30))(a-h,o-z),integer(i-n) 400 integer, parameter :: rk=selected_real_kind(30) 401#else 402 implicit real(kind=selected_real_kind(15))(a-h,o-z),integer(i-n) 403 integer, parameter :: rk=selected_real_kind(15) 404#endif 405#include "nwxc_param.fh" 406 double precision param(*) !< [Input] Parameters of functional 407 double precision tol_rho !< [Input] The lower limit on the density 408 integer ipol !< [Input] The number of spin channels 409 integer nq !< [Input] The number of points 410 double precision wght !< [Input] The weight of the functional 411 double precision rho(nq,*) !< [Input] The density 412 double precision rgamma(nq,*) !< [Input] The norm of the density 413 !< gradients 414 double precision fnc(nq) !< [Output] The value of the functional 415c 416c Sampling Matrices for the XC Kernel 417c 418 double precision Amat(nq,*) !< [Output] The derivative wrt rho 419 double precision Cmat(nq,*) !< [Output] The derivative wrt rgamma 420c 421c Sampling Matrices for the XC Kernel 422c 423 double precision Amat2(nq,*) !< [Output] The 2nd derivative wrt rho 424 double precision Cmat2(nq,*) !< [Output] The 2nd derivative wrt rgamma 425 !< and possibly rho 426 integer iq 427 double precision tmp 428 double precision rhoa,rhob 429 double precision gammaaa,gammaab,gammabb 430 double precision taua,taub 431 double precision nwxcm_heaviside 432 external nwxcm_heaviside 433CDIR$ NOVECTOR 434 do iq = 1, nq 435 if (ipol.eq.1) then 436 rhoa = 0.5d0*rho(iq,R_T) 437 gammaaa = 0.25d0*rgamma(iq,G_TT) 438 if (rhoa.gt.tol_rho) then 439 t1 = 1.0d+0*rhoa**3.333333333333333d-1 440 t2 = t1**4.0d+0 441 t3 = param(1) 442 t4 = 5.0d-1*t3 443 t5 = gammaaa**t4 444 t6 = 1/t2**t3 445 t7 = param(2) 446 t8 = gammaaa**5.0d-1 447 t9 = 1/t2 448 t10 = asinh(t8*t9) 449 t11 = 6.0d+0*t10*t7*t8*t9+1.0746613026776465d-6*t5*t6+1.0d+0 450 t12 = 1/t11 451 t13 = 1/t1**8.0d+0 452 t14 = t7-1.890381166699926d-3 453 t15 = exp(-1.6455d+0*gammaaa*t13) 454 t16 = -gammaaa*t13*t7+1.0d-6*t5*t6+gammaaa*t13*t14*t15 455 t17 = t1**3.0d+0 456 t18 = 1/rhoa**6.666666666666666d-1 457 t19 = 1/t11**2 458 t20 = 1/rhoa 459 t21 = (gammaaa*t13+1)**5.0d-1 460 t22 = 1/t21 461 t23 = 1/t1**9.0d+0 462 t24 = 1/t1**5.0d+0 463 t25 = -8.0d+0*t10*t18*t24*t7*t8-8.0d+0*gammaaa*t18*t22*t23*t 464 1 7-1.4328817369035288d-6*t20*t3*t5*t6 465 t26 = gammaaa**2 466 t27 = 1/t1**1.7d+1 467 t28 = 2.6666666666666666d+0*gammaaa*t18*t23*t7-1.33333333333 468 1 33333d-6*t20*t3*t5*t6+4.3879999999999997d+0*t14*t15*t18*t 469 2 26*t27-2.6666666666666666d+0*gammaaa*t14*t15*t18*t23 470 t29 = t4-1 471 t30 = gammaaa**t29 472 t31 = 1/t1**1.6d+1 473 t32 = -t13*t7+5.0d-7*t3*t30*t6-1.6455d+0*gammaaa*t14*t15*t31 474 1 +t13*t14*t15 475 t33 = 1/t8 476 t34 = 3.0d+0*t10*t33*t7*t9+3.0d+0*t13*t22*t7+5.3733065133882 477 1 33d-7*t3*t30*t6 478 t35 = 1/rhoa**1.6666666666666669d+0 479 t36 = t1**2.0d+0 480 t37 = 1/rhoa**1.3333333333333333d+0 481 t38 = 1/t11**3 482 t39 = 1/rhoa**2 483 t40 = t3**2 484 t41 = 1/t1**1.0d+1 485 t42 = 1/t1**1.8d+1 486 t43 = 1/t21**3 487 t44 = gammaaa**(t4-2) 488 fnc(iq) = (2.0d+0*t12*t16*t2-1.8610514726982d+0*t2)*wght+fnc 489 1 (iq) 490 Amat(iq,D1_RA) = (1.0d+0*t12*t2*t28-1.0d+0*t16*t19*t2*t25+1. 491 1 3333333333333333d+0*t12*t16*t17*t18-1.2407009817987999d+0 492 2 *t17*t18)*wght+Amat(iq,D1_RA) 493 Cmat(iq,D1_GAA) = (1.0d+0*t12*t2*t32-1.0d+0*t16*t19*t2*t34)* 494 1 wght+Cmat(iq,D1_GAA) 495 Cmat(iq,D1_GAB) = Cmat(iq,D1_GAB) 496 Amat2(iq,D2_RA_RA) = (t19*(-1.0d+0*t16*t2*(1.333333333333333 497 1 3d+1*t10*t37*t7*t8/t1**6.0d+0+5.333333333333333d+0*t10*t2 498 2 4*t35*t7*t8-1.0666666666666666d+1*t26*t37*t42*t43*t7+3.46 499 3 6666666666666d+1*gammaaa*t22*t37*t41*t7+5.333333333333333 500 4 d+0*gammaaa*t22*t23*t35*t7+1.9105089825380384d-6*t39*t40* 501 5 t5*t6+1.4328817369035288d-6*t3*t39*t5*t6)-2.0d+0*t2*t25*t 502 6 28)+1.0d+0*t12*t2*(-8.0d+0*gammaaa*t37*t41*t7-1.777777777 503 7 7777776d+0*gammaaa*t23*t35*t7+1.7777777777777776d-6*t39*t 504 8 40*t5*t6+1.3333333333333333d-6*t3*t39*t5*t6-3.65666666666 505 9 66664d+1*t14*t15*t26*t37*t42+8.0d+0*gammaaa*t14*t15*t37*t 506 : 41+1.9254543999999998d+1*gammaaa**3*t14*t15*t37/t1**2.6d+ 507 ; 1-2.925333333333333d+0*t14*t15*t26*t27*t35+1.777777777777 508 < 7776d+0*gammaaa*t14*t15*t23*t35)+2.0d+0*t16*t2*t25**2*t38 509 = +1.3333333333333333d+0*t12*t16*t36*t37-1.2407009817987999 510 > d+0*t36*t37-8.888888888888888d-1*t12*t16*t17*t35+8.271339 511 ? 878658666d-1*t17*t35+2.6666666666666666d+0*t12*t17*t18*t2 512 @ 8-2.6666666666666666d+0*t16*t17*t18*t19*t25)*wght+Amat2(i 513 1 q,D2_RA_RA) 514 Amat2(iq,D2_RA_RB) = Amat2(iq,D2_RA_RB) 515 Cmat2(iq,D2_RA_GAA) = (t19*(-1.0d+0*t16*t2*(4.0d+0*gammaaa*t 516 1 18*t27*t43*t7-4.0d+0*t10*t18*t24*t33*t7-1.2d+1*t18*t22*t2 517 2 3*t7-7.164408684517644d-7*t20*t30*t40*t6)-1.0d+0*t2*t28*t 518 3 34-1.0d+0*t2*t25*t32)+1.0d+0*t12*t2*(2.6666666666666666d+ 519 4 0*t18*t23*t7-6.666666666666666d-7*t20*t30*t40*t6+1.3164d+ 520 5 1*gammaaa*t14*t15*t18*t27-7.220454d+0*t14*t15*t18*t26/t1* 521 6 *2.5d+1-2.6666666666666666d+0*t14*t15*t18*t23)+2.0d+0*t16 522 7 *t2*t25*t34*t38-1.3333333333333333d+0*t16*t17*t18*t19*t34 523 8 +1.3333333333333333d+0*t12*t17*t18*t32)*wght+Cmat2(iq,D2_ 524 9 RA_GAA) 525 Cmat2(iq,D2_RA_GAB) = Cmat2(iq,D2_RA_GAB) 526 Cmat2(iq,D2_RA_GBB) = Cmat2(iq,D2_RA_GBB) 527 Cmat2(iq,D2_GAA_GAA) = (t19*(-1.0d+0*t16*t2*(-1.5d+0*t10*t7* 528 1 t9/t8**3-1.5d+0*t31*t43*t7+1.5d+0*t13*t22*t7/gammaaa+5.37 529 2 3306513388233d-7*t29*t3*t44*t6)-2.0d+0*t2*t32*t34)+1.0d+0 530 3 *t12*t2*(5.0d-7*t29*t3*t44*t6-3.291d+0*t14*t15*t31+2.7076 531 4 7025d+0*gammaaa*t14*t15/t1**2.4d+1)+2.0d+0*t16*t2*t34**2* 532 5 t38)*wght+Cmat2(iq,D2_GAA_GAA) 533 Cmat2(iq,D2_GAA_GAB) = Cmat2(iq,D2_GAA_GAB) 534 Cmat2(iq,D2_GAA_GBB) = Cmat2(iq,D2_GAA_GBB) 535 Cmat2(iq,D2_GAB_GAB) = Cmat2(iq,D2_GAB_GAB) 536 endif ! rhoa.gt.tol_rho 537 else ! ipol.eq.1 538 rhoa = rho(iq,R_A) 539 rhob = rho(iq,R_B) 540 gammaaa = rgamma(iq,G_AA) 541 gammaab = rgamma(iq,G_AB) 542 gammabb = rgamma(iq,G_BB) 543 if (rhoa.gt.tol_rho.and.rhob.gt.tol_rho) then 544 t1 = rhoa**3.333333333333333d-1 545 t2 = 1.0d+0*t1 546 t3 = t2**4.0d+0 547 t4 = param(1) 548 t5 = 5.0d-1*t4 549 t6 = gammaaa**t5 550 t7 = 1/t3**t4 551 t8 = param(2) 552 t9 = gammaaa**5.0d-1 553 t10 = 1/t3 554 t11 = asinh(t10*t9) 555 t12 = 6.0d+0*t10*t11*t8*t9+1.0746613026776465d-6*t6*t7+1.0d+ 556 1 0 557 t13 = 1/t12 558 t14 = 1/t2**8.0d+0 559 t15 = t8-1.890381166699926d-3 560 t16 = exp(-1.6455d+0*gammaaa*t14) 561 t17 = -gammaaa*t14*t8+1.0d-6*t6*t7+gammaaa*t14*t15*t16 562 t18 = rhob**3.333333333333333d-1 563 t19 = 1.0d+0*t18 564 t20 = t19**4.0d+0 565 t21 = gammabb**t5 566 t22 = 1/t20**t4 567 t23 = gammabb**5.0d-1 568 t24 = 1/t20 569 t25 = asinh(t23*t24) 570 t26 = 6.0d+0*t23*t24*t25*t8+1.0746613026776465d-6*t21*t22+1. 571 1 0d+0 572 t27 = 1/t26 573 t28 = 1/t19**8.0d+0 574 t29 = exp(-1.6455d+0*gammabb*t28) 575 t30 = -gammabb*t28*t8+gammabb*t15*t28*t29+1.0d-6*t21*t22 576 t31 = t2**3.0d+0 577 t32 = 1/rhoa**6.666666666666666d-1 578 t33 = 1/t12**2 579 t34 = (gammaaa*t14+1)**5.0d-1 580 t35 = 1/t34 581 t36 = 1/t2**9.0d+0 582 t37 = -8.0d+0*gammaaa*t32*t35*t36*t8 583 t38 = 1/t2**5.0d+0 584 t39 = -8.0d+0*t11*t32*t38*t8*t9 585 t40 = 1.0d+0/t1 586 t41 = -1.4328817369035288d-6*t32*t4*t40*t6*t7+t39+t37 587 t42 = 2.6666666666666666d+0*gammaaa*t32*t36*t8 588 t43 = gammaaa**2 589 t44 = 1/t2**1.7d+1 590 t45 = 4.3879999999999997d+0*t15*t16*t32*t43*t44 591 t46 = -2.6666666666666666d+0*gammaaa*t15*t16*t32*t36 592 t47 = -1.3333333333333333d-6*t32*t4*t40*t6*t7+t46+t45+t42 593 t48 = t19**3.0d+0 594 t49 = 1/rhob**6.666666666666666d-1 595 t50 = 1/t26**2 596 t51 = (gammabb*t28+1)**5.0d-1 597 t52 = 1/t51 598 t53 = 1/t19**9.0d+0 599 t54 = -8.0d+0*gammabb*t49*t52*t53*t8 600 t55 = 1/t19**5.0d+0 601 t56 = -8.0d+0*t23*t25*t49*t55*t8 602 t57 = 1.0d+0/t18 603 t58 = -1.4328817369035288d-6*t21*t22*t4*t49*t57+t56+t54 604 t59 = 2.6666666666666666d+0*gammabb*t49*t53*t8 605 t60 = gammabb**2 606 t61 = 1/t19**1.7d+1 607 t62 = 4.3879999999999997d+0*t15*t29*t49*t60*t61 608 t63 = -2.6666666666666666d+0*gammabb*t15*t29*t49*t53 609 t64 = t63+t62+t59-1.3333333333333333d-6*t21*t22*t4*t49*t57 610 t65 = t5-1 611 t66 = gammaaa**t65 612 t67 = 1/t2**1.6d+1 613 t68 = -t14*t8+5.0d-7*t4*t66*t7-1.6455d+0*gammaaa*t15*t16*t67 614 1 +t14*t15*t16 615 t69 = 1/t9 616 t70 = 3.0d+0*t10*t11*t69*t8+3.0d+0*t14*t35*t8+5.373306513388 617 1 233d-7*t4*t66*t7 618 t71 = gammabb**t65 619 t72 = 1/t19**1.6d+1 620 t73 = -t28*t8-1.6455d+0*gammabb*t15*t29*t72+5.0d-7*t22*t4*t7 621 1 1+t15*t28*t29 622 t74 = 1/t23 623 t75 = 3.0d+0*t24*t25*t74*t8+3.0d+0*t28*t52*t8+5.373306513388 624 1 233d-7*t22*t4*t71 625 t76 = 1/rhoa**1.6666666666666669d+0 626 t77 = t2**2.0d+0 627 t78 = 1/rhoa**1.3333333333333333d+0 628 t79 = 1/t12**3 629 t80 = 1/rhoa 630 t81 = -1.4328817369035288d-6*t4*t6*t7*t80+t39+t37 631 t82 = 1/rhoa**2 632 t83 = t4**2 633 t84 = 1/t2**1.0d+1 634 t85 = 1/t2**1.8d+1 635 t86 = -1.3333333333333333d-6*t4*t6*t7*t80+t46+t45+t42 636 t87 = 1/t34**3 637 t88 = 1/rhob**1.6666666666666669d+0 638 t89 = t19**2.0d+0 639 t90 = 1/rhob**1.3333333333333333d+0 640 t91 = 1/t26**3 641 t92 = 1/rhob 642 t93 = -1.4328817369035288d-6*t21*t22*t4*t92+t56+t54 643 t94 = 1/rhob**2 644 t95 = 1/t19**1.0d+1 645 t96 = 1/t19**1.8d+1 646 t97 = -1.3333333333333333d-6*t21*t22*t4*t92+t63+t62+t59 647 t98 = 1/t51**3 648 t99 = t5-2 649 t100 = gammaaa**t99 650 t101 = gammabb**t99 651 fnc(iq) = (1.0d+0*t20*t27*t30+1.0d+0*t13*t17*t3-9.3052573634 652 1 91d-1*t3-9.305257363491d-1*t20)*wght+fnc(iq) 653 Amat(iq,D1_RA) = (1.0d+0*t13*t3*t47-1.0d+0*t17*t3*t33*t41+1. 654 1 3333333333333333d+0*t13*t17*t31*t32-1.2407009817987999d+0 655 2 *t31*t32)*wght+Amat(iq,D1_RA) 656 Amat(iq,D1_RB) = (1.0d+0*t20*t27*t64-1.0d+0*t20*t30*t50*t58+ 657 1 1.3333333333333333d+0*t27*t30*t48*t49-1.2407009817987999d 658 2 +0*t48*t49)*wght+Amat(iq,D1_RB) 659 Cmat(iq,D1_GAA) = (1.0d+0*t13*t3*t68-1.0d+0*t17*t3*t33*t70)* 660 1 wght+Cmat(iq,D1_GAA) 661 Cmat(iq,D1_GAB) = Cmat(iq,D1_GAB) 662 Cmat(iq,D1_GBB) = (1.0d+0*t20*t27*t73-1.0d+0*t20*t30*t50*t75 663 1 )*wght+Cmat(iq,D1_GBB) 664 Amat2(iq,D2_RA_RA) = (t33*(-1.0d+0*t17*t3*(1.333333333333333 665 1 3d+1*t11*t78*t8*t9/t2**6.0d+0+5.333333333333333d+0*t11*t3 666 2 8*t76*t8*t9-1.0666666666666666d+1*t43*t78*t8*t85*t87+3.46 667 3 6666666666666d+1*gammaaa*t35*t78*t8*t84+1.910508982538038 668 4 4d-6*t40*t6*t7*t76*t83+1.4328817369035288d-6*t4*t6*t7*t82 669 5 +5.333333333333333d+0*gammaaa*t35*t36*t76*t8)-1.0d+0*t3*t 670 6 41*t86-1.0d+0*t3*t47*t81)+t13*t32*(1.3333333333333333d+0* 671 7 t31*t86+1.3333333333333333d+0*t31*t47)+1.0d+0*t13*t3*(-3. 672 8 6566666666666664d+1*t15*t16*t43*t78*t85-8.0d+0*gammaaa*t7 673 9 8*t8*t84+8.0d+0*gammaaa*t15*t16*t78*t84+1.777777777777777 674 : 6d-6*t40*t6*t7*t76*t83+1.3333333333333333d-6*t4*t6*t7*t82 675 ; -1.7777777777777776d+0*gammaaa*t36*t76*t8+1.9254543999999 676 < 998d+1*gammaaa**3*t15*t16*t78/t2**2.6d+1-2.92533333333333 677 = 3d+0*t15*t16*t43*t44*t76+1.7777777777777776d+0*gammaaa*t1 678 > 5*t16*t36*t76)+t32*t33*(-1.3333333333333333d+0*t17*t31*t8 679 ? 1-1.3333333333333333d+0*t17*t31*t41)+2.0d+0*t17*t3*t41*t7 680 @ 9*t81+1.3333333333333333d+0*t13*t17*t77*t78-1.24070098179 681 1 87999d+0*t77*t78-8.888888888888888d-1*t13*t17*t31*t76+8.2 682 2 71339878658666d-1*t31*t76)*wght+Amat2(iq,D2_RA_RA) 683 Amat2(iq,D2_RA_RB) = Amat2(iq,D2_RA_RB) 684 Amat2(iq,D2_RB_RB) = (t50*(-1.0d+0*t20*t30*(-1.0666666666666 685 1 666d+1*t60*t8*t90*t96*t98+3.466666666666666d+1*gammabb*t5 686 2 2*t8*t90*t95+1.4328817369035288d-6*t21*t22*t4*t94+1.33333 687 3 33333333333d+1*t23*t25*t8*t90/t19**6.0d+0+1.9105089825380 688 4 384d-6*t21*t22*t57*t83*t88+5.333333333333333d+0*t23*t25*t 689 5 55*t8*t88+5.333333333333333d+0*gammabb*t52*t53*t8*t88)-1. 690 6 0d+0*t20*t58*t97-1.0d+0*t20*t64*t93)+t27*t49*(1.333333333 691 7 3333333d+0*t48*t97+1.3333333333333333d+0*t48*t64)+1.0d+0* 692 8 t20*t27*(-3.6566666666666664d+1*t15*t29*t60*t90*t96-8.0d+ 693 9 0*gammabb*t8*t90*t95+8.0d+0*gammabb*t15*t29*t90*t95+1.333 694 : 3333333333333d-6*t21*t22*t4*t94+1.9254543999999998d+1*gam 695 ; mabb**3*t15*t29*t90/t19**2.6d+1+1.7777777777777776d-6*t21 696 < *t22*t57*t83*t88-1.7777777777777776d+0*gammabb*t53*t8*t88 697 = -2.925333333333333d+0*t15*t29*t60*t61*t88+1.7777777777777 698 > 776d+0*gammabb*t15*t29*t53*t88)+t49*t50*(-1.3333333333333 699 ? 333d+0*t30*t48*t93-1.3333333333333333d+0*t30*t48*t58)+2.0 700 @ d+0*t20*t30*t58*t91*t93+1.3333333333333333d+0*t27*t30*t89 701 1 *t90-1.2407009817987999d+0*t89*t90-8.888888888888888d-1*t 702 2 27*t30*t48*t88+8.271339878658666d-1*t48*t88)*wght+Amat2(i 703 3 q,D2_RB_RB) 704 Cmat2(iq,D2_RA_GAA) = (t33*(-1.0d+0*t17*t3*(4.0d+0*gammaaa*t 705 1 32*t44*t8*t87-7.164408684517644d-7*t66*t7*t80*t83-4.0d+0* 706 2 t11*t32*t38*t69*t8-1.2d+1*t32*t35*t36*t8)-1.0d+0*t3*t70*t 707 3 86-1.0d+0*t3*t68*t81)+1.0d+0*t13*t3*(-6.666666666666666d- 708 4 7*t66*t7*t80*t83+2.6666666666666666d+0*t32*t36*t8+1.3164d 709 5 +1*gammaaa*t15*t16*t32*t44-7.220454d+0*t15*t16*t32*t43/t2 710 6 **2.5d+1-2.6666666666666666d+0*t15*t16*t32*t36)+2.0d+0*t1 711 7 7*t3*t70*t79*t81-1.3333333333333333d+0*t17*t31*t32*t33*t7 712 8 0+1.3333333333333333d+0*t13*t31*t32*t68)*wght+Cmat2(iq,D2 713 9 _RA_GAA) 714 Cmat2(iq,D2_RA_GAB) = Cmat2(iq,D2_RA_GAB) 715 Cmat2(iq,D2_RA_GBB) = Cmat2(iq,D2_RA_GBB) 716 Cmat2(iq,D2_RB_GAA) = Cmat2(iq,D2_RB_GAA) 717 Cmat2(iq,D2_RB_GAB) = Cmat2(iq,D2_RB_GAB) 718 Cmat2(iq,D2_RB_GBB) = (t50*(-1.0d+0*t20*t30*(4.0d+0*gammabb* 719 1 t49*t61*t8*t98-7.164408684517644d-7*t22*t71*t83*t92-4.0d+ 720 2 0*t25*t49*t55*t74*t8-1.2d+1*t49*t52*t53*t8)-1.0d+0*t20*t7 721 3 5*t97-1.0d+0*t20*t73*t93)+2.0d+0*t20*t30*t75*t91*t93+1.0d 722 4 +0*t20*t27*(-6.666666666666666d-7*t22*t71*t83*t92+2.66666 723 5 66666666666d+0*t49*t53*t8+1.3164d+1*gammabb*t15*t29*t49*t 724 6 61-7.220454d+0*t15*t29*t49*t60/t19**2.5d+1-2.666666666666 725 7 6666d+0*t15*t29*t49*t53)-1.3333333333333333d+0*t30*t48*t4 726 8 9*t50*t75+1.3333333333333333d+0*t27*t48*t49*t73)*wght+Cma 727 9 t2(iq,D2_RB_GBB) 728 Cmat2(iq,D2_GAA_GAA) = (t33*(-1.0d+0*t17*t3*(-1.5d+0*t10*t11 729 1 *t8/t9**3-1.5d+0*t67*t8*t87+1.5d+0*t14*t35*t8/gammaaa+5.3 730 2 73306513388233d-7*t100*t4*t65*t7)-2.0d+0*t3*t68*t70)+2.0d 731 3 +0*t17*t3*t70**2*t79+1.0d+0*t13*t3*(5.0d-7*t100*t4*t65*t7 732 4 -3.291d+0*t15*t16*t67+2.70767025d+0*gammaaa*t15*t16/t2**2 733 5 .4d+1))*wght+Cmat2(iq,D2_GAA_GAA) 734 Cmat2(iq,D2_GAA_GAB) = Cmat2(iq,D2_GAA_GAB) 735 Cmat2(iq,D2_GAA_GBB) = Cmat2(iq,D2_GAA_GBB) 736 Cmat2(iq,D2_GAB_GAB) = Cmat2(iq,D2_GAB_GAB) 737 Cmat2(iq,D2_GAB_GBB) = Cmat2(iq,D2_GAB_GBB) 738 Cmat2(iq,D2_GBB_GBB) = (t50*(-1.0d+0*t20*t30*(-1.5d+0*t72*t8 739 1 *t98+1.5d+0*t28*t52*t8/gammabb-1.5d+0*t24*t25*t8/t23**3+5 740 2 .373306513388233d-7*t101*t22*t4*t65)-2.0d+0*t20*t73*t75)+ 741 3 2.0d+0*t20*t30*t75**2*t91+1.0d+0*t20*t27*(-3.291d+0*t15*t 742 4 29*t72+5.0d-7*t101*t22*t4*t65+2.70767025d+0*gammabb*t15*t 743 5 29/t19**2.4d+1))*wght+Cmat2(iq,D2_GBB_GBB) 744 elseif (rhoa.gt.tol_rho.and.rhob.le.tol_rho) then 745 t1 = rhoa**3.333333333333333d-1 746 t2 = 1.0d+0*t1 747 t3 = t2**4.0d+0 748 t4 = param(1) 749 t5 = 5.0d-1*t4 750 t6 = gammaaa**t5 751 t7 = 1/t3**t4 752 t8 = param(2) 753 t9 = gammaaa**5.0d-1 754 t10 = 1/t3 755 t11 = asinh(t10*t9) 756 t12 = 6.0d+0*t10*t11*t8*t9+1.0746613026776465d-6*t6*t7+1.0d+ 757 1 0 758 t13 = 1/t12 759 t14 = 1/t2**8.0d+0 760 t15 = t8-1.890381166699926d-3 761 t16 = exp(-1.6455d+0*gammaaa*t14) 762 t17 = -gammaaa*t14*t8+1.0d-6*t6*t7+gammaaa*t14*t15*t16 763 t18 = t2**3.0d+0 764 t19 = 1/rhoa**6.666666666666666d-1 765 t20 = 1/t12**2 766 t21 = (gammaaa*t14+1)**5.0d-1 767 t22 = 1/t21 768 t23 = 1/t2**9.0d+0 769 t24 = -8.0d+0*gammaaa*t19*t22*t23*t8 770 t25 = 1/t2**5.0d+0 771 t26 = -8.0d+0*t11*t19*t25*t8*t9 772 t27 = 1.0d+0/t1 773 t28 = -1.4328817369035288d-6*t19*t27*t4*t6*t7+t26+t24 774 t29 = 2.6666666666666666d+0*gammaaa*t19*t23*t8 775 t30 = gammaaa**2 776 t31 = 1/t2**1.7d+1 777 t32 = 4.3879999999999997d+0*t15*t16*t19*t30*t31 778 t33 = -2.6666666666666666d+0*gammaaa*t15*t16*t19*t23 779 t34 = -1.3333333333333333d-6*t19*t27*t4*t6*t7+t33+t32+t29 780 t35 = t5-1 781 t36 = gammaaa**t35 782 t37 = 1/t2**1.6d+1 783 t38 = -t14*t8+5.0d-7*t36*t4*t7-1.6455d+0*gammaaa*t15*t16*t37 784 1 +t14*t15*t16 785 t39 = 1/t9 786 t40 = 3.0d+0*t10*t11*t39*t8+3.0d+0*t14*t22*t8+5.373306513388 787 1 233d-7*t36*t4*t7 788 t41 = 1/rhoa**1.6666666666666669d+0 789 t42 = t2**2.0d+0 790 t43 = 1/rhoa**1.3333333333333333d+0 791 t44 = 1/t12**3 792 t45 = 1/rhoa 793 t46 = -1.4328817369035288d-6*t4*t45*t6*t7+t26+t24 794 t47 = 1/rhoa**2 795 t48 = t4**2 796 t49 = 1/t2**1.0d+1 797 t50 = 1/t2**1.8d+1 798 t51 = -1.3333333333333333d-6*t4*t45*t6*t7+t33+t32+t29 799 t52 = 1/t21**3 800 t53 = gammaaa**(t5-2) 801 fnc(iq) = (1.0d+0*t13*t17*t3-9.305257363491d-1*t3)*wght+fnc( 802 1 iq) 803 Amat(iq,D1_RA) = (1.0d+0*t13*t3*t34-1.0d+0*t17*t20*t28*t3+1. 804 1 3333333333333333d+0*t13*t17*t18*t19-1.2407009817987999d+0 805 2 *t18*t19)*wght+Amat(iq,D1_RA) 806 Cmat(iq,D1_GAA) = (1.0d+0*t13*t3*t38-1.0d+0*t17*t20*t3*t40)* 807 1 wght+Cmat(iq,D1_GAA) 808 Amat2(iq,D2_RA_RA) = (t20*(-1.0d+0*t17*t3*(1.333333333333333 809 1 3d+1*t11*t43*t8*t9/t2**6.0d+0+5.333333333333333d+0*t11*t2 810 2 5*t41*t8*t9-1.0666666666666666d+1*t30*t43*t50*t52*t8+3.46 811 3 6666666666666d+1*gammaaa*t22*t43*t49*t8+5.333333333333333 812 4 d+0*gammaaa*t22*t23*t41*t8+1.9105089825380384d-6*t27*t41* 813 5 t48*t6*t7+1.4328817369035288d-6*t4*t47*t6*t7)-1.0d+0*t28* 814 6 t3*t51-1.0d+0*t3*t34*t46)+1.0d+0*t13*t3*(-8.0d+0*gammaaa* 815 7 t43*t49*t8-1.7777777777777776d+0*gammaaa*t23*t41*t8+1.777 816 8 7777777777776d-6*t27*t41*t48*t6*t7+1.3333333333333333d-6* 817 9 t4*t47*t6*t7-3.6566666666666664d+1*t15*t16*t30*t43*t50+8. 818 : 0d+0*gammaaa*t15*t16*t43*t49+1.9254543999999998d+1*gammaa 819 ; a**3*t15*t16*t43/t2**2.6d+1-2.925333333333333d+0*t15*t16* 820 < t30*t31*t41+1.7777777777777776d+0*gammaaa*t15*t16*t23*t41 821 = )+t13*t19*(1.3333333333333333d+0*t18*t51+1.33333333333333 822 > 33d+0*t18*t34)+t19*t20*(-1.3333333333333333d+0*t17*t18*t4 823 ? 6-1.3333333333333333d+0*t17*t18*t28)+2.0d+0*t17*t28*t3*t4 824 @ 4*t46+1.3333333333333333d+0*t13*t17*t42*t43-1.24070098179 825 1 87999d+0*t42*t43-8.888888888888888d-1*t13*t17*t18*t41+8.2 826 2 71339878658666d-1*t18*t41)*wght+Amat2(iq,D2_RA_RA) 827 Cmat2(iq,D2_RA_GAA) = (t20*(-1.0d+0*t17*t3*(4.0d+0*gammaaa*t 828 1 19*t31*t52*t8-4.0d+0*t11*t19*t25*t39*t8-1.2d+1*t19*t22*t2 829 2 3*t8-7.164408684517644d-7*t36*t45*t48*t7)-1.0d+0*t3*t40*t 830 3 51-1.0d+0*t3*t38*t46)+1.0d+0*t13*t3*(2.6666666666666666d+ 831 4 0*t19*t23*t8-6.666666666666666d-7*t36*t45*t48*t7+1.3164d+ 832 5 1*gammaaa*t15*t16*t19*t31-7.220454d+0*t15*t16*t19*t30/t2* 833 6 *2.5d+1-2.6666666666666666d+0*t15*t16*t19*t23)+2.0d+0*t17 834 7 *t3*t40*t44*t46-1.3333333333333333d+0*t17*t18*t19*t20*t40 835 8 +1.3333333333333333d+0*t13*t18*t19*t38)*wght+Cmat2(iq,D2_ 836 9 RA_GAA) 837 Cmat2(iq,D2_GAA_GAA) = (t20*(-1.0d+0*t17*t3*(-1.5d+0*t10*t11 838 1 *t8/t9**3-1.5d+0*t37*t52*t8+1.5d+0*t14*t22*t8/gammaaa+5.3 839 2 73306513388233d-7*t35*t4*t53*t7)-2.0d+0*t3*t38*t40)+1.0d+ 840 3 0*t13*t3*(5.0d-7*t35*t4*t53*t7-3.291d+0*t15*t16*t37+2.707 841 4 67025d+0*gammaaa*t15*t16/t2**2.4d+1)+2.0d+0*t17*t3*t40**2 842 5 *t44)*wght+Cmat2(iq,D2_GAA_GAA) 843 elseif (rhoa.le.tol_rho.and.rhob.gt.tol_rho) then 844 t1 = rhob**3.333333333333333d-1 845 t2 = 1.0d+0*t1 846 t3 = t2**4.0d+0 847 t4 = param(1) 848 t5 = 5.0d-1*t4 849 t6 = gammabb**t5 850 t7 = 1/t3**t4 851 t8 = param(2) 852 t9 = gammabb**5.0d-1 853 t10 = 1/t3 854 t11 = asinh(t10*t9) 855 t12 = 6.0d+0*t10*t11*t8*t9+1.0746613026776465d-6*t6*t7+1.0d+ 856 1 0 857 t13 = 1/t12 858 t14 = 1/t2**8.0d+0 859 t15 = t8-1.890381166699926d-3 860 t16 = exp(-1.6455d+0*gammabb*t14) 861 t17 = -gammabb*t14*t8+1.0d-6*t6*t7+gammabb*t14*t15*t16 862 t18 = t2**3.0d+0 863 t19 = 1/rhob**6.666666666666666d-1 864 t20 = 1/t12**2 865 t21 = (gammabb*t14+1)**5.0d-1 866 t22 = 1/t21 867 t23 = 1/t2**9.0d+0 868 t24 = -8.0d+0*gammabb*t19*t22*t23*t8 869 t25 = 1/t2**5.0d+0 870 t26 = -8.0d+0*t11*t19*t25*t8*t9 871 t27 = 1.0d+0/t1 872 t28 = -1.4328817369035288d-6*t19*t27*t4*t6*t7+t26+t24 873 t29 = 2.6666666666666666d+0*gammabb*t19*t23*t8 874 t30 = gammabb**2 875 t31 = 1/t2**1.7d+1 876 t32 = 4.3879999999999997d+0*t15*t16*t19*t30*t31 877 t33 = -2.6666666666666666d+0*gammabb*t15*t16*t19*t23 878 t34 = -1.3333333333333333d-6*t19*t27*t4*t6*t7+t33+t32+t29 879 t35 = t5-1 880 t36 = gammabb**t35 881 t37 = 1/t2**1.6d+1 882 t38 = -t14*t8+5.0d-7*t36*t4*t7-1.6455d+0*gammabb*t15*t16*t37 883 1 +t14*t15*t16 884 t39 = 1/t9 885 t40 = 3.0d+0*t10*t11*t39*t8+3.0d+0*t14*t22*t8+5.373306513388 886 1 233d-7*t36*t4*t7 887 t41 = 1/rhob**1.6666666666666669d+0 888 t42 = t2**2.0d+0 889 t43 = 1/rhob**1.3333333333333333d+0 890 t44 = 1/t12**3 891 t45 = 1/rhob 892 t46 = -1.4328817369035288d-6*t4*t45*t6*t7+t26+t24 893 t47 = 1/rhob**2 894 t48 = t4**2 895 t49 = 1/t2**1.0d+1 896 t50 = 1/t2**1.8d+1 897 t51 = -1.3333333333333333d-6*t4*t45*t6*t7+t33+t32+t29 898 t52 = 1/t21**3 899 t53 = gammabb**(t5-2) 900 fnc(iq) = (1.0d+0*t13*t17*t3-9.305257363491d-1*t3)*wght+fnc( 901 1 iq) 902 Amat(iq,D1_RB) = (1.0d+0*t13*t3*t34-1.0d+0*t17*t20*t28*t3+1. 903 1 3333333333333333d+0*t13*t17*t18*t19-1.2407009817987999d+0 904 2 *t18*t19)*wght+Amat(iq,D1_RB) 905 Cmat(iq,D1_GBB) = (1.0d+0*t13*t3*t38-1.0d+0*t17*t20*t3*t40)* 906 1 wght+Cmat(iq,D1_GBB) 907 Amat2(iq,D2_RB_RB) = (t20*(-1.0d+0*t17*t3*(1.333333333333333 908 1 3d+1*t11*t43*t8*t9/t2**6.0d+0+5.333333333333333d+0*t11*t2 909 2 5*t41*t8*t9-1.0666666666666666d+1*t30*t43*t50*t52*t8+3.46 910 3 6666666666666d+1*gammabb*t22*t43*t49*t8+5.333333333333333 911 4 d+0*gammabb*t22*t23*t41*t8+1.9105089825380384d-6*t27*t41* 912 5 t48*t6*t7+1.4328817369035288d-6*t4*t47*t6*t7)-1.0d+0*t28* 913 6 t3*t51-1.0d+0*t3*t34*t46)+1.0d+0*t13*t3*(-8.0d+0*gammabb* 914 7 t43*t49*t8-1.7777777777777776d+0*gammabb*t23*t41*t8+1.777 915 8 7777777777776d-6*t27*t41*t48*t6*t7+1.3333333333333333d-6* 916 9 t4*t47*t6*t7-3.6566666666666664d+1*t15*t16*t30*t43*t50+8. 917 : 0d+0*gammabb*t15*t16*t43*t49+1.9254543999999998d+1*gammab 918 ; b**3*t15*t16*t43/t2**2.6d+1-2.925333333333333d+0*t15*t16* 919 < t30*t31*t41+1.7777777777777776d+0*gammabb*t15*t16*t23*t41 920 = )+t13*t19*(1.3333333333333333d+0*t18*t51+1.33333333333333 921 > 33d+0*t18*t34)+t19*t20*(-1.3333333333333333d+0*t17*t18*t4 922 ? 6-1.3333333333333333d+0*t17*t18*t28)+2.0d+0*t17*t28*t3*t4 923 @ 4*t46+1.3333333333333333d+0*t13*t17*t42*t43-1.24070098179 924 1 87999d+0*t42*t43-8.888888888888888d-1*t13*t17*t18*t41+8.2 925 2 71339878658666d-1*t18*t41)*wght+Amat2(iq,D2_RB_RB) 926 Cmat2(iq,D2_RB_GBB) = (t20*(-1.0d+0*t17*t3*(4.0d+0*gammabb*t 927 1 19*t31*t52*t8-4.0d+0*t11*t19*t25*t39*t8-1.2d+1*t19*t22*t2 928 2 3*t8-7.164408684517644d-7*t36*t45*t48*t7)-1.0d+0*t3*t40*t 929 3 51-1.0d+0*t3*t38*t46)+1.0d+0*t13*t3*(2.6666666666666666d+ 930 4 0*t19*t23*t8-6.666666666666666d-7*t36*t45*t48*t7+1.3164d+ 931 5 1*gammabb*t15*t16*t19*t31-7.220454d+0*t15*t16*t19*t30/t2* 932 6 *2.5d+1-2.6666666666666666d+0*t15*t16*t19*t23)+2.0d+0*t17 933 7 *t3*t40*t44*t46-1.3333333333333333d+0*t17*t18*t19*t20*t40 934 8 +1.3333333333333333d+0*t13*t18*t19*t38)*wght+Cmat2(iq,D2_ 935 9 RB_GBB) 936 Cmat2(iq,D2_GBB_GBB) = (t20*(-1.0d+0*t17*t3*(-1.5d+0*t10*t11 937 1 *t8/t9**3-1.5d+0*t37*t52*t8+1.5d+0*t14*t22*t8/gammabb+5.3 938 2 73306513388233d-7*t35*t4*t53*t7)-2.0d+0*t3*t38*t40)+1.0d+ 939 3 0*t13*t3*(5.0d-7*t35*t4*t53*t7-3.291d+0*t15*t16*t37+2.707 940 4 67025d+0*gammabb*t15*t16/t2**2.4d+1)+2.0d+0*t17*t3*t40**2 941 5 *t44)*wght+Cmat2(iq,D2_GBB_GBB) 942 endif ! rhoa.gt.tol_rho.and.rhob.gt.tol_rho 943 endif ! ipol.eq.1 944 enddo ! iq 945 end 946C> 947C> \brief Evaluate the nwxcm_x_pw91 functional [1] 948C> 949C> \f{eqnarray*}{ 950C> {\it t_1} &=& 1.0\,\rho_\alpha^{{{1}\over{3}}}\\\\ 951C> {\it t_2} &=& {\it t_1}^{4.0}\\\\ 952C> {\it t_3} &=& 1.0\,\rho_\beta^{{{1}\over{3}}}\\\\ 953C> {\it t_4} &=& {\it t_3}^{4.0}\\\\ 954C> {\it t_5} &=& {\it param}\left(2\right)\\\\ 955C> {\it t_6} &=& {{1}\over{{\it t_2}}}\\\\ 956C> {\it t_7} &=& \sqrt{\sigma_{\alpha\alpha}}\\\\ 957C> {\it t_8} &=& {\it param}\left(1\right)\\\\ 958C> {\it t_9} &=& {{1}\over{{\it t_2}^{{\it t_8}}}}\\\\ 959C> {\it t_{10}} &=& {{{\it t_8}}\over{2}}\\\\ 960C> {\it t_{11}} &=& \sigma_{\alpha\alpha}^{{\it t_{10}}}\\\\ 961C> {\it t_{12}} &=& {{1}\over{{\it t_1}^{8.0}}}\\\\ 962C> {\it t_{13}} &=& {\it t_5}-0.001890381166699926\\\\ 963C> {\it t_{14}} &=& {{1}\over{{\it t_4}}}\\\\ 964C> {\it t_{15}} &=& \sqrt{\sigma_{\beta\beta}}\\\\ 965C> {\it t_{16}} &=& {{1}\over{{\it t_4}^{{\it t_8}}}}\\\\ 966C> {\it t_{17}} &=& \sigma_{\beta\beta}^{{\it t_{10}}}\\\\ 967C> {\it t_{18}} &=& {{1}\over{{\it t_3}^{8.0}}}\\\\ 968C> {\it t_{19}} &=& 1.0\,\rho_s^{{{1}\over{3}}}\\\\ 969C> {\it t_{20}} &=& {\it t_{19}}^{4.0}\\\\ 970C> {\it t_{21}} &=& {{1}\over{{\it t_{20}}}}\\\\ 971C> {\it t_{22}} &=& \sqrt{\sigma_{ss}}\\\\ 972C> {\it t_{23}} &=& {{1}\over{{\it t_{20}}^{{\it t_8}}}}\\\\ 973C> {\it t_{24}} &=& \sigma_{ss}^{{\it t_{10}}}\\\\ 974C> {\it t_{25}} &=& {{1}\over{{\it t_{19}}^{8.0}}}\\\\ 975C> f &=& {{1.0\,{\it t_4}\,\left({{{\it t_{13}}\,{\it t_{18}} 976C> \,\sigma_{\beta\beta}}\over{e^{1.6455\,{\it t_{18}}\, 977C> \sigma_{\beta\beta}}}}+1.0 \times 10^{-6}\,{\it t_{16}}\,{ 978C> \it t_{17}}-{\it t_5}\,{\it t_{18}}\, 979C> \sigma_{\beta\beta}\right)}\over{1.074661302677647 \times 10^{ 980C> -6}\,{\it t_{16}}\,{\it t_{17}}+6.0\,{\it t_5}\,{\it t_{14}} 981C> \,{\rm asinh}\; \left({\it t_{14}}\,{\it t_{15}}\right)\,{ 982C> \it t_{15}}+1.0}}+{{1.0\,{\it t_2}\,\left({{{\it t_{13}}\,{ 983C> \it t_{12}}\,\sigma_{\alpha\alpha}}\over{e^{1.6455\,{ 984C> \it t_{12}}\,\sigma_{\alpha\alpha}}}}+1.0 \times 10^{-6}\,{ 985C> \it t_9}\,{\it t_{11}}-{\it t_5}\,{\it t_{12}}\, 986C> \sigma_{\alpha\alpha}\right)}\over{1.074661302677647 \times 10^{ 987C> -6}\,{\it t_9}\,{\it t_{11}}+6.0\,{\it t_5}\,{\it t_6} 988C> \,{\rm asinh}\; \left({\it t_6}\,{\it t_7}\right)\,{\it t_7} 989C> +1.0}}-0.9305257363491\,{\it t_4}-0.9305257363491\,{\it t_2}\\\\ 990C> g &=& 0\\\\ 991C> G &=& {{1.0\,{\it t_{20}}\,\left({{{\it t_{13}}\,{\it t_{25}} 992C> \,\sigma_{ss}}\over{e^{1.6455\,{\it t_{25}}\,\sigma_{ss}}}} 993C> +1.0 \times 10^{-6}\,{\it t_{23}}\,{\it t_{24}}-{\it t_5}\,{ 994C> \it t_{25}}\,\sigma_{ss}\right)} 995C> \over{1.074661302677647 \times 10^{-6}\,{\it t_{23}}\,{ 996C> \it t_{24}}+6.0\,{\it t_5}\,{\it t_{21}}\,{\rm asinh}\; 997C> \left({\it t_{21}}\,{\it t_{22}}\right)\,{\it t_{22}}+1.0}} 998C> -0.9305257363491\,{\it t_{20}}\\\\ 999C> \f} 1000C> 1001C> Code generated with Maxima 5.34.0 [2] 1002C> driven by autoxc [3]. 1003C> 1004C> ### References ### 1005C> 1006C> [1] JP Perdew, JA Chevary, SH Vosko, KA Jackson, MR Pederson 1007C> , DJ Singh, C. Fiolhais, Phys.Rev. B 46, 6671 (1992) , DOI: 1008C> <a href="https://doi.org/10.1021/jp050536c "> 1009C> 10.1021/jp050536c </a> 1010C> 1011C> [2] Maxima, a computer algebra system, 1012C> <a href="http://maxima.sourceforge.net/"> 1013C> http://maxima.sourceforge.net/</a> 1014C> 1015C> [3] autoxc, revision 27097 2015-05-08 1016C> 1017 subroutine nwxcm_x_pw91_d3(param,tol_rho,ipol,nq,wght, 1018 +rho,rgamma,fnc,Amat,Amat2,Amat3, 1019 +Cmat,Cmat2,Cmat3) 1020c $Id: $ 1021#ifdef NWXC_QUAD_PREC 1022 implicit real(kind=selected_real_kind(30))(a-h,o-z),integer(i-n) 1023 integer, parameter :: rk=selected_real_kind(30) 1024#else 1025 implicit real(kind=selected_real_kind(15))(a-h,o-z),integer(i-n) 1026 integer, parameter :: rk=selected_real_kind(15) 1027#endif 1028#include "nwxc_param.fh" 1029 double precision param(*) !< [Input] Parameters of functional 1030 double precision tol_rho !< [Input] The lower limit on the density 1031 integer ipol !< [Input] The number of spin channels 1032 integer nq !< [Input] The number of points 1033 double precision wght !< [Input] The weight of the functional 1034 double precision rho(nq,*) !< [Input] The density 1035 double precision rgamma(nq,*) !< [Input] The norm of the density 1036 !< gradients 1037 double precision fnc(nq) !< [Output] The value of the functional 1038c 1039c Sampling Matrices for the XC Kernel 1040c 1041 double precision Amat(nq,*) !< [Output] The derivative wrt rho 1042 double precision Cmat(nq,*) !< [Output] The derivative wrt rgamma 1043c 1044c Sampling Matrices for the XC Kernel 1045c 1046 double precision Amat2(nq,*) !< [Output] The 2nd derivative wrt rho 1047 double precision Cmat2(nq,*) !< [Output] The 2nd derivative wrt rgamma 1048 !< and possibly rho 1049c 1050c Sampling Matrices for the XC Kernel 1051c 1052 double precision Amat3(nq,*) !< [Output] The 3rd derivative wrt rho 1053 double precision Cmat3(nq,*) !< [Output] The 3rd derivative wrt rgamma 1054 !< and possibly rho 1055 integer iq 1056 double precision tmp 1057 double precision rhoa,rhob 1058 double precision gammaaa,gammaab,gammabb 1059 double precision taua,taub 1060 double precision nwxcm_heaviside 1061 external nwxcm_heaviside 1062CDIR$ NOVECTOR 1063 do iq = 1, nq 1064 if (ipol.eq.1) then 1065 rhoa = 0.5d0*rho(iq,R_T) 1066 gammaaa = 0.25d0*rgamma(iq,G_TT) 1067 if (rhoa.gt.tol_rho) then 1068 t1 = 1.0d+0*rhoa**3.333333333333333d-1 1069 t2 = t1**4.0d+0 1070 t3 = param(1) 1071 t4 = 5.0d-1*t3 1072 t5 = gammaaa**t4 1073 t6 = 1/t2**t3 1074 t7 = param(2) 1075 t8 = gammaaa**5.0d-1 1076 t9 = 1/t2 1077 t10 = asinh(t8*t9) 1078 t11 = 6.0d+0*t10*t7*t8*t9+1.0746613026776465d-6*t5*t6+1.0d+0 1079 t12 = 1/t11 1080 t13 = 1/t1**8.0d+0 1081 t14 = t7-1.890381166699926d-3 1082 t15 = exp(-1.6455d+0*gammaaa*t13) 1083 t16 = -gammaaa*t13*t7+1.0d-6*t5*t6+gammaaa*t13*t14*t15 1084 t17 = t1**3.0d+0 1085 t18 = 1/rhoa**6.666666666666666d-1 1086 t19 = 1/t11**2 1087 t20 = 1/rhoa 1088 t21 = (gammaaa*t13+1)**5.0d-1 1089 t22 = 1/t21 1090 t23 = 1/t1**9.0d+0 1091 t24 = 1/t1**5.0d+0 1092 t25 = -8.0d+0*t10*t18*t24*t7*t8-8.0d+0*gammaaa*t18*t22*t23*t 1093 1 7-1.4328817369035288d-6*t20*t3*t5*t6 1094 t26 = gammaaa**2 1095 t27 = 1/t1**1.7d+1 1096 t28 = 2.6666666666666666d+0*gammaaa*t18*t23*t7-1.33333333333 1097 1 33333d-6*t20*t3*t5*t6+4.3879999999999997d+0*t14*t15*t18*t 1098 2 26*t27-2.6666666666666666d+0*gammaaa*t14*t15*t18*t23 1099 t29 = t4-1 1100 t30 = gammaaa**t29 1101 t31 = 1/t1**1.6d+1 1102 t32 = -t13*t7+5.0d-7*t3*t30*t6-1.6455d+0*gammaaa*t14*t15*t31 1103 1 +t13*t14*t15 1104 t33 = 1/t8 1105 t34 = 3.0d+0*t10*t33*t7*t9+3.0d+0*t13*t22*t7+5.3733065133882 1106 1 33d-7*t3*t30*t6 1107 t35 = 1/rhoa**1.6666666666666669d+0 1108 t36 = t1**2.0d+0 1109 t37 = 1/rhoa**1.3333333333333333d+0 1110 t38 = 1/t11**3 1111 t39 = t25**2 1112 t40 = 1/rhoa**2 1113 t41 = t3**2 1114 t42 = 1/t1**1.0d+1 1115 t43 = gammaaa**3 1116 t44 = 1/t1**2.6d+1 1117 t45 = 1/t1**1.8d+1 1118 t46 = -8.0d+0*gammaaa*t37*t42*t7-1.7777777777777776d+0*gamma 1119 1 aa*t23*t35*t7+1.7777777777777776d-6*t40*t41*t5*t6+1.33333 1120 2 33333333333d-6*t3*t40*t5*t6-3.6566666666666664d+1*t14*t15 1121 3 *t26*t37*t45+1.9254543999999998d+1*t14*t15*t37*t43*t44+8. 1122 4 0d+0*gammaaa*t14*t15*t37*t42-2.925333333333333d+0*t14*t15 1123 5 *t26*t27*t35+1.7777777777777776d+0*gammaaa*t14*t15*t23*t3 1124 6 5 1125 t47 = 1/t21**3 1126 t48 = 1/t1**6.0d+0 1127 t49 = 1.3333333333333333d+1*t10*t37*t48*t7*t8+5.333333333333 1128 1 333d+0*t10*t24*t35*t7*t8-1.0666666666666666d+1*t26*t37*t4 1129 2 5*t47*t7+3.466666666666666d+1*gammaaa*t22*t37*t42*t7+5.33 1130 3 3333333333333d+0*gammaaa*t22*t23*t35*t7+1.910508982538038 1131 4 4d-6*t40*t41*t5*t6+1.4328817369035288d-6*t3*t40*t5*t6 1132 t50 = -1.0d+0*t16*t2*t49-2.0d+0*t2*t25*t28 1133 t51 = 1/t1**2.5d+1 1134 t52 = 2.6666666666666666d+0*t18*t23*t7-6.666666666666666d-7* 1135 1 t20*t30*t41*t6-7.220454d+0*t14*t15*t18*t26*t51+1.3164d+1* 1136 2 gammaaa*t14*t15*t18*t27-2.6666666666666666d+0*t14*t15*t18 1137 3 *t23 1138 t53 = 4.0d+0*gammaaa*t18*t27*t47*t7-4.0d+0*t10*t18*t24*t33*t 1139 1 7-1.2d+1*t18*t22*t23*t7-7.164408684517644d-7*t20*t30*t41* 1140 2 t6 1141 t54 = -1.0d+0*t16*t2*t53-1.0d+0*t2*t28*t34-1.0d+0*t2*t25*t32 1142 t55 = t4-2 1143 t56 = gammaaa**t55 1144 t57 = 1/t1**2.4d+1 1145 t58 = 5.0d-7*t29*t3*t56*t6+2.70767025d+0*gammaaa*t14*t15*t57 1146 1 -3.291d+0*t14*t15*t31 1147 t59 = t34**2 1148 t60 = 1/gammaaa 1149 t61 = 1/t8**3 1150 t62 = -1.5d+0*t10*t61*t7*t9+1.5d+0*t13*t22*t60*t7-1.5d+0*t31 1151 1 *t47*t7+5.373306513388233d-7*t29*t3*t56*t6 1152 t63 = -1.0d+0*t16*t2*t62-2.0d+0*t2*t32*t34 1153 t64 = 1/rhoa**2.6666666666666666d+0 1154 t65 = 1/rhoa**2.3333333333333334d+0 1155 t66 = 1/t11**4 1156 t67 = 1/rhoa**3 1157 t68 = t3**3 1158 t69 = 1/t1**1.1d+1 1159 t70 = 1/t1**2.7d+1 1160 t71 = 1/t1**1.9d+1 1161 t72 = 1/t21**5 1162 t73 = gammaaa**(t4-3) 1163 fnc(iq) = (2.0d+0*t12*t16*t2-1.8610514726982d+0*t2)*wght+fnc 1164 1 (iq) 1165 Amat(iq,D1_RA) = (1.0d+0*t12*t2*t28-1.0d+0*t16*t19*t2*t25+1. 1166 1 3333333333333333d+0*t12*t16*t17*t18-1.2407009817987999d+0 1167 2 *t17*t18)*wght+Amat(iq,D1_RA) 1168 Cmat(iq,D1_GAA) = (1.0d+0*t12*t2*t32-1.0d+0*t16*t19*t2*t34)* 1169 1 wght+Cmat(iq,D1_GAA) 1170 Cmat(iq,D1_GAB) = Cmat(iq,D1_GAB) 1171 Amat2(iq,D2_RA_RA) = (t19*t50+1.0d+0*t12*t2*t46+2.0d+0*t16*t 1172 1 2*t38*t39+1.3333333333333333d+0*t12*t16*t36*t37-1.2407009 1173 2 817987999d+0*t36*t37-8.888888888888888d-1*t12*t16*t17*t35 1174 3 +8.271339878658666d-1*t17*t35+2.6666666666666666d+0*t12*t 1175 4 17*t18*t28-2.6666666666666666d+0*t16*t17*t18*t19*t25)*wgh 1176 5 t+Amat2(iq,D2_RA_RA) 1177 Amat2(iq,D2_RA_RB) = Amat2(iq,D2_RA_RB) 1178 Cmat2(iq,D2_RA_GAA) = (t19*t54+1.0d+0*t12*t2*t52+2.0d+0*t16* 1179 1 t2*t25*t34*t38-1.3333333333333333d+0*t16*t17*t18*t19*t34+ 1180 2 1.3333333333333333d+0*t12*t17*t18*t32)*wght+Cmat2(iq,D2_R 1181 3 A_GAA) 1182 Cmat2(iq,D2_RA_GAB) = Cmat2(iq,D2_RA_GAB) 1183 Cmat2(iq,D2_RA_GBB) = Cmat2(iq,D2_RA_GBB) 1184 Cmat2(iq,D2_GAA_GAA) = (t19*t63+2.0d+0*t16*t2*t38*t59+1.0d+0 1185 1 *t12*t2*t58)*wght+Cmat2(iq,D2_GAA_GAA) 1186 Cmat2(iq,D2_GAA_GAB) = Cmat2(iq,D2_GAA_GAB) 1187 Cmat2(iq,D2_GAA_GBB) = Cmat2(iq,D2_GAA_GBB) 1188 Cmat2(iq,D2_GAB_GAB) = Cmat2(iq,D2_GAB_GAB) 1189 Amat3(iq,D3_RA_RA_RA) = (t19*(-1.0d+0*t16*t2*(-2.66666666666 1190 1 66666d+1*t10*t48*t65*t7*t8-8.88888888888889d+0*t10*t24*t6 1191 2 4*t7*t8-2.6666666666666666d+1*t10*t40*t7*t8/t1**7.0d+0-4. 1192 3 2666666666666664d+1*t40*t43*t7*t70*t72+1.1022222222222221 1193 4 d+2*t26*t40*t47*t7*t71-1.333333333333333d+2*gammaaa*t22*t 1194 5 40*t69*t7+2.1333333333333332d+1*t26*t45*t47*t65*t7-6.9333 1195 6 33333333332d+1*gammaaa*t22*t42*t65*t7-8.88888888888889d+0 1196 7 *gammaaa*t22*t23*t64*t7-2.5473453100507176d-6*t5*t6*t67*t 1197 8 68-5.731526947614115d-6*t41*t5*t6*t67-2.8657634738070575d 1198 9 -6*t3*t5*t6*t67)-3.0d+0*t2*t28*t49-1.3333333333333333d+0* 1199 : t16*t17*t18*t49-3.0d+0*t2*t25*t46-2.6666666666666666d+0*t 1200 ; 17*t18*t25*t28)+1.0d+0*t12*t2*(2.5450399999999995d+2*t14* 1201 < t15*t26*t40*t71-3.2732724799999996d+2*t14*t15*t40*t43*t70 1202 = +2.6666666666666666d+1*gammaaa*t40*t69*t7+1.6d+1*gammaaa* 1203 > t42*t65*t7+2.962962962962963d+0*gammaaa*t23*t64*t7-2.6666 1204 ? 666666666666d+1*gammaaa*t14*t15*t40*t69-2.37037037037037d 1205 @ -6*t5*t6*t67*t68-5.333333333333333d-6*t41*t5*t6*t67-2.666 1206 1 6666666666666d-6*t3*t5*t6*t67+7.313333333333333d+1*t14*t1 1207 2 5*t26*t45*t65-3.8509087999999997d+1*t14*t15*t43*t44*t65-1 1208 3 .6d+1*gammaaa*t14*t15*t42*t65+4.875555555555556d+0*t14*t1 1209 4 5*t26*t27*t64-2.962962962962963d+0*gammaaa*t14*t15*t23*t6 1210 5 4+8.448893907199999d+1*gammaaa**4*t14*t15*t40/t1**3.5d+1) 1211 6 -6.0d+0*t16*t2*t25**3*t66-2.6666666666666666d+0*t12*t16*t 1212 7 36*t65+2.4814019635975998d+0*t36*t65+1.4814814814814814d+ 1213 8 0*t12*t16*t17*t64-1.378556646443111d+0*t17*t64+t38*(-2*t2 1214 9 5*t50+4.0d+0*t16*t2*t25*t49+2.0d+0*t2*t28*t39)+t18*t19*(- 1215 : 2.6666666666666666d+0*t16*t17*t49-2.6666666666666666d+0*t 1216 ; 16*t18*t25*t36-5.333333333333333d+0*t17*t25*t28)+t12*t18* 1217 < (4.0d+0*t17*t46+2.6666666666666666d+0*t18*t28*t36)+8.0d+0 1218 = *t16*t17*t18*t38*t39+1.3333333333333333d+0*t12*t28*t36*t3 1219 > 7-1.3333333333333333d+0*t16*t19*t25*t36*t37-2.66666666666 1220 ? 66666d+0*t12*t17*t28*t35+2.6666666666666666d+0*t16*t17*t1 1221 @ 9*t25*t35+8.888888888888888d-1*t12*t16*t35-8.271339878658 1222 1 666d-1*t35)*wght+Amat3(iq,D3_RA_RA_RA) 1223 Amat3(iq,D3_RA_RA_RB) = Amat3(iq,D3_RA_RA_RB) 1224 Cmat3(iq,D3_RA_RA_GAA) = (t19*(-1.0d+0*t16*t2*(1.6d+1*t26*t3 1225 1 7*t44*t7*t72+6.666666666666666d+0*t10*t33*t37*t48*t7-3.86 1226 2 66666666666666d+1*gammaaa*t37*t45*t47*t7-2.66666666666666 1227 3 66d+0*gammaaa*t27*t35*t47*t7+4.133333333333333d+1*t22*t37 1228 4 *t42*t7+2.6666666666666666d+0*t10*t24*t33*t35*t7+8.0d+0*t 1229 5 22*t23*t35*t7+9.552544912690191d-7*t30*t40*t6*t68+7.16440 1230 6 8684517644d-7*t30*t40*t41*t6)-2.0d+0*t2*t28*t53-2.0d+0*t2 1231 7 *t25*t52-1.0d+0*t2*t32*t49-1.0d+0*t2*t34*t46)+1.0d+0*t12* 1232 8 t2*(-8.0d+0*t37*t42*t7-1.7777777777777776d+0*t23*t35*t7+8 1233 9 .888888888888887d-7*t30*t40*t6*t68+6.666666666666666d-7*t 1234 : 30*t40*t41*t6+4.813636d+0*t14*t15*t26*t35*t51-8.629733333 1235 ; 333333d+1*gammaaa*t14*t15*t37*t45+1.17934082d+2*t14*t15*t 1236 < 26*t37*t44-3.1683352152d+1*t14*t15*t37*t43/t1**3.4d+1+8.0 1237 = d+0*t14*t15*t37*t42-8.775999999999999d+0*gammaaa*t14*t15* 1238 > t27*t35+1.7777777777777776d+0*t14*t15*t23*t35)-6.0d+0*t16 1239 ? *t2*t34*t39*t66+t38*(4.0d+0*t16*t2*t25*t53-2*t34*t50+2.0d 1240 @ +0*t2*t32*t39)+t18*t19*(-2.6666666666666666d+0*t16*t17*t5 1241 1 3-2.6666666666666666d+0*t17*t28*t34-2.6666666666666666d+0 1242 2 *t17*t25*t32)+2.6666666666666666d+0*t12*t17*t18*t52+5.333 1243 3 333333333333d+0*t16*t17*t18*t25*t34*t38-1.333333333333333 1244 4 3d+0*t16*t19*t34*t36*t37+1.3333333333333333d+0*t12*t32*t3 1245 5 6*t37+8.888888888888888d-1*t16*t17*t19*t34*t35-8.88888888 1246 6 8888888d-1*t12*t17*t32*t35)*wght+Cmat3(iq,D3_RA_RA_GAA) 1247 Cmat3(iq,D3_RA_RA_GAB) = Cmat3(iq,D3_RA_RA_GAB) 1248 Cmat3(iq,D3_RA_RA_GBB) = Cmat3(iq,D3_RA_RA_GBB) 1249 Cmat3(iq,D3_RA_RB_GAA) = Cmat3(iq,D3_RA_RB_GAA) 1250 Cmat3(iq,D3_RA_RB_GAB) = Cmat3(iq,D3_RA_RB_GAB) 1251 Cmat3(iq,D3_RA_GAA_GAA) = (t19*(-1.0d+0*t16*t2*(-6.0d+0*gamm 1252 1 aaa*t18*t51*t7*t72+2.0d+0*t10*t18*t24*t61*t7-2.0d+0*t18*t 1253 2 22*t23*t60*t7+1.0d+1*t18*t27*t47*t7-7.164408684517644d-7* 1254 3 t20*t29*t41*t56*t6)-1.0d+0*t2*t28*t62-1.0d+0*t2*t25*t58-2 1255 4 .0d+0*t2*t32*t53-2.0d+0*t2*t34*t52)-6.0d+0*t16*t2*t25*t59 1256 5 *t66+t38*(2.0d+0*t16*t2*t25*t62-2*t34*t54+2.0d+0*t16*t2*t 1257 6 34*t53+2.0d+0*t2*t25*t32*t34)+t18*t19*(-1.333333333333333 1258 7 3d+0*t16*t17*t62-2.6666666666666666d+0*t17*t32*t34)+1.0d+ 1259 8 0*t12*t2*(-6.666666666666666d-7*t20*t29*t41*t56*t6-3.6102 1260 9 27d+1*gammaaa*t14*t15*t18*t51+1.7552d+1*t14*t15*t18*t27+1 1261 : .1881257057d+1*t14*t15*t18*t26/t1**3.3d+1)+2.666666666666 1262 ; 6666d+0*t16*t17*t18*t38*t59+1.3333333333333333d+0*t12*t17 1263 < *t18*t58)*wght+Cmat3(iq,D3_RA_GAA_GAA) 1264 Cmat3(iq,D3_RA_GAA_GAB) = Cmat3(iq,D3_RA_GAA_GAB) 1265 Cmat3(iq,D3_RA_GAA_GBB) = Cmat3(iq,D3_RA_GAA_GBB) 1266 Cmat3(iq,D3_RA_GAB_GAB) = Cmat3(iq,D3_RA_GAB_GAB) 1267 Cmat3(iq,D3_RA_GAB_GBB) = Cmat3(iq,D3_RA_GAB_GBB) 1268 Cmat3(iq,D3_RA_GBB_GBB) = Cmat3(iq,D3_RA_GBB_GBB) 1269 Cmat3(iq,D3_GAA_GAA_GAA) = (t19*(-1.0d+0*t16*t2*(2.25d+0*t10 1270 1 *t7*t9/t8**5+5.373306513388233d-7*t29*t3*t55*t6*t73+2.25d 1271 2 +0*t57*t7*t72-7.5d-1*t31*t47*t60*t7-2.25d+0*t13*t22*t7/t2 1272 3 6)-3.0d+0*t2*t32*t62-3.0d+0*t2*t34*t58)+1.0d+0*t12*t2*(5. 1273 4 0d-7*t29*t3*t55*t6*t73+8.12301075d+0*t14*t15*t57-4.455471 1274 5 396375d+0*gammaaa*t14*t15/t1**3.2d+1)-6.0d+0*t16*t2*t34** 1275 6 3*t66+t38*(-2*t34*t63+4.0d+0*t16*t2*t34*t62+2.0d+0*t2*t32 1276 7 *t59))*wght+Cmat3(iq,D3_GAA_GAA_GAA) 1277 Cmat3(iq,D3_GAA_GAA_GAB) = Cmat3(iq,D3_GAA_GAA_GAB) 1278 Cmat3(iq,D3_GAA_GAA_GBB) = Cmat3(iq,D3_GAA_GAA_GBB) 1279 Cmat3(iq,D3_GAA_GAB_GAB) = Cmat3(iq,D3_GAA_GAB_GAB) 1280 Cmat3(iq,D3_GAA_GAB_GBB) = Cmat3(iq,D3_GAA_GAB_GBB) 1281 Cmat3(iq,D3_GAB_GAB_GAB) = Cmat3(iq,D3_GAB_GAB_GAB) 1282 endif ! rhoa.gt.tol_rho 1283 else ! ipol.eq.1 1284 rhoa = rho(iq,R_A) 1285 rhob = rho(iq,R_B) 1286 gammaaa = rgamma(iq,G_AA) 1287 gammaab = rgamma(iq,G_AB) 1288 gammabb = rgamma(iq,G_BB) 1289 if (rhoa.gt.tol_rho.and.rhob.gt.tol_rho) then 1290 t1 = rhoa**3.333333333333333d-1 1291 t2 = 1.0d+0*t1 1292 t3 = t2**4.0d+0 1293 t4 = param(1) 1294 t5 = 5.0d-1*t4 1295 t6 = gammaaa**t5 1296 t7 = 1/t3**t4 1297 t8 = param(2) 1298 t9 = gammaaa**5.0d-1 1299 t10 = 1/t3 1300 t11 = asinh(t10*t9) 1301 t12 = 6.0d+0*t10*t11*t8*t9+1.0746613026776465d-6*t6*t7+1.0d+ 1302 1 0 1303 t13 = 1/t12 1304 t14 = 1/t2**8.0d+0 1305 t15 = t8-1.890381166699926d-3 1306 t16 = exp(-1.6455d+0*gammaaa*t14) 1307 t17 = -gammaaa*t14*t8+1.0d-6*t6*t7+gammaaa*t14*t15*t16 1308 t18 = rhob**3.333333333333333d-1 1309 t19 = 1.0d+0*t18 1310 t20 = t19**4.0d+0 1311 t21 = gammabb**t5 1312 t22 = 1/t20**t4 1313 t23 = gammabb**5.0d-1 1314 t24 = 1/t20 1315 t25 = asinh(t23*t24) 1316 t26 = 6.0d+0*t23*t24*t25*t8+1.0746613026776465d-6*t21*t22+1. 1317 1 0d+0 1318 t27 = 1/t26 1319 t28 = 1/t19**8.0d+0 1320 t29 = exp(-1.6455d+0*gammabb*t28) 1321 t30 = -gammabb*t28*t8+gammabb*t15*t28*t29+1.0d-6*t21*t22 1322 t31 = t2**3.0d+0 1323 t32 = 1/rhoa**6.666666666666666d-1 1324 t33 = 1/t12**2 1325 t34 = (gammaaa*t14+1)**5.0d-1 1326 t35 = 1/t34 1327 t36 = 1/t2**9.0d+0 1328 t37 = -8.0d+0*gammaaa*t32*t35*t36*t8 1329 t38 = 1/t2**5.0d+0 1330 t39 = -8.0d+0*t11*t32*t38*t8*t9 1331 t40 = 1.0d+0/t1 1332 t41 = -1.4328817369035288d-6*t32*t4*t40*t6*t7+t39+t37 1333 t42 = 2.6666666666666666d+0*gammaaa*t32*t36*t8 1334 t43 = gammaaa**2 1335 t44 = 1/t2**1.7d+1 1336 t45 = 4.3879999999999997d+0*t15*t16*t32*t43*t44 1337 t46 = -2.6666666666666666d+0*gammaaa*t15*t16*t32*t36 1338 t47 = -1.3333333333333333d-6*t32*t4*t40*t6*t7+t46+t45+t42 1339 t48 = t19**3.0d+0 1340 t49 = 1/rhob**6.666666666666666d-1 1341 t50 = 1/t26**2 1342 t51 = (gammabb*t28+1)**5.0d-1 1343 t52 = 1/t51 1344 t53 = 1/t19**9.0d+0 1345 t54 = -8.0d+0*gammabb*t49*t52*t53*t8 1346 t55 = 1/t19**5.0d+0 1347 t56 = -8.0d+0*t23*t25*t49*t55*t8 1348 t57 = 1.0d+0/t18 1349 t58 = -1.4328817369035288d-6*t21*t22*t4*t49*t57+t56+t54 1350 t59 = 2.6666666666666666d+0*gammabb*t49*t53*t8 1351 t60 = gammabb**2 1352 t61 = 1/t19**1.7d+1 1353 t62 = 4.3879999999999997d+0*t15*t29*t49*t60*t61 1354 t63 = -2.6666666666666666d+0*gammabb*t15*t29*t49*t53 1355 t64 = t63+t62+t59-1.3333333333333333d-6*t21*t22*t4*t49*t57 1356 t65 = t5-1 1357 t66 = gammaaa**t65 1358 t67 = 1/t2**1.6d+1 1359 t68 = -t14*t8+5.0d-7*t4*t66*t7-1.6455d+0*gammaaa*t15*t16*t67 1360 1 +t14*t15*t16 1361 t69 = 1/t9 1362 t70 = 3.0d+0*t10*t11*t69*t8+3.0d+0*t14*t35*t8+5.373306513388 1363 1 233d-7*t4*t66*t7 1364 t71 = gammabb**t65 1365 t72 = 1/t19**1.6d+1 1366 t73 = -t28*t8-1.6455d+0*gammabb*t15*t29*t72+5.0d-7*t22*t4*t7 1367 1 1+t15*t28*t29 1368 t74 = 1/t23 1369 t75 = 3.0d+0*t24*t25*t74*t8+3.0d+0*t28*t52*t8+5.373306513388 1370 1 233d-7*t22*t4*t71 1371 t76 = 1/rhoa**1.6666666666666669d+0 1372 t77 = t2**2.0d+0 1373 t78 = 1/rhoa**1.3333333333333333d+0 1374 t79 = 1/t12**3 1375 t80 = 1/rhoa 1376 t81 = -1.4328817369035288d-6*t4*t6*t7*t80+t39+t37 1377 t82 = 1/rhoa**2 1378 t83 = 1.3333333333333333d-6*t4*t6*t7*t82 1379 t84 = -1.7777777777777776d+0*gammaaa*t36*t76*t8 1380 t85 = t4**2 1381 t86 = 1/t2**1.0d+1 1382 t87 = -8.0d+0*gammaaa*t78*t8*t86 1383 t88 = -2.925333333333333d+0*t15*t16*t43*t44*t76 1384 t89 = 1.7777777777777776d+0*gammaaa*t15*t16*t36*t76 1385 t90 = gammaaa**3 1386 t91 = 1/t2**2.6d+1 1387 t92 = 1.9254543999999998d+1*t15*t16*t78*t90*t91 1388 t93 = 1/t2**1.8d+1 1389 t94 = -3.6566666666666664d+1*t15*t16*t43*t78*t93 1390 t95 = 8.0d+0*gammaaa*t15*t16*t78*t86 1391 t96 = t95+t94+t92+t89+t88+t87+1.7777777777777776d-6*t40*t6*t 1392 1 7*t76*t85+t84+t83 1393 t97 = -1.3333333333333333d-6*t4*t6*t7*t80+t46+t45+t42 1394 t98 = 1.4328817369035288d-6*t4*t6*t7*t82 1395 t99 = 5.333333333333333d+0*gammaaa*t35*t36*t76*t8 1396 t100 = 5.333333333333333d+0*t11*t38*t76*t8*t9 1397 t101 = 1/t34**3 1398 t102 = -1.0666666666666666d+1*t101*t43*t78*t8*t93 1399 t103 = 3.466666666666666d+1*gammaaa*t35*t78*t8*t86 1400 t104 = 1/t2**6.0d+0 1401 t105 = 1.3333333333333333d+1*t104*t11*t78*t8*t9 1402 t106 = t99+t98+1.9105089825380384d-6*t40*t6*t7*t76*t85+t105+ 1403 1 t103+t102+t100 1404 t107 = 1/rhob**1.6666666666666669d+0 1405 t108 = t19**2.0d+0 1406 t109 = 1/rhob**1.3333333333333333d+0 1407 t110 = 1/t26**3 1408 t111 = 1/rhob 1409 t112 = t56+t54-1.4328817369035288d-6*t111*t21*t22*t4 1410 t113 = 1/rhob**2 1411 t114 = 1.3333333333333333d-6*t113*t21*t22*t4 1412 t115 = -1.7777777777777776d+0*gammabb*t107*t53*t8 1413 t116 = 1/t19**1.0d+1 1414 t117 = -8.0d+0*gammabb*t109*t116*t8 1415 t118 = -2.925333333333333d+0*t107*t15*t29*t60*t61 1416 t119 = 1.7777777777777776d+0*gammabb*t107*t15*t29*t53 1417 t120 = gammabb**3 1418 t121 = 1/t19**2.6d+1 1419 t122 = 1.9254543999999998d+1*t109*t120*t121*t15*t29 1420 t123 = 1/t19**1.8d+1 1421 t124 = -3.6566666666666664d+1*t109*t123*t15*t29*t60 1422 t125 = 8.0d+0*gammabb*t109*t116*t15*t29 1423 t126 = 1.7777777777777776d-6*t107*t21*t22*t57*t85+t125+t124+ 1424 1 t122+t119+t118+t117+t115+t114 1425 t127 = t63+t62+t59-1.3333333333333333d-6*t111*t21*t22*t4 1426 t128 = 1.4328817369035288d-6*t113*t21*t22*t4 1427 t129 = 5.333333333333333d+0*gammabb*t107*t52*t53*t8 1428 t130 = 5.333333333333333d+0*t107*t23*t25*t55*t8 1429 t131 = 1/t51**3 1430 t132 = -1.0666666666666666d+1*t109*t123*t131*t60*t8 1431 t133 = 3.466666666666666d+1*gammabb*t109*t116*t52*t8 1432 t134 = 1/t19**6.0d+0 1433 t135 = 1.3333333333333333d+1*t109*t134*t23*t25*t8 1434 t136 = 1.9105089825380384d-6*t107*t21*t22*t57*t85+t135+t133+ 1435 1 t132+t130+t129+t128 1436 t137 = 1/t2**2.5d+1 1437 t138 = -6.666666666666666d-7*t66*t7*t80*t85+2.66666666666666 1438 1 66d+0*t32*t36*t8+1.3164d+1*gammaaa*t15*t16*t32*t44-7.2204 1439 2 54d+0*t137*t15*t16*t32*t43-2.6666666666666666d+0*t15*t16* 1440 3 t32*t36 1441 t139 = -7.164408684517644d-7*t66*t7*t80*t85-4.0d+0*t11*t32*t 1442 1 38*t69*t8+4.0d+0*gammaaa*t101*t32*t44*t8-1.2d+1*t32*t35*t 1443 2 36*t8 1444 t140 = -1.0d+0*t3*t70*t97-1.0d+0*t3*t68*t81-1.0d+0*t139*t17* 1445 1 t3 1446 t141 = 1/t19**2.5d+1 1447 t142 = -6.666666666666666d-7*t111*t22*t71*t85+2.666666666666 1448 1 6666d+0*t49*t53*t8+1.3164d+1*gammabb*t15*t29*t49*t61-7.22 1449 2 0454d+0*t141*t15*t29*t49*t60-2.6666666666666666d+0*t15*t2 1450 3 9*t49*t53 1451 t143 = -7.164408684517644d-7*t111*t22*t71*t85-4.0d+0*t25*t49 1452 1 *t55*t74*t8+4.0d+0*gammabb*t131*t49*t61*t8-1.2d+1*t49*t52 1453 2 *t53*t8 1454 t144 = -1.0d+0*t127*t20*t75-1.0d+0*t112*t20*t73-1.0d+0*t143* 1455 1 t20*t30 1456 t145 = t5-2 1457 t146 = gammaaa**t145 1458 t147 = 1/t2**2.4d+1 1459 t148 = 5.0d-7*t146*t4*t65*t7-3.291d+0*t15*t16*t67+2.70767025 1460 1 d+0*gammaaa*t147*t15*t16 1461 t149 = t70**2 1462 t150 = 1/gammaaa 1463 t151 = 1/t9**3 1464 t152 = -1.5d+0*t101*t67*t8+1.5d+0*t14*t150*t35*t8-1.5d+0*t10 1465 1 *t11*t151*t8+5.373306513388233d-7*t146*t4*t65*t7 1466 t153 = -2.0d+0*t3*t68*t70-1.0d+0*t152*t17*t3 1467 t154 = gammabb**t145 1468 t155 = 1/t19**2.4d+1 1469 t156 = -3.291d+0*t15*t29*t72+5.0d-7*t154*t22*t4*t65+2.707670 1470 1 25d+0*gammabb*t15*t155*t29 1471 t157 = t75**2 1472 t158 = 1/gammabb 1473 t159 = 1/t23**3 1474 t160 = -1.5d+0*t131*t72*t8+1.5d+0*t158*t28*t52*t8-1.5d+0*t15 1475 1 9*t24*t25*t8+5.373306513388233d-7*t154*t22*t4*t65 1476 t161 = -2.0d+0*t20*t73*t75-1.0d+0*t160*t20*t30 1477 t162 = 1/rhoa**2.6666666666666666d+0 1478 t163 = 1/rhoa**2.3333333333333334d+0 1479 t164 = 1/t12**4 1480 t165 = t81**2 1481 t166 = 1/rhoa**3 1482 t167 = t4**3 1483 t168 = 1/t2**1.1d+1 1484 t169 = 1/t2**2.7d+1 1485 t170 = 1/t2**1.9d+1 1486 t171 = t95+t94+t92+t89+t88+t87+1.7777777777777776d-6*t6*t7*t 1487 1 82*t85+t84+t83 1488 t172 = 1/t34**5 1489 t173 = t99+t98+1.9105089825380384d-6*t6*t7*t82*t85+t105+t103 1490 1 +t102+t100 1491 t174 = -2.0d+0*t3*t81*t97-1.0d+0*t17*t173*t3 1492 t175 = 1/rhob**2.6666666666666666d+0 1493 t176 = 1/rhob**2.3333333333333334d+0 1494 t177 = 1/t26**4 1495 t178 = t112**2 1496 t179 = 1/rhob**3 1497 t180 = 1/t19**1.1d+1 1498 t181 = 1/t19**2.7d+1 1499 t182 = 1/t19**1.9d+1 1500 t183 = 1.7777777777777776d-6*t113*t21*t22*t85+t125+t124+t122 1501 1 +t119+t118+t117+t115+t114 1502 t184 = 1/t51**5 1503 t185 = 1.9105089825380384d-6*t113*t21*t22*t85+t135+t133+t132 1504 1 +t130+t129+t128 1505 t186 = -1.0d+0*t185*t20*t30-2.0d+0*t112*t127*t20 1506 t187 = t5-3 1507 t188 = gammaaa**t187 1508 t189 = gammabb**t187 1509 fnc(iq) = (1.0d+0*t20*t27*t30+1.0d+0*t13*t17*t3-9.3052573634 1510 1 91d-1*t3-9.305257363491d-1*t20)*wght+fnc(iq) 1511 Amat(iq,D1_RA) = (1.0d+0*t13*t3*t47-1.0d+0*t17*t3*t33*t41+1. 1512 1 3333333333333333d+0*t13*t17*t31*t32-1.2407009817987999d+0 1513 2 *t31*t32)*wght+Amat(iq,D1_RA) 1514 Amat(iq,D1_RB) = (1.0d+0*t20*t27*t64-1.0d+0*t20*t30*t50*t58+ 1515 1 1.3333333333333333d+0*t27*t30*t48*t49-1.2407009817987999d 1516 2 +0*t48*t49)*wght+Amat(iq,D1_RB) 1517 Cmat(iq,D1_GAA) = (1.0d+0*t13*t3*t68-1.0d+0*t17*t3*t33*t70)* 1518 1 wght+Cmat(iq,D1_GAA) 1519 Cmat(iq,D1_GAB) = Cmat(iq,D1_GAB) 1520 Cmat(iq,D1_GBB) = (1.0d+0*t20*t27*t73-1.0d+0*t20*t30*t50*t75 1521 1 )*wght+Cmat(iq,D1_GBB) 1522 Amat2(iq,D2_RA_RA) = (t33*(-1.0d+0*t3*t41*t97-1.0d+0*t3*t47* 1523 1 t81-1.0d+0*t106*t17*t3)+t13*t32*(1.3333333333333333d+0*t3 1524 2 1*t97+1.3333333333333333d+0*t31*t47)+1.0d+0*t13*t3*t96+t3 1525 3 2*t33*(-1.3333333333333333d+0*t17*t31*t81-1.3333333333333 1526 4 333d+0*t17*t31*t41)+2.0d+0*t17*t3*t41*t79*t81+1.333333333 1527 5 3333333d+0*t13*t17*t77*t78-1.2407009817987999d+0*t77*t78- 1528 6 8.888888888888888d-1*t13*t17*t31*t76+8.271339878658666d-1 1529 7 *t31*t76)*wght+Amat2(iq,D2_RA_RA) 1530 Amat2(iq,D2_RA_RB) = Amat2(iq,D2_RA_RB) 1531 Amat2(iq,D2_RB_RB) = (t27*t49*(1.3333333333333333d+0*t48*t64 1532 1 +1.3333333333333333d+0*t127*t48)+t50*(-1.0d+0*t112*t20*t6 1533 2 4-1.0d+0*t127*t20*t58-1.0d+0*t136*t20*t30)+t49*t50*(-1.33 1534 3 33333333333333d+0*t30*t48*t58-1.3333333333333333d+0*t112* 1535 4 t30*t48)+2.0d+0*t110*t112*t20*t30*t58-8.888888888888888d- 1536 5 1*t107*t27*t30*t48+8.271339878658666d-1*t107*t48+1.333333 1537 6 3333333333d+0*t108*t109*t27*t30+1.0d+0*t126*t20*t27-1.240 1538 7 7009817987999d+0*t108*t109)*wght+Amat2(iq,D2_RB_RB) 1539 Cmat2(iq,D2_RA_GAA) = (2.0d+0*t17*t3*t70*t79*t81-1.333333333 1540 1 3333333d+0*t17*t31*t32*t33*t70+1.3333333333333333d+0*t13* 1541 2 t31*t32*t68+t140*t33+1.0d+0*t13*t138*t3)*wght+Cmat2(iq,D2 1542 3 _RA_GAA) 1543 Cmat2(iq,D2_RA_GAB) = Cmat2(iq,D2_RA_GAB) 1544 Cmat2(iq,D2_RA_GBB) = Cmat2(iq,D2_RA_GBB) 1545 Cmat2(iq,D2_RB_GAA) = Cmat2(iq,D2_RB_GAA) 1546 Cmat2(iq,D2_RB_GAB) = Cmat2(iq,D2_RB_GAB) 1547 Cmat2(iq,D2_RB_GBB) = (-1.3333333333333333d+0*t30*t48*t49*t5 1548 1 0*t75+2.0d+0*t110*t112*t20*t30*t75+1.3333333333333333d+0* 1549 2 t27*t48*t49*t73+t144*t50+1.0d+0*t142*t20*t27)*wght+Cmat2( 1550 3 iq,D2_RB_GBB) 1551 Cmat2(iq,D2_GAA_GAA) = (2.0d+0*t149*t17*t3*t79+t153*t33+1.0d 1552 1 +0*t13*t148*t3)*wght+Cmat2(iq,D2_GAA_GAA) 1553 Cmat2(iq,D2_GAA_GAB) = Cmat2(iq,D2_GAA_GAB) 1554 Cmat2(iq,D2_GAA_GBB) = Cmat2(iq,D2_GAA_GBB) 1555 Cmat2(iq,D2_GAB_GAB) = Cmat2(iq,D2_GAB_GAB) 1556 Cmat2(iq,D2_GAB_GBB) = Cmat2(iq,D2_GAB_GBB) 1557 Cmat2(iq,D2_GBB_GBB) = (t161*t50+2.0d+0*t110*t157*t20*t30+1. 1558 1 0d+0*t156*t20*t27)*wght+Cmat2(iq,D2_GBB_GBB) 1559 Amat3(iq,D3_RA_RA_RA) = (t33*(-2.6666666666666666d+0*t31*t32 1560 1 *t81*t97-2.0d+0*t106*t3*t97-2.0d+0*t3*t81*t96-1.0d+0*t17* 1561 2 t3*(2.1333333333333332d+1*t101*t163*t43*t8*t93-4.26666666 1562 3 66666664d+1*t169*t172*t8*t82*t90-2.6666666666666666d+1*t1 1563 4 1*t8*t82*t9/t2**7.0d+0-8.88888888888889d+0*t11*t162*t38*t 1564 5 8*t9-2.6666666666666666d+1*t104*t11*t163*t8*t9-6.93333333 1565 6 3333332d+1*gammaaa*t163*t35*t8*t86-1.9105089825380384d-6* 1566 7 t162*t40*t6*t7*t85-3.8210179650760767d-6*t166*t6*t7*t85+1 1567 8 .1022222222222221d+2*t101*t170*t43*t8*t82-1.3333333333333 1568 9 33d+2*gammaaa*t168*t35*t8*t82-8.88888888888889d+0*gammaaa 1569 : *t162*t35*t36*t8-2.5473453100507176d-6*t162*t167*t40*t6*t 1570 ; 7-2.8657634738070575d-6*t166*t4*t6*t7)-1.0d+0*t173*t3*t47 1571 < -1.0d+0*t171*t3*t41-1.3333333333333333d+0*t17*t173*t31*t3 1572 = 2)+t13*t32*(2.6666666666666666d+0*t32*t77*t97+2.666666666 1573 > 6666666d+0*t31*t96+1.3333333333333333d+0*t171*t31)+t32*t3 1574 ? 3*(-2.6666666666666666d+0*t31*t41*t97-2.6666666666666666d 1575 @ +0*t17*t32*t77*t81-2.6666666666666666d+0*t31*t47*t81-2.66 1576 1 66666666666666d+0*t106*t17*t31)-1.7777777777777776d+0*t13 1577 2 *t31*t76*t97+1.0d+0*t13*t3*(7.313333333333333d+1*t15*t16* 1578 3 t163*t43*t93-3.8509087999999997d+1*t15*t16*t163*t90*t91-3 1579 4 .2732724799999996d+2*t15*t16*t169*t82*t90+1.6d+1*gammaaa* 1580 5 t163*t8*t86-1.6d+1*gammaaa*t15*t16*t163*t86-1.77777777777 1581 6 77776d-6*t162*t40*t6*t7*t85-3.555555555555555d-6*t166*t6* 1582 7 t7*t85+2.6666666666666666d+1*gammaaa*t168*t8*t82+2.545039 1583 8 9999999995d+2*t15*t16*t170*t43*t82+8.448893907199999d+1*g 1584 9 ammaaa**4*t15*t16*t82/t2**3.5d+1-2.6666666666666666d+1*ga 1585 : mmaaa*t15*t16*t168*t82+2.962962962962963d+0*gammaaa*t162* 1586 ; t36*t8-2.37037037037037d-6*t162*t167*t40*t6*t7-2.66666666 1587 < 66666666d-6*t166*t4*t6*t7+4.875555555555556d+0*t15*t16*t1 1588 = 62*t43*t44-2.962962962962963d+0*gammaaa*t15*t16*t162*t36) 1589 > +t32*t79*(5.333333333333333d+0*t17*t31*t41*t81+2.66666666 1590 ? 66666666d+0*t165*t17*t31)+t79*(4.0d+0*t106*t17*t3*t81+2.0 1591 @ d+0*t165*t3*t47-2*t174*t41)+1.7777777777777776d+0*t17*t31 1592 1 *t33*t76*t81+1.3333333333333333d+0*t13*t47*t77*t78-1.3333 1593 2 333333333333d+0*t17*t33*t41*t77*t78-2.6666666666666666d+0 1594 3 *t13*t163*t17*t77+2.4814019635975998d+0*t163*t77-8.888888 1595 4 888888888d-1*t13*t31*t47*t76+8.888888888888888d-1*t17*t31 1596 5 *t33*t41*t76+8.888888888888888d-1*t13*t17*t76-8.271339878 1597 6 658666d-1*t76-6.0d+0*t164*t165*t17*t3*t41+1.4814814814814 1598 7 814d+0*t13*t162*t17*t31-1.378556646443111d+0*t162*t31)*wg 1599 8 ht+Amat3(iq,D3_RA_RA_RA) 1600 Amat3(iq,D3_RA_RA_RB) = Amat3(iq,D3_RA_RA_RB) 1601 Amat3(iq,D3_RA_RB_RB) = Amat3(iq,D3_RA_RB_RB) 1602 Amat3(iq,D3_RB_RB_RB) = (t50*(-1.0d+0*t20*t30*(-1.9105089825 1603 1 380384d-6*t175*t21*t22*t57*t85-3.8210179650760767d-6*t179 1604 2 *t21*t22*t85+1.1022222222222221d+2*t113*t131*t182*t60*t8+ 1605 3 2.1333333333333332d+1*t123*t131*t176*t60*t8-8.88888888888 1606 4 889d+0*t175*t23*t25*t55*t8-8.88888888888889d+0*gammabb*t1 1607 5 75*t52*t53*t8-1.333333333333333d+2*gammabb*t113*t180*t52* 1608 6 t8-6.933333333333332d+1*gammabb*t116*t176*t52*t8-2.666666 1609 7 6666666666d+1*t113*t23*t25*t8/t19**7.0d+0-2.6666666666666 1610 8 666d+1*t134*t176*t23*t25*t8-4.2666666666666664d+1*t113*t1 1611 9 20*t181*t184*t8-2.5473453100507176d-6*t167*t175*t21*t22*t 1612 : 57-2.8657634738070575d-6*t179*t21*t22*t4)-1.0d+0*t185*t20 1613 ; *t64-1.0d+0*t183*t20*t58-1.3333333333333333d+0*t185*t30*t 1614 < 48*t49-2.6666666666666666d+0*t112*t127*t48*t49-2.0d+0*t12 1615 = 7*t136*t20-2.0d+0*t112*t126*t20)+1.0d+0*t20*t27*(-1.77777 1616 > 77777777776d-6*t175*t21*t22*t57*t85-3.555555555555555d-6* 1617 ? t179*t21*t22*t85+2.962962962962963d+0*gammabb*t175*t53*t8 1618 @ +2.6666666666666666d+1*gammabb*t113*t180*t8+1.6d+1*gammab 1619 1 b*t116*t176*t8+4.875555555555556d+0*t15*t175*t29*t60*t61+ 1620 2 2.5450399999999995d+2*t113*t15*t182*t29*t60+7.31333333333 1621 3 3333d+1*t123*t15*t176*t29*t60-2.37037037037037d-6*t167*t1 1622 4 75*t21*t22*t57-2.962962962962963d+0*gammabb*t15*t175*t29* 1623 5 t53-2.6666666666666666d-6*t179*t21*t22*t4+8.4488939071999 1624 6 99d+1*gammabb**4*t113*t15*t29/t19**3.5d+1-3.2732724799999 1625 7 996d+2*t113*t120*t15*t181*t29-2.6666666666666666d+1*gamma 1626 8 bb*t113*t15*t180*t29-3.8509087999999997d+1*t120*t121*t15* 1627 9 t176*t29-1.6d+1*gammabb*t116*t15*t176*t29)+t49*t50*(-2.66 1628 : 66666666666666d+0*t112*t48*t64-2.6666666666666666d+0*t127 1629 ; *t48*t58-2.6666666666666666d+0*t108*t112*t30*t49-2.666666 1630 < 6666666666d+0*t136*t30*t48)+t110*(2.0d+0*t178*t20*t64-2*t 1631 = 186*t58+4.0d+0*t112*t136*t20*t30)-8.888888888888888d-1*t1 1632 > 07*t27*t48*t64+1.3333333333333333d+0*t108*t109*t27*t64+t1 1633 ? 10*t49*(5.333333333333333d+0*t112*t30*t48*t58+2.666666666 1634 @ 6666666d+0*t178*t30*t48)+8.888888888888888d-1*t107*t30*t4 1635 1 8*t50*t58-1.3333333333333333d+0*t108*t109*t30*t50*t58-6.0 1636 2 d+0*t177*t178*t20*t30*t58+1.7777777777777776d+0*t107*t112 1637 3 *t30*t48*t50+t27*t49*(2.6666666666666666d+0*t108*t127*t49 1638 4 +1.3333333333333333d+0*t183*t48+2.6666666666666666d+0*t12 1639 5 6*t48)+1.4814814814814814d+0*t175*t27*t30*t48-1.777777777 1640 6 7777776d+0*t107*t127*t27*t48-1.378556646443111d+0*t175*t4 1641 7 8-2.6666666666666666d+0*t108*t176*t27*t30+8.8888888888888 1642 8 88d-1*t107*t27*t30+2.4814019635975998d+0*t108*t176-8.2713 1643 9 39878658666d-1*t107)*wght+Amat3(iq,D3_RB_RB_RB) 1644 Cmat3(iq,D3_RA_RA_GAA) = (t32*t33*(-2.6666666666666666d+0*t3 1645 1 1*t70*t97-2.6666666666666666d+0*t31*t68*t81-2.66666666666 1646 2 66666d+0*t139*t17*t31)+t33*(-2.0d+0*t139*t3*t97-1.0d+0*t1 1647 3 7*t3*(-3.8666666666666666d+1*gammaaa*t101*t78*t8*t93+1.6d 1648 4 +1*t172*t43*t78*t8*t91+4.133333333333333d+1*t35*t78*t8*t8 1649 5 6+7.164408684517644d-7*t66*t7*t82*t85+9.552544912690191d- 1650 6 7*t167*t66*t7*t82+6.666666666666666d+0*t104*t11*t69*t78*t 1651 7 8+2.6666666666666666d+0*t11*t38*t69*t76*t8-2.666666666666 1652 8 6666d+0*gammaaa*t101*t44*t76*t8+8.0d+0*t35*t36*t76*t8)-2. 1653 9 0d+0*t138*t3*t81-1.0d+0*t171*t3*t70-1.0d+0*t173*t3*t68)+1 1654 : .0d+0*t13*t3*(-8.629733333333333d+1*gammaaa*t15*t16*t78*t 1655 ; 93+1.17934082d+2*t15*t16*t43*t78*t91-3.1683352152d+1*t15* 1656 < t16*t78*t90/t2**3.4d+1-8.0d+0*t78*t8*t86+8.0d+0*t15*t16*t 1657 = 78*t86+6.666666666666666d-7*t66*t7*t82*t85+8.888888888888 1658 > 887d-7*t167*t66*t7*t82-1.7777777777777776d+0*t36*t76*t8-8 1659 ? .775999999999999d+0*gammaaa*t15*t16*t44*t76+4.813636d+0*t 1660 @ 137*t15*t16*t43*t76+1.7777777777777776d+0*t15*t16*t36*t76 1661 1 )+t79*(4.0d+0*t139*t17*t3*t81-2*t174*t70+2.0d+0*t165*t3*t 1662 2 68)+5.333333333333333d+0*t17*t31*t32*t70*t79*t81-1.333333 1663 3 3333333333d+0*t17*t33*t70*t77*t78+1.3333333333333333d+0*t 1664 4 13*t68*t77*t78+8.888888888888888d-1*t17*t31*t33*t70*t76-8 1665 5 .888888888888888d-1*t13*t31*t68*t76-6.0d+0*t164*t165*t17* 1666 6 t3*t70+2.6666666666666666d+0*t13*t138*t31*t32)*wght+Cmat3 1667 7 (iq,D3_RA_RA_GAA) 1668 Cmat3(iq,D3_RA_RA_GAB) = Cmat3(iq,D3_RA_RA_GAB) 1669 Cmat3(iq,D3_RA_RA_GBB) = Cmat3(iq,D3_RA_RA_GBB) 1670 Cmat3(iq,D3_RA_RB_GAA) = Cmat3(iq,D3_RA_RB_GAA) 1671 Cmat3(iq,D3_RA_RB_GAB) = Cmat3(iq,D3_RA_RB_GAB) 1672 Cmat3(iq,D3_RA_RB_GBB) = Cmat3(iq,D3_RA_RB_GBB) 1673 Cmat3(iq,D3_RB_RB_GAA) = Cmat3(iq,D3_RB_RB_GAA) 1674 Cmat3(iq,D3_RB_RB_GAB) = Cmat3(iq,D3_RB_RB_GAB) 1675 Cmat3(iq,D3_RB_RB_GBB) = (t50*(-1.0d+0*t20*t30*(7.1644086845 1676 1 17644d-7*t113*t22*t71*t85+2.6666666666666666d+0*t107*t25* 1677 2 t55*t74*t8+6.666666666666666d+0*t109*t134*t25*t74*t8-2.66 1678 3 66666666666666d+0*gammabb*t107*t131*t61*t8+1.6d+1*t109*t1 1679 4 21*t184*t60*t8+8.0d+0*t107*t52*t53*t8+4.133333333333333d+ 1680 5 1*t109*t116*t52*t8-3.8666666666666666d+1*gammabb*t109*t12 1681 6 3*t131*t8+9.552544912690191d-7*t113*t167*t22*t71)-1.0d+0* 1682 7 t183*t20*t75-1.0d+0*t185*t20*t73-2.0d+0*t127*t143*t20-2.0 1683 8 d+0*t112*t142*t20)+1.0d+0*t20*t27*(6.666666666666666d-7*t 1684 9 113*t22*t71*t85-1.7777777777777776d+0*t107*t53*t8-8.0d+0* 1685 : t109*t116*t8+8.888888888888887d-7*t113*t167*t22*t71-8.775 1686 ; 999999999999d+0*gammabb*t107*t15*t29*t61+4.813636d+0*t107 1687 < *t141*t15*t29*t60+1.17934082d+2*t109*t121*t15*t29*t60+1.7 1688 = 777777777777776d+0*t107*t15*t29*t53-3.1683352152d+1*t109* 1689 > t120*t15*t29/t19**3.4d+1-8.629733333333333d+1*gammabb*t10 1690 ? 9*t123*t15*t29+8.0d+0*t109*t116*t15*t29)+t49*t50*(-2.6666 1691 @ 666666666666d+0*t127*t48*t75-2.6666666666666666d+0*t112*t 1692 1 48*t73-2.6666666666666666d+0*t143*t30*t48)+t110*(-2*t186* 1693 2 t75+2.0d+0*t178*t20*t73+4.0d+0*t112*t143*t20*t30)+8.88888 1694 3 8888888888d-1*t107*t30*t48*t50*t75-1.3333333333333333d+0* 1695 4 t108*t109*t30*t50*t75+5.333333333333333d+0*t110*t112*t30* 1696 5 t48*t49*t75-6.0d+0*t177*t178*t20*t30*t75-8.88888888888888 1697 6 8d-1*t107*t27*t48*t73+1.3333333333333333d+0*t108*t109*t27 1698 7 *t73+2.6666666666666666d+0*t142*t27*t48*t49)*wght+Cmat3(i 1699 8 q,D3_RB_RB_GBB) 1700 Cmat3(iq,D3_RA_GAA_GAA) = (t33*(-1.0d+0*t152*t3*t97-1.0d+0*t 1701 1 17*t3*(-7.164408684517644d-7*t146*t65*t7*t80*t85+1.0d+1*t 1702 2 101*t32*t44*t8+2.0d+0*t11*t151*t32*t38*t8-2.0d+0*t150*t32 1703 3 *t35*t36*t8-6.0d+0*gammaaa*t137*t172*t32*t8)-1.0d+0*t148* 1704 4 t3*t81-2.0d+0*t138*t3*t70-2.0d+0*t139*t3*t68)+1.0d+0*t13* 1705 5 t3*(-6.666666666666666d-7*t146*t65*t7*t80*t85+1.7552d+1*t 1706 6 15*t16*t32*t44+1.1881257057d+1*t15*t16*t32*t43/t2**3.3d+1 1707 7 -3.610227d+1*gammaaa*t137*t15*t16*t32)+t79*(2.0d+0*t3*t68 1708 8 *t70*t81+2.0d+0*t152*t17*t3*t81+2.0d+0*t139*t17*t3*t70-2* 1709 9 t140*t70)-6.0d+0*t149*t164*t17*t3*t81+2.6666666666666666d 1710 : +0*t149*t17*t31*t32*t79+t32*t33*(-2.6666666666666666d+0*t 1711 ; 31*t68*t70-1.3333333333333333d+0*t152*t17*t31)+1.33333333 1712 < 33333333d+0*t13*t148*t31*t32)*wght+Cmat3(iq,D3_RA_GAA_GAA 1713 = ) 1714 Cmat3(iq,D3_RA_GAA_GAB) = Cmat3(iq,D3_RA_GAA_GAB) 1715 Cmat3(iq,D3_RA_GAA_GBB) = Cmat3(iq,D3_RA_GAA_GBB) 1716 Cmat3(iq,D3_RA_GAB_GAB) = Cmat3(iq,D3_RA_GAB_GAB) 1717 Cmat3(iq,D3_RA_GAB_GBB) = Cmat3(iq,D3_RA_GAB_GBB) 1718 Cmat3(iq,D3_RA_GBB_GBB) = Cmat3(iq,D3_RA_GBB_GBB) 1719 Cmat3(iq,D3_RB_GAA_GAA) = Cmat3(iq,D3_RB_GAA_GAA) 1720 Cmat3(iq,D3_RB_GAA_GAB) = Cmat3(iq,D3_RB_GAA_GAB) 1721 Cmat3(iq,D3_RB_GAA_GBB) = Cmat3(iq,D3_RB_GAA_GBB) 1722 Cmat3(iq,D3_RB_GAB_GAB) = Cmat3(iq,D3_RB_GAB_GAB) 1723 Cmat3(iq,D3_RB_GAB_GBB) = Cmat3(iq,D3_RB_GAB_GBB) 1724 Cmat3(iq,D3_RB_GBB_GBB) = (t50*(-1.0d+0*t20*t30*(-7.16440868 1725 1 4517644d-7*t111*t154*t22*t65*t85+1.0d+1*t131*t49*t61*t8+2 1726 2 .0d+0*t159*t25*t49*t55*t8-2.0d+0*t158*t49*t52*t53*t8-6.0d 1727 3 +0*gammabb*t141*t184*t49*t8)-2.0d+0*t142*t20*t75-2.0d+0*t 1728 4 143*t20*t73-1.0d+0*t127*t160*t20-1.0d+0*t112*t156*t20)+1. 1729 5 0d+0*t20*t27*(-6.666666666666666d-7*t111*t154*t22*t65*t85 1730 6 +1.7552d+1*t15*t29*t49*t61+1.1881257057d+1*t15*t29*t49*t6 1731 7 0/t19**3.3d+1-3.610227d+1*gammabb*t141*t15*t29*t49)+t49*t 1732 8 50*(-2.6666666666666666d+0*t48*t73*t75-1.3333333333333333 1733 9 d+0*t160*t30*t48)+t110*(2.0d+0*t112*t20*t73*t75+2.0d+0*t1 1734 : 43*t20*t30*t75-2*t144*t75+2.0d+0*t112*t160*t20*t30)+2.666 1735 ; 6666666666666d+0*t110*t157*t30*t48*t49+1.3333333333333333 1736 < d+0*t156*t27*t48*t49-6.0d+0*t112*t157*t177*t20*t30)*wght+ 1737 = Cmat3(iq,D3_RB_GBB_GBB) 1738 Cmat3(iq,D3_GAA_GAA_GAA) = (t33*(-1.0d+0*t17*t3*(2.25d+0*t10 1739 1 *t11*t8/t9**5-7.5d-1*t101*t150*t67*t8-2.25d+0*t14*t35*t8/ 1740 2 t43+2.25d+0*t147*t172*t8+5.373306513388233d-7*t145*t188*t 1741 3 4*t65*t7)-3.0d+0*t148*t3*t70-3.0d+0*t152*t3*t68)+(4.0d+0* 1742 4 t152*t17*t3*t70-2*t153*t70+2.0d+0*t149*t3*t68)*t79-6.0d+0 1743 5 *t164*t17*t3*t70**3+1.0d+0*t13*t3*(5.0d-7*t145*t188*t4*t6 1744 6 5*t7-4.455471396375d+0*gammaaa*t15*t16/t2**3.2d+1+8.12301 1745 7 075d+0*t147*t15*t16))*wght+Cmat3(iq,D3_GAA_GAA_GAA) 1746 Cmat3(iq,D3_GAA_GAA_GAB) = Cmat3(iq,D3_GAA_GAA_GAB) 1747 Cmat3(iq,D3_GAA_GAA_GBB) = Cmat3(iq,D3_GAA_GAA_GBB) 1748 Cmat3(iq,D3_GAA_GAB_GAB) = Cmat3(iq,D3_GAA_GAB_GAB) 1749 Cmat3(iq,D3_GAA_GAB_GBB) = Cmat3(iq,D3_GAA_GAB_GBB) 1750 Cmat3(iq,D3_GAA_GBB_GBB) = Cmat3(iq,D3_GAA_GBB_GBB) 1751 Cmat3(iq,D3_GAB_GAB_GAB) = Cmat3(iq,D3_GAB_GAB_GAB) 1752 Cmat3(iq,D3_GAB_GAB_GBB) = Cmat3(iq,D3_GAB_GAB_GBB) 1753 Cmat3(iq,D3_GAB_GBB_GBB) = Cmat3(iq,D3_GAB_GBB_GBB) 1754 Cmat3(iq,D3_GBB_GBB_GBB) = (t50*(-1.0d+0*t20*t30*(-7.5d-1*t1 1755 1 31*t158*t72*t8-2.25d+0*t28*t52*t8/t60+2.25d+0*t24*t25*t8/ 1756 2 t23**5+2.25d+0*t155*t184*t8+5.373306513388233d-7*t145*t18 1757 3 9*t22*t4*t65)-3.0d+0*t156*t20*t75-3.0d+0*t160*t20*t73)-6. 1758 4 0d+0*t177*t20*t30*t75**3+t110*(4.0d+0*t160*t20*t30*t75-2* 1759 5 t161*t75+2.0d+0*t157*t20*t73)+1.0d+0*t20*t27*(5.0d-7*t145 1760 6 *t189*t22*t4*t65-4.455471396375d+0*gammabb*t15*t29/t19**3 1761 7 .2d+1+8.12301075d+0*t15*t155*t29))*wght+Cmat3(iq,D3_GBB_G 1762 8 BB_GBB) 1763 elseif (rhoa.gt.tol_rho.and.rhob.le.tol_rho) then 1764 t1 = rhoa**3.333333333333333d-1 1765 t2 = 1.0d+0*t1 1766 t3 = t2**4.0d+0 1767 t4 = param(1) 1768 t5 = 5.0d-1*t4 1769 t6 = gammaaa**t5 1770 t7 = 1/t3**t4 1771 t8 = param(2) 1772 t9 = gammaaa**5.0d-1 1773 t10 = 1/t3 1774 t11 = asinh(t10*t9) 1775 t12 = 6.0d+0*t10*t11*t8*t9+1.0746613026776465d-6*t6*t7+1.0d+ 1776 1 0 1777 t13 = 1/t12 1778 t14 = 1/t2**8.0d+0 1779 t15 = t8-1.890381166699926d-3 1780 t16 = exp(-1.6455d+0*gammaaa*t14) 1781 t17 = -gammaaa*t14*t8+1.0d-6*t6*t7+gammaaa*t14*t15*t16 1782 t18 = t2**3.0d+0 1783 t19 = 1/rhoa**6.666666666666666d-1 1784 t20 = 1/t12**2 1785 t21 = (gammaaa*t14+1)**5.0d-1 1786 t22 = 1/t21 1787 t23 = 1/t2**9.0d+0 1788 t24 = -8.0d+0*gammaaa*t19*t22*t23*t8 1789 t25 = 1/t2**5.0d+0 1790 t26 = -8.0d+0*t11*t19*t25*t8*t9 1791 t27 = 1.0d+0/t1 1792 t28 = -1.4328817369035288d-6*t19*t27*t4*t6*t7+t26+t24 1793 t29 = 2.6666666666666666d+0*gammaaa*t19*t23*t8 1794 t30 = gammaaa**2 1795 t31 = 1/t2**1.7d+1 1796 t32 = 4.3879999999999997d+0*t15*t16*t19*t30*t31 1797 t33 = -2.6666666666666666d+0*gammaaa*t15*t16*t19*t23 1798 t34 = -1.3333333333333333d-6*t19*t27*t4*t6*t7+t33+t32+t29 1799 t35 = t5-1 1800 t36 = gammaaa**t35 1801 t37 = 1/t2**1.6d+1 1802 t38 = -t14*t8+5.0d-7*t36*t4*t7-1.6455d+0*gammaaa*t15*t16*t37 1803 1 +t14*t15*t16 1804 t39 = 1/t9 1805 t40 = 3.0d+0*t10*t11*t39*t8+3.0d+0*t14*t22*t8+5.373306513388 1806 1 233d-7*t36*t4*t7 1807 t41 = 1/rhoa**1.6666666666666669d+0 1808 t42 = t2**2.0d+0 1809 t43 = 1/rhoa**1.3333333333333333d+0 1810 t44 = 1/t12**3 1811 t45 = 1/rhoa 1812 t46 = -1.4328817369035288d-6*t4*t45*t6*t7+t26+t24 1813 t47 = 1/rhoa**2 1814 t48 = 1.3333333333333333d-6*t4*t47*t6*t7 1815 t49 = -1.7777777777777776d+0*gammaaa*t23*t41*t8 1816 t50 = t4**2 1817 t51 = 1/t2**1.0d+1 1818 t52 = -8.0d+0*gammaaa*t43*t51*t8 1819 t53 = -2.925333333333333d+0*t15*t16*t30*t31*t41 1820 t54 = 1.7777777777777776d+0*gammaaa*t15*t16*t23*t41 1821 t55 = gammaaa**3 1822 t56 = 1/t2**2.6d+1 1823 t57 = 1.9254543999999998d+1*t15*t16*t43*t55*t56 1824 t58 = 1/t2**1.8d+1 1825 t59 = -3.6566666666666664d+1*t15*t16*t30*t43*t58 1826 t60 = 8.0d+0*gammaaa*t15*t16*t43*t51 1827 t61 = 1.7777777777777776d-6*t27*t41*t50*t6*t7+t60+t59+t57+t5 1828 1 4+t53+t52+t49+t48 1829 t62 = -1.3333333333333333d-6*t4*t45*t6*t7+t33+t32+t29 1830 t63 = 1.4328817369035288d-6*t4*t47*t6*t7 1831 t64 = 5.333333333333333d+0*gammaaa*t22*t23*t41*t8 1832 t65 = 5.333333333333333d+0*t11*t25*t41*t8*t9 1833 t66 = 1/t21**3 1834 t67 = -1.0666666666666666d+1*t30*t43*t58*t66*t8 1835 t68 = 3.466666666666666d+1*gammaaa*t22*t43*t51*t8 1836 t69 = 1/t2**6.0d+0 1837 t70 = 1.3333333333333333d+1*t11*t43*t69*t8*t9 1838 t71 = t70+1.9105089825380384d-6*t27*t41*t50*t6*t7+t68+t67+t6 1839 1 5+t64+t63 1840 t72 = 1/t2**2.5d+1 1841 t73 = 2.6666666666666666d+0*t19*t23*t8-7.220454d+0*t15*t16*t 1842 1 19*t30*t72-6.666666666666666d-7*t36*t45*t50*t7+1.3164d+1* 1843 2 gammaaa*t15*t16*t19*t31-2.6666666666666666d+0*t15*t16*t19 1844 3 *t23 1845 t74 = 4.0d+0*gammaaa*t19*t31*t66*t8-4.0d+0*t11*t19*t25*t39*t 1846 1 8-1.2d+1*t19*t22*t23*t8-7.164408684517644d-7*t36*t45*t50* 1847 2 t7 1848 t75 = -1.0d+0*t17*t3*t74-1.0d+0*t3*t40*t62-1.0d+0*t3*t38*t46 1849 t76 = t5-2 1850 t77 = gammaaa**t76 1851 t78 = 1/t2**2.4d+1 1852 t79 = 2.70767025d+0*gammaaa*t15*t16*t78+5.0d-7*t35*t4*t7*t77 1853 1 -3.291d+0*t15*t16*t37 1854 t80 = t40**2 1855 t81 = 1/gammaaa 1856 t82 = 1/t9**3 1857 t83 = -1.5d+0*t10*t11*t8*t82+1.5d+0*t14*t22*t8*t81-1.5d+0*t3 1858 1 7*t66*t8+5.373306513388233d-7*t35*t4*t7*t77 1859 t84 = -1.0d+0*t17*t3*t83-2.0d+0*t3*t38*t40 1860 t85 = 1/rhoa**2.6666666666666666d+0 1861 t86 = 1/rhoa**2.3333333333333334d+0 1862 t87 = 1/t12**4 1863 t88 = t46**2 1864 t89 = 1/rhoa**3 1865 t90 = t4**3 1866 t91 = 1/t2**1.1d+1 1867 t92 = 1/t2**2.7d+1 1868 t93 = 1/t2**1.9d+1 1869 t94 = 1.7777777777777776d-6*t47*t50*t6*t7+t60+t59+t57+t54+t5 1870 1 3+t52+t49+t48 1871 t95 = 1/t21**5 1872 t96 = t70+1.9105089825380384d-6*t47*t50*t6*t7+t68+t67+t65+t6 1873 1 4+t63 1874 t97 = -1.0d+0*t17*t3*t96-2.0d+0*t3*t46*t62 1875 t98 = gammaaa**(t5-3) 1876 fnc(iq) = (1.0d+0*t13*t17*t3-9.305257363491d-1*t3)*wght+fnc( 1877 1 iq) 1878 Amat(iq,D1_RA) = (1.0d+0*t13*t3*t34-1.0d+0*t17*t20*t28*t3+1. 1879 1 3333333333333333d+0*t13*t17*t18*t19-1.2407009817987999d+0 1880 2 *t18*t19)*wght+Amat(iq,D1_RA) 1881 Cmat(iq,D1_GAA) = (1.0d+0*t13*t3*t38-1.0d+0*t17*t20*t3*t40)* 1882 1 wght+Cmat(iq,D1_GAA) 1883 Amat2(iq,D2_RA_RA) = (t20*(-1.0d+0*t17*t3*t71-1.0d+0*t28*t3* 1884 1 t62-1.0d+0*t3*t34*t46)+t13*t19*(1.3333333333333333d+0*t18 1885 2 *t62+1.3333333333333333d+0*t18*t34)+1.0d+0*t13*t3*t61+t19 1886 3 *t20*(-1.3333333333333333d+0*t17*t18*t46-1.33333333333333 1887 4 33d+0*t17*t18*t28)+2.0d+0*t17*t28*t3*t44*t46+1.3333333333 1888 5 333333d+0*t13*t17*t42*t43-1.2407009817987999d+0*t42*t43-8 1889 6 .888888888888888d-1*t13*t17*t18*t41+8.271339878658666d-1* 1890 7 t18*t41)*wght+Amat2(iq,D2_RA_RA) 1891 Cmat2(iq,D2_RA_GAA) = (t20*t75+1.0d+0*t13*t3*t73+2.0d+0*t17* 1892 1 t3*t40*t44*t46-1.3333333333333333d+0*t17*t18*t19*t20*t40+ 1893 2 1.3333333333333333d+0*t13*t18*t19*t38)*wght+Cmat2(iq,D2_R 1894 3 A_GAA) 1895 Cmat2(iq,D2_GAA_GAA) = (t20*t84+2.0d+0*t17*t3*t44*t80+1.0d+0 1896 1 *t13*t3*t79)*wght+Cmat2(iq,D2_GAA_GAA) 1897 Amat3(iq,D3_RA_RA_RA) = (t44*(-2*t28*t97+2.0d+0*t3*t34*t88+4 1898 1 .0d+0*t17*t3*t46*t71)+t20*(-1.0d+0*t3*t34*t96-1.333333333 1899 2 3333333d+0*t17*t18*t19*t96-1.0d+0*t17*t3*(-4.266666666666 1900 3 6664d+1*t47*t55*t8*t92*t95+1.1022222222222221d+2*t30*t47* 1901 4 t66*t8*t93-1.333333333333333d+2*gammaaa*t22*t47*t8*t91-2. 1902 5 5473453100507176d-6*t27*t6*t7*t85*t90-2.6666666666666666d 1903 6 +1*t11*t69*t8*t86*t9-8.88888888888889d+0*t11*t25*t8*t85*t 1904 7 9-2.6666666666666666d+1*t11*t47*t8*t9/t2**7.0d+0-3.821017 1905 8 9650760767d-6*t50*t6*t7*t89-2.8657634738070575d-6*t4*t6*t 1906 9 7*t89+2.1333333333333332d+1*t30*t58*t66*t8*t86-6.93333333 1907 : 3333332d+1*gammaaa*t22*t51*t8*t86-8.88888888888889d+0*gam 1908 ; maaa*t22*t23*t8*t85-1.9105089825380384d-6*t27*t50*t6*t7*t 1909 < 85)-1.0d+0*t28*t3*t94-2.0d+0*t3*t62*t71-2.666666666666666 1910 = 6d+0*t18*t19*t46*t62-2.0d+0*t3*t46*t61)+t13*t19*(1.333333 1911 > 3333333333d+0*t18*t94+2.6666666666666666d+0*t19*t42*t62+2 1912 ? .6666666666666666d+0*t18*t61)+1.0d+0*t13*t3*(2.5450399999 1913 @ 999995d+2*t15*t16*t30*t47*t93-3.2732724799999996d+2*t15*t 1914 1 16*t47*t55*t92+2.6666666666666666d+1*gammaaa*t47*t8*t91-2 1915 2 .6666666666666666d+1*gammaaa*t15*t16*t47*t91-2.3703703703 1916 3 7037d-6*t27*t6*t7*t85*t90-3.555555555555555d-6*t50*t6*t7* 1917 4 t89-2.6666666666666666d-6*t4*t6*t7*t89+1.6d+1*gammaaa*t51 1918 5 *t8*t86+7.313333333333333d+1*t15*t16*t30*t58*t86-3.850908 1919 6 7999999997d+1*t15*t16*t55*t56*t86-1.6d+1*gammaaa*t15*t16* 1920 7 t51*t86+2.962962962962963d+0*gammaaa*t23*t8*t85-1.7777777 1921 8 777777776d-6*t27*t50*t6*t7*t85+4.875555555555556d+0*t15*t 1922 9 16*t30*t31*t85-2.962962962962963d+0*gammaaa*t15*t16*t23*t 1923 : 85+8.448893907199999d+1*gammaaa**4*t15*t16*t47/t2**3.5d+1 1924 ; )+t19*t44*(2.6666666666666666d+0*t17*t18*t88+5.3333333333 1925 < 33333d+0*t17*t18*t28*t46)-6.0d+0*t17*t28*t3*t87*t88-2.666 1926 = 6666666666666d+0*t13*t17*t42*t86+2.4814019635975998d+0*t4 1927 > 2*t86+1.4814814814814814d+0*t13*t17*t18*t85-1.37855664644 1928 ? 3111d+0*t18*t85+t19*t20*(-2.6666666666666666d+0*t17*t18*t 1929 @ 71-2.6666666666666666d+0*t18*t28*t62-2.6666666666666666d+ 1930 1 0*t17*t19*t42*t46-2.6666666666666666d+0*t18*t34*t46)-1.77 1931 2 77777777777776d+0*t13*t18*t41*t62+1.7777777777777776d+0*t 1932 3 17*t18*t20*t41*t46+1.3333333333333333d+0*t13*t34*t42*t43- 1933 4 1.3333333333333333d+0*t17*t20*t28*t42*t43-8.8888888888888 1934 5 88d-1*t13*t18*t34*t41+8.888888888888888d-1*t17*t18*t20*t2 1935 6 8*t41+8.888888888888888d-1*t13*t17*t41-8.271339878658666d 1936 7 -1*t41)*wght+Amat3(iq,D3_RA_RA_RA) 1937 Cmat3(iq,D3_RA_RA_GAA) = (t44*(-2*t40*t97+2.0d+0*t3*t38*t88+ 1938 1 4.0d+0*t17*t3*t46*t74)+t20*(-1.0d+0*t3*t38*t96-1.0d+0*t17 1939 2 *t3*(1.6d+1*t30*t43*t56*t8*t95+9.552544912690191d-7*t36*t 1940 3 47*t7*t90+6.666666666666666d+0*t11*t39*t43*t69*t8-3.86666 1941 4 66666666666d+1*gammaaa*t43*t58*t66*t8-2.6666666666666666d 1942 5 +0*gammaaa*t31*t41*t66*t8+4.133333333333333d+1*t22*t43*t5 1943 6 1*t8+2.6666666666666666d+0*t11*t25*t39*t41*t8+8.0d+0*t22* 1944 7 t23*t41*t8+7.164408684517644d-7*t36*t47*t50*t7)-1.0d+0*t3 1945 8 *t40*t94-2.0d+0*t3*t62*t74-2.0d+0*t3*t46*t73)+1.0d+0*t13* 1946 9 t3*(8.888888888888887d-7*t36*t47*t7*t90-8.0d+0*t43*t51*t8 1947 : -1.7777777777777776d+0*t23*t41*t8+4.813636d+0*t15*t16*t30 1948 ; *t41*t72+6.666666666666666d-7*t36*t47*t50*t7-8.6297333333 1949 < 33333d+1*gammaaa*t15*t16*t43*t58+1.17934082d+2*t15*t16*t3 1950 = 0*t43*t56-3.1683352152d+1*t15*t16*t43*t55/t2**3.4d+1+8.0d 1951 > +0*t15*t16*t43*t51-8.775999999999999d+0*gammaaa*t15*t16*t 1952 ? 31*t41+1.7777777777777776d+0*t15*t16*t23*t41)-6.0d+0*t17* 1953 @ t3*t40*t87*t88+t19*t20*(-2.6666666666666666d+0*t17*t18*t7 1954 1 4-2.6666666666666666d+0*t18*t40*t62-2.6666666666666666d+0 1955 2 *t18*t38*t46)+2.6666666666666666d+0*t13*t18*t19*t73+5.333 1956 3 333333333333d+0*t17*t18*t19*t40*t44*t46-1.333333333333333 1957 4 3d+0*t17*t20*t40*t42*t43+1.3333333333333333d+0*t13*t38*t4 1958 5 2*t43+8.888888888888888d-1*t17*t18*t20*t40*t41-8.88888888 1959 6 8888888d-1*t13*t18*t38*t41)*wght+Cmat3(iq,D3_RA_RA_GAA) 1960 Cmat3(iq,D3_RA_GAA_GAA) = (t20*(-1.0d+0*t17*t3*(-6.0d+0*gamm 1961 1 aaa*t19*t72*t8*t95+2.0d+0*t11*t19*t25*t8*t82-2.0d+0*t19*t 1962 2 22*t23*t8*t81+1.0d+1*t19*t31*t66*t8-7.164408684517644d-7* 1963 3 t35*t45*t50*t7*t77)-1.0d+0*t3*t62*t83-1.0d+0*t3*t46*t79-2 1964 4 .0d+0*t3*t38*t74-2.0d+0*t3*t40*t73)-6.0d+0*t17*t3*t46*t80 1965 5 *t87+t44*(2.0d+0*t17*t3*t46*t83-2*t40*t75+2.0d+0*t17*t3*t 1966 6 40*t74+2.0d+0*t3*t38*t40*t46)+t19*t20*(-1.333333333333333 1967 7 3d+0*t17*t18*t83-2.6666666666666666d+0*t18*t38*t40)+2.666 1968 8 6666666666666d+0*t17*t18*t19*t44*t80+1.3333333333333333d+ 1969 9 0*t13*t18*t19*t79+1.0d+0*t13*t3*(-6.666666666666666d-7*t3 1970 : 5*t45*t50*t7*t77-3.610227d+1*gammaaa*t15*t16*t19*t72+1.75 1971 ; 52d+1*t15*t16*t19*t31+1.1881257057d+1*t15*t16*t19*t30/t2* 1972 < *3.3d+1))*wght+Cmat3(iq,D3_RA_GAA_GAA) 1973 Cmat3(iq,D3_GAA_GAA_GAA) = (t20*(-1.0d+0*t17*t3*(5.373306513 1974 1 388233d-7*t35*t4*t7*t76*t98+2.25d+0*t78*t8*t95+2.25d+0*t1 1975 2 0*t11*t8/t9**5-7.5d-1*t37*t66*t8*t81-2.25d+0*t14*t22*t8/t 1976 3 30)-3.0d+0*t3*t38*t83-3.0d+0*t3*t40*t79)+1.0d+0*t13*t3*(5 1977 4 .0d-7*t35*t4*t7*t76*t98+8.12301075d+0*t15*t16*t78-4.45547 1978 5 1396375d+0*gammaaa*t15*t16/t2**3.2d+1)-6.0d+0*t17*t3*t40* 1979 6 *3*t87+t44*(-2*t40*t84+4.0d+0*t17*t3*t40*t83+2.0d+0*t3*t3 1980 7 8*t80))*wght+Cmat3(iq,D3_GAA_GAA_GAA) 1981 elseif (rhoa.le.tol_rho.and.rhob.gt.tol_rho) then 1982 t1 = rhob**3.333333333333333d-1 1983 t2 = 1.0d+0*t1 1984 t3 = t2**4.0d+0 1985 t4 = param(1) 1986 t5 = 5.0d-1*t4 1987 t6 = gammabb**t5 1988 t7 = 1/t3**t4 1989 t8 = param(2) 1990 t9 = gammabb**5.0d-1 1991 t10 = 1/t3 1992 t11 = asinh(t10*t9) 1993 t12 = 6.0d+0*t10*t11*t8*t9+1.0746613026776465d-6*t6*t7+1.0d+ 1994 1 0 1995 t13 = 1/t12 1996 t14 = 1/t2**8.0d+0 1997 t15 = t8-1.890381166699926d-3 1998 t16 = exp(-1.6455d+0*gammabb*t14) 1999 t17 = -gammabb*t14*t8+1.0d-6*t6*t7+gammabb*t14*t15*t16 2000 t18 = t2**3.0d+0 2001 t19 = 1/rhob**6.666666666666666d-1 2002 t20 = 1/t12**2 2003 t21 = (gammabb*t14+1)**5.0d-1 2004 t22 = 1/t21 2005 t23 = 1/t2**9.0d+0 2006 t24 = -8.0d+0*gammabb*t19*t22*t23*t8 2007 t25 = 1/t2**5.0d+0 2008 t26 = -8.0d+0*t11*t19*t25*t8*t9 2009 t27 = 1.0d+0/t1 2010 t28 = -1.4328817369035288d-6*t19*t27*t4*t6*t7+t26+t24 2011 t29 = 2.6666666666666666d+0*gammabb*t19*t23*t8 2012 t30 = gammabb**2 2013 t31 = 1/t2**1.7d+1 2014 t32 = 4.3879999999999997d+0*t15*t16*t19*t30*t31 2015 t33 = -2.6666666666666666d+0*gammabb*t15*t16*t19*t23 2016 t34 = -1.3333333333333333d-6*t19*t27*t4*t6*t7+t33+t32+t29 2017 t35 = t5-1 2018 t36 = gammabb**t35 2019 t37 = 1/t2**1.6d+1 2020 t38 = -t14*t8+5.0d-7*t36*t4*t7-1.6455d+0*gammabb*t15*t16*t37 2021 1 +t14*t15*t16 2022 t39 = 1/t9 2023 t40 = 3.0d+0*t10*t11*t39*t8+3.0d+0*t14*t22*t8+5.373306513388 2024 1 233d-7*t36*t4*t7 2025 t41 = 1/rhob**1.6666666666666669d+0 2026 t42 = t2**2.0d+0 2027 t43 = 1/rhob**1.3333333333333333d+0 2028 t44 = 1/t12**3 2029 t45 = 1/rhob 2030 t46 = -1.4328817369035288d-6*t4*t45*t6*t7+t26+t24 2031 t47 = 1/rhob**2 2032 t48 = 1.3333333333333333d-6*t4*t47*t6*t7 2033 t49 = -1.7777777777777776d+0*gammabb*t23*t41*t8 2034 t50 = t4**2 2035 t51 = 1/t2**1.0d+1 2036 t52 = -8.0d+0*gammabb*t43*t51*t8 2037 t53 = -2.925333333333333d+0*t15*t16*t30*t31*t41 2038 t54 = 1.7777777777777776d+0*gammabb*t15*t16*t23*t41 2039 t55 = gammabb**3 2040 t56 = 1/t2**2.6d+1 2041 t57 = 1.9254543999999998d+1*t15*t16*t43*t55*t56 2042 t58 = 1/t2**1.8d+1 2043 t59 = -3.6566666666666664d+1*t15*t16*t30*t43*t58 2044 t60 = 8.0d+0*gammabb*t15*t16*t43*t51 2045 t61 = 1.7777777777777776d-6*t27*t41*t50*t6*t7+t60+t59+t57+t5 2046 1 4+t53+t52+t49+t48 2047 t62 = -1.3333333333333333d-6*t4*t45*t6*t7+t33+t32+t29 2048 t63 = 1.4328817369035288d-6*t4*t47*t6*t7 2049 t64 = 5.333333333333333d+0*gammabb*t22*t23*t41*t8 2050 t65 = 5.333333333333333d+0*t11*t25*t41*t8*t9 2051 t66 = 1/t21**3 2052 t67 = -1.0666666666666666d+1*t30*t43*t58*t66*t8 2053 t68 = 3.466666666666666d+1*gammabb*t22*t43*t51*t8 2054 t69 = 1/t2**6.0d+0 2055 t70 = 1.3333333333333333d+1*t11*t43*t69*t8*t9 2056 t71 = t70+1.9105089825380384d-6*t27*t41*t50*t6*t7+t68+t67+t6 2057 1 5+t64+t63 2058 t72 = 1/t2**2.5d+1 2059 t73 = 2.6666666666666666d+0*t19*t23*t8-7.220454d+0*t15*t16*t 2060 1 19*t30*t72-6.666666666666666d-7*t36*t45*t50*t7+1.3164d+1* 2061 2 gammabb*t15*t16*t19*t31-2.6666666666666666d+0*t15*t16*t19 2062 3 *t23 2063 t74 = 4.0d+0*gammabb*t19*t31*t66*t8-4.0d+0*t11*t19*t25*t39*t 2064 1 8-1.2d+1*t19*t22*t23*t8-7.164408684517644d-7*t36*t45*t50* 2065 2 t7 2066 t75 = -1.0d+0*t17*t3*t74-1.0d+0*t3*t40*t62-1.0d+0*t3*t38*t46 2067 t76 = t5-2 2068 t77 = gammabb**t76 2069 t78 = 1/t2**2.4d+1 2070 t79 = 2.70767025d+0*gammabb*t15*t16*t78+5.0d-7*t35*t4*t7*t77 2071 1 -3.291d+0*t15*t16*t37 2072 t80 = t40**2 2073 t81 = 1/gammabb 2074 t82 = 1/t9**3 2075 t83 = -1.5d+0*t10*t11*t8*t82+1.5d+0*t14*t22*t8*t81-1.5d+0*t3 2076 1 7*t66*t8+5.373306513388233d-7*t35*t4*t7*t77 2077 t84 = -1.0d+0*t17*t3*t83-2.0d+0*t3*t38*t40 2078 t85 = 1/rhob**2.6666666666666666d+0 2079 t86 = 1/rhob**2.3333333333333334d+0 2080 t87 = 1/t12**4 2081 t88 = t46**2 2082 t89 = 1/rhob**3 2083 t90 = t4**3 2084 t91 = 1/t2**1.1d+1 2085 t92 = 1/t2**2.7d+1 2086 t93 = 1/t2**1.9d+1 2087 t94 = 1.7777777777777776d-6*t47*t50*t6*t7+t60+t59+t57+t54+t5 2088 1 3+t52+t49+t48 2089 t95 = 1/t21**5 2090 t96 = t70+1.9105089825380384d-6*t47*t50*t6*t7+t68+t67+t65+t6 2091 1 4+t63 2092 t97 = -1.0d+0*t17*t3*t96-2.0d+0*t3*t46*t62 2093 t98 = gammabb**(t5-3) 2094 fnc(iq) = (1.0d+0*t13*t17*t3-9.305257363491d-1*t3)*wght+fnc( 2095 1 iq) 2096 Amat(iq,D1_RB) = (1.0d+0*t13*t3*t34-1.0d+0*t17*t20*t28*t3+1. 2097 1 3333333333333333d+0*t13*t17*t18*t19-1.2407009817987999d+0 2098 2 *t18*t19)*wght+Amat(iq,D1_RB) 2099 Cmat(iq,D1_GBB) = (1.0d+0*t13*t3*t38-1.0d+0*t17*t20*t3*t40)* 2100 1 wght+Cmat(iq,D1_GBB) 2101 Amat2(iq,D2_RB_RB) = (t20*(-1.0d+0*t17*t3*t71-1.0d+0*t28*t3* 2102 1 t62-1.0d+0*t3*t34*t46)+t13*t19*(1.3333333333333333d+0*t18 2103 2 *t62+1.3333333333333333d+0*t18*t34)+1.0d+0*t13*t3*t61+t19 2104 3 *t20*(-1.3333333333333333d+0*t17*t18*t46-1.33333333333333 2105 4 33d+0*t17*t18*t28)+2.0d+0*t17*t28*t3*t44*t46+1.3333333333 2106 5 333333d+0*t13*t17*t42*t43-1.2407009817987999d+0*t42*t43-8 2107 6 .888888888888888d-1*t13*t17*t18*t41+8.271339878658666d-1* 2108 7 t18*t41)*wght+Amat2(iq,D2_RB_RB) 2109 Cmat2(iq,D2_RB_GBB) = (t20*t75+1.0d+0*t13*t3*t73+2.0d+0*t17* 2110 1 t3*t40*t44*t46-1.3333333333333333d+0*t17*t18*t19*t20*t40+ 2111 2 1.3333333333333333d+0*t13*t18*t19*t38)*wght+Cmat2(iq,D2_R 2112 3 B_GBB) 2113 Cmat2(iq,D2_GBB_GBB) = (t20*t84+2.0d+0*t17*t3*t44*t80+1.0d+0 2114 1 *t13*t3*t79)*wght+Cmat2(iq,D2_GBB_GBB) 2115 Amat3(iq,D3_RB_RB_RB) = (t44*(-2*t28*t97+2.0d+0*t3*t34*t88+4 2116 1 .0d+0*t17*t3*t46*t71)+t20*(-1.0d+0*t3*t34*t96-1.333333333 2117 2 3333333d+0*t17*t18*t19*t96-1.0d+0*t17*t3*(-4.266666666666 2118 3 6664d+1*t47*t55*t8*t92*t95+1.1022222222222221d+2*t30*t47* 2119 4 t66*t8*t93-1.333333333333333d+2*gammabb*t22*t47*t8*t91-2. 2120 5 5473453100507176d-6*t27*t6*t7*t85*t90-2.6666666666666666d 2121 6 +1*t11*t69*t8*t86*t9-8.88888888888889d+0*t11*t25*t8*t85*t 2122 7 9-2.6666666666666666d+1*t11*t47*t8*t9/t2**7.0d+0-3.821017 2123 8 9650760767d-6*t50*t6*t7*t89-2.8657634738070575d-6*t4*t6*t 2124 9 7*t89+2.1333333333333332d+1*t30*t58*t66*t8*t86-6.93333333 2125 : 3333332d+1*gammabb*t22*t51*t8*t86-8.88888888888889d+0*gam 2126 ; mabb*t22*t23*t8*t85-1.9105089825380384d-6*t27*t50*t6*t7*t 2127 < 85)-1.0d+0*t28*t3*t94-2.0d+0*t3*t62*t71-2.666666666666666 2128 = 6d+0*t18*t19*t46*t62-2.0d+0*t3*t46*t61)+t13*t19*(1.333333 2129 > 3333333333d+0*t18*t94+2.6666666666666666d+0*t19*t42*t62+2 2130 ? .6666666666666666d+0*t18*t61)+1.0d+0*t13*t3*(2.5450399999 2131 @ 999995d+2*t15*t16*t30*t47*t93-3.2732724799999996d+2*t15*t 2132 1 16*t47*t55*t92+2.6666666666666666d+1*gammabb*t47*t8*t91-2 2133 2 .6666666666666666d+1*gammabb*t15*t16*t47*t91-2.3703703703 2134 3 7037d-6*t27*t6*t7*t85*t90-3.555555555555555d-6*t50*t6*t7* 2135 4 t89-2.6666666666666666d-6*t4*t6*t7*t89+1.6d+1*gammabb*t51 2136 5 *t8*t86+7.313333333333333d+1*t15*t16*t30*t58*t86-3.850908 2137 6 7999999997d+1*t15*t16*t55*t56*t86-1.6d+1*gammabb*t15*t16* 2138 7 t51*t86+2.962962962962963d+0*gammabb*t23*t8*t85-1.7777777 2139 8 777777776d-6*t27*t50*t6*t7*t85+4.875555555555556d+0*t15*t 2140 9 16*t30*t31*t85-2.962962962962963d+0*gammabb*t15*t16*t23*t 2141 : 85+8.448893907199999d+1*gammabb**4*t15*t16*t47/t2**3.5d+1 2142 ; )+t19*t44*(2.6666666666666666d+0*t17*t18*t88+5.3333333333 2143 < 33333d+0*t17*t18*t28*t46)-6.0d+0*t17*t28*t3*t87*t88-2.666 2144 = 6666666666666d+0*t13*t17*t42*t86+2.4814019635975998d+0*t4 2145 > 2*t86+1.4814814814814814d+0*t13*t17*t18*t85-1.37855664644 2146 ? 3111d+0*t18*t85+t19*t20*(-2.6666666666666666d+0*t17*t18*t 2147 @ 71-2.6666666666666666d+0*t18*t28*t62-2.6666666666666666d+ 2148 1 0*t17*t19*t42*t46-2.6666666666666666d+0*t18*t34*t46)-1.77 2149 2 77777777777776d+0*t13*t18*t41*t62+1.7777777777777776d+0*t 2150 3 17*t18*t20*t41*t46+1.3333333333333333d+0*t13*t34*t42*t43- 2151 4 1.3333333333333333d+0*t17*t20*t28*t42*t43-8.8888888888888 2152 5 88d-1*t13*t18*t34*t41+8.888888888888888d-1*t17*t18*t20*t2 2153 6 8*t41+8.888888888888888d-1*t13*t17*t41-8.271339878658666d 2154 7 -1*t41)*wght+Amat3(iq,D3_RB_RB_RB) 2155 Cmat3(iq,D3_RB_RB_GBB) = (t44*(-2*t40*t97+2.0d+0*t3*t38*t88+ 2156 1 4.0d+0*t17*t3*t46*t74)+t20*(-1.0d+0*t3*t38*t96-1.0d+0*t17 2157 2 *t3*(1.6d+1*t30*t43*t56*t8*t95+9.552544912690191d-7*t36*t 2158 3 47*t7*t90+6.666666666666666d+0*t11*t39*t43*t69*t8-3.86666 2159 4 66666666666d+1*gammabb*t43*t58*t66*t8-2.6666666666666666d 2160 5 +0*gammabb*t31*t41*t66*t8+4.133333333333333d+1*t22*t43*t5 2161 6 1*t8+2.6666666666666666d+0*t11*t25*t39*t41*t8+8.0d+0*t22* 2162 7 t23*t41*t8+7.164408684517644d-7*t36*t47*t50*t7)-1.0d+0*t3 2163 8 *t40*t94-2.0d+0*t3*t62*t74-2.0d+0*t3*t46*t73)+1.0d+0*t13* 2164 9 t3*(8.888888888888887d-7*t36*t47*t7*t90-8.0d+0*t43*t51*t8 2165 : -1.7777777777777776d+0*t23*t41*t8+4.813636d+0*t15*t16*t30 2166 ; *t41*t72+6.666666666666666d-7*t36*t47*t50*t7-8.6297333333 2167 < 33333d+1*gammabb*t15*t16*t43*t58+1.17934082d+2*t15*t16*t3 2168 = 0*t43*t56-3.1683352152d+1*t15*t16*t43*t55/t2**3.4d+1+8.0d 2169 > +0*t15*t16*t43*t51-8.775999999999999d+0*gammabb*t15*t16*t 2170 ? 31*t41+1.7777777777777776d+0*t15*t16*t23*t41)-6.0d+0*t17* 2171 @ t3*t40*t87*t88+t19*t20*(-2.6666666666666666d+0*t17*t18*t7 2172 1 4-2.6666666666666666d+0*t18*t40*t62-2.6666666666666666d+0 2173 2 *t18*t38*t46)+2.6666666666666666d+0*t13*t18*t19*t73+5.333 2174 3 333333333333d+0*t17*t18*t19*t40*t44*t46-1.333333333333333 2175 4 3d+0*t17*t20*t40*t42*t43+1.3333333333333333d+0*t13*t38*t4 2176 5 2*t43+8.888888888888888d-1*t17*t18*t20*t40*t41-8.88888888 2177 6 8888888d-1*t13*t18*t38*t41)*wght+Cmat3(iq,D3_RB_RB_GBB) 2178 Cmat3(iq,D3_RB_GBB_GBB) = (t20*(-1.0d+0*t17*t3*(-6.0d+0*gamm 2179 1 abb*t19*t72*t8*t95+2.0d+0*t11*t19*t25*t8*t82-2.0d+0*t19*t 2180 2 22*t23*t8*t81+1.0d+1*t19*t31*t66*t8-7.164408684517644d-7* 2181 3 t35*t45*t50*t7*t77)-1.0d+0*t3*t62*t83-1.0d+0*t3*t46*t79-2 2182 4 .0d+0*t3*t38*t74-2.0d+0*t3*t40*t73)-6.0d+0*t17*t3*t46*t80 2183 5 *t87+t44*(2.0d+0*t17*t3*t46*t83-2*t40*t75+2.0d+0*t17*t3*t 2184 6 40*t74+2.0d+0*t3*t38*t40*t46)+t19*t20*(-1.333333333333333 2185 7 3d+0*t17*t18*t83-2.6666666666666666d+0*t18*t38*t40)+2.666 2186 8 6666666666666d+0*t17*t18*t19*t44*t80+1.3333333333333333d+ 2187 9 0*t13*t18*t19*t79+1.0d+0*t13*t3*(-6.666666666666666d-7*t3 2188 : 5*t45*t50*t7*t77-3.610227d+1*gammabb*t15*t16*t19*t72+1.75 2189 ; 52d+1*t15*t16*t19*t31+1.1881257057d+1*t15*t16*t19*t30/t2* 2190 < *3.3d+1))*wght+Cmat3(iq,D3_RB_GBB_GBB) 2191 Cmat3(iq,D3_GBB_GBB_GBB) = (t20*(-1.0d+0*t17*t3*(5.373306513 2192 1 388233d-7*t35*t4*t7*t76*t98+2.25d+0*t78*t8*t95+2.25d+0*t1 2193 2 0*t11*t8/t9**5-7.5d-1*t37*t66*t8*t81-2.25d+0*t14*t22*t8/t 2194 3 30)-3.0d+0*t3*t38*t83-3.0d+0*t3*t40*t79)+1.0d+0*t13*t3*(5 2195 4 .0d-7*t35*t4*t7*t76*t98+8.12301075d+0*t15*t16*t78-4.45547 2196 5 1396375d+0*gammabb*t15*t16/t2**3.2d+1)-6.0d+0*t17*t3*t40* 2197 6 *3*t87+t44*(-2*t40*t84+4.0d+0*t17*t3*t40*t83+2.0d+0*t3*t3 2198 7 8*t80))*wght+Cmat3(iq,D3_GBB_GBB_GBB) 2199 endif ! rhoa.gt.tol_rho.and.rhob.gt.tol_rho 2200 endif ! ipol.eq.1 2201 enddo ! iq 2202 end 2203C> @} 2204