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