1C> \ingroup nwxc 2C> @{ 3C> 4C> \file nwxcm_x_b86b.F 5C> The nwxcm_x_b86b functional 6C> 7C> @} 8C> 9C> \ingroup nwxc_priv 10C> @{ 11C> 12C> \brief Evaluate the nwxcm_x_b86b functional [1] 13C> 14C> \f{eqnarray*}{ 15C> f &=& -{{0.00375\,\sigma_{\beta\beta}} 16C> \over{\rho_\beta^{{{4}\over{3}}}\,\left({{0.007\, 17C> \sigma_{\beta\beta}}\over{\rho_\beta^{{{8}\over{3}}}}} 18C> +1.0\right)^{0.8}}}-{{0.00375\,\sigma_{\alpha\alpha}} 19C> \over{\rho_\alpha^{{{4}\over{3}}}\,\left({{0.007\, 20C> \sigma_{\alpha\alpha}}\over{\rho_\alpha^{{{8}\over{3}}}}} 21C> +1.0\right)^{0.8}}}\\\\ 22C> g &=& 0\\\\ 23C> G &=& -{{0.00375\,\sigma_{ss}}\over{\rho_s^{{{4}\over{3}}} 24C> \,\left({{0.007\,\sigma_{ss}}\over{\rho_s^{{{8}\over{3}}}}} 25C> +1.0\right)^{0.8}}}\\\\ 26C> \f} 27C> 28C> Code generated with Maxima 5.34.0 [2] 29C> driven by autoxc [3]. 30C> 31C> ### References ### 32C> 33C> [1] AD Becke, J.Chem.Phys. 85, 7184 (1986) , DOI: 34C> <a href="https://doi.org/10.1063/1.451353 "> 35C> 10.1063/1.451353 </a> 36C> 37C> [2] Maxima, a computer algebra system, 38C> <a href="http://maxima.sourceforge.net/"> 39C> http://maxima.sourceforge.net/</a> 40C> 41C> [3] autoxc, revision 27097 2015-05-08 42C> 43 subroutine nwxcm_x_b86b(param,tol_rho,ipol,nq,wght, 44 +rho,rgamma,fnc,Amat,Cmat) 45c $Id: $ 46#ifdef NWXC_QUAD_PREC 47 implicit real(kind=selected_real_kind(30))(a-h,o-z),integer(i-n) 48 integer, parameter :: rk=selected_real_kind(30) 49#else 50 implicit real(kind=selected_real_kind(15))(a-h,o-z),integer(i-n) 51 integer, parameter :: rk=selected_real_kind(15) 52#endif 53#include "nwxc_param.fh" 54 double precision param(*) !< [Input] Parameters of functional 55 double precision tol_rho !< [Input] The lower limit on the density 56 integer ipol !< [Input] The number of spin channels 57 integer nq !< [Input] The number of points 58 double precision wght !< [Input] The weight of the functional 59 double precision rho(nq,*) !< [Input] The density 60 double precision rgamma(nq,*) !< [Input] The norm of the density 61 !< gradients 62 double precision fnc(nq) !< [Output] The value of the functional 63c 64c Sampling Matrices for the XC Kernel 65c 66 double precision Amat(nq,*) !< [Output] The derivative wrt rho 67 double precision Cmat(nq,*) !< [Output] The derivative wrt rgamma 68 integer iq 69 double precision tmp 70 double precision rhoa,rhob 71 double precision gammaaa,gammaab,gammabb 72 double precision taua,taub 73 double precision nwxcm_heaviside 74 external nwxcm_heaviside 75CDIR$ NOVECTOR 76 do iq = 1, nq 77 if (ipol.eq.1) then 78 rhoa = 0.5d0*rho(iq,R_T) 79 gammaaa = 0.25d0*rgamma(iq,G_TT) 80 if (rhoa.gt.tol_rho) then 81 t1 = 7.000000000000001d-3*gammaaa/rhoa**2.6666666666666666d+ 82 1 0+1.0d+0 83 t2 = 1/t1**8.0d-1 84 t3 = 1/rhoa**1.3333333333333333d+0 85 t4 = 1/t1**1.8d+0 86 fnc(iq) = fnc(iq)-7.5d-3*gammaaa*t2*t3*wght 87 Amat(iq,D1_RA) = (4.9999999999999994d-3*gammaaa*t2/rhoa**2.3 88 1 333333333333334d+0-5.599999999999999d-5*gammaaa**2*t4/rho 89 2 a**5)*wght+Amat(iq,D1_RA) 90 Cmat(iq,D1_GAA) = (2.1000000000000002d-5*gammaaa*t4/rhoa**4- 91 1 3.75d-3*t2*t3)*wght+Cmat(iq,D1_GAA) 92 Cmat(iq,D1_GAB) = Cmat(iq,D1_GAB) 93 endif ! rhoa.gt.tol_rho 94 else ! ipol.eq.1 95 rhoa = rho(iq,R_A) 96 rhob = rho(iq,R_B) 97 gammaaa = rgamma(iq,G_AA) 98 gammaab = rgamma(iq,G_AB) 99 gammabb = rgamma(iq,G_BB) 100 if (rhoa.gt.tol_rho.and.rhob.gt.tol_rho) then 101 t1 = 7.000000000000001d-3*gammaaa/rhoa**2.6666666666666666d+ 102 1 0+1.0d+0 103 t2 = 1/t1**8.0d-1 104 t3 = 1/rhoa**1.3333333333333333d+0 105 t4 = 7.000000000000001d-3*gammabb/rhob**2.6666666666666666d+ 106 1 0+1.0d+0 107 t5 = 1/t4**8.0d-1 108 t6 = 1/rhob**1.3333333333333333d+0 109 t7 = 1/t1**1.8d+0 110 t8 = 1/t4**1.8d+0 111 fnc(iq) = (-3.75d-3*gammabb*t5*t6-3.75d-3*gammaaa*t2*t3)*wgh 112 1 t+fnc(iq) 113 Amat(iq,D1_RA) = (4.9999999999999994d-3*gammaaa*t2/rhoa**2.3 114 1 333333333333334d+0-5.599999999999999d-5*gammaaa**2*t7/rho 115 2 a**5)*wght+Amat(iq,D1_RA) 116 Amat(iq,D1_RB) = (4.9999999999999994d-3*gammabb*t5/rhob**2.3 117 1 333333333333334d+0-5.599999999999999d-5*gammabb**2*t8/rho 118 2 b**5)*wght+Amat(iq,D1_RB) 119 Cmat(iq,D1_GAA) = (2.1000000000000002d-5*gammaaa*t7/rhoa**4- 120 1 3.75d-3*t2*t3)*wght+Cmat(iq,D1_GAA) 121 Cmat(iq,D1_GAB) = Cmat(iq,D1_GAB) 122 Cmat(iq,D1_GBB) = (2.1000000000000002d-5*gammabb*t8/rhob**4- 123 1 3.75d-3*t5*t6)*wght+Cmat(iq,D1_GBB) 124 elseif (rhoa.gt.tol_rho.and.rhob.le.tol_rho) then 125 t1 = 7.000000000000001d-3*gammaaa/rhoa**2.6666666666666666d+ 126 1 0+1.0d+0 127 t2 = 1/t1**8.0d-1 128 t3 = 1/rhoa**1.3333333333333333d+0 129 t4 = 1/t1**1.8d+0 130 fnc(iq) = fnc(iq)-3.75d-3*gammaaa*t2*t3*wght 131 Amat(iq,D1_RA) = -5.599999999999999d-5*gammaaa**2*t4*wght/rh 132 1 oa**5+4.9999999999999994d-3*gammaaa*t2*wght/rhoa**2.33333 133 2 33333333334d+0+Amat(iq,D1_RA) 134 Cmat(iq,D1_GAA) = 2.1000000000000002d-5*gammaaa*t4*wght/rhoa 135 1 **4-3.75d-3*t2*t3*wght+Cmat(iq,D1_GAA) 136 elseif (rhoa.le.tol_rho.and.rhob.gt.tol_rho) then 137 t1 = 7.000000000000001d-3*gammabb/rhob**2.6666666666666666d+ 138 1 0+1.0d+0 139 t2 = 1/t1**8.0d-1 140 t3 = 1/rhob**1.3333333333333333d+0 141 t4 = 1/t1**1.8d+0 142 fnc(iq) = fnc(iq)-3.75d-3*gammabb*t2*t3*wght 143 Amat(iq,D1_RB) = -5.599999999999999d-5*gammabb**2*t4*wght/rh 144 1 ob**5+4.9999999999999994d-3*gammabb*t2*wght/rhob**2.33333 145 2 33333333334d+0+Amat(iq,D1_RB) 146 Cmat(iq,D1_GBB) = 2.1000000000000002d-5*gammabb*t4*wght/rhob 147 1 **4-3.75d-3*t2*t3*wght+Cmat(iq,D1_GBB) 148 endif ! rhoa.gt.tol_rho.and.rhob.gt.tol_rho 149 endif ! ipol.eq.1 150 enddo ! iq 151 end 152C> 153C> \brief Evaluate the nwxcm_x_b86b functional [1] 154C> 155C> \f{eqnarray*}{ 156C> f &=& -{{0.00375\,\sigma_{\beta\beta}} 157C> \over{\rho_\beta^{{{4}\over{3}}}\,\left({{0.007\, 158C> \sigma_{\beta\beta}}\over{\rho_\beta^{{{8}\over{3}}}}} 159C> +1.0\right)^{0.8}}}-{{0.00375\,\sigma_{\alpha\alpha}} 160C> \over{\rho_\alpha^{{{4}\over{3}}}\,\left({{0.007\, 161C> \sigma_{\alpha\alpha}}\over{\rho_\alpha^{{{8}\over{3}}}}} 162C> +1.0\right)^{0.8}}}\\\\ 163C> g &=& 0\\\\ 164C> G &=& -{{0.00375\,\sigma_{ss}}\over{\rho_s^{{{4}\over{3}}} 165C> \,\left({{0.007\,\sigma_{ss}}\over{\rho_s^{{{8}\over{3}}}}} 166C> +1.0\right)^{0.8}}}\\\\ 167C> \f} 168C> 169C> Code generated with Maxima 5.34.0 [2] 170C> driven by autoxc [3]. 171C> 172C> ### References ### 173C> 174C> [1] AD Becke, J.Chem.Phys. 85, 7184 (1986) , DOI: 175C> <a href="https://doi.org/10.1063/1.451353 "> 176C> 10.1063/1.451353 </a> 177C> 178C> [2] Maxima, a computer algebra system, 179C> <a href="http://maxima.sourceforge.net/"> 180C> http://maxima.sourceforge.net/</a> 181C> 182C> [3] autoxc, revision 27097 2015-05-08 183C> 184 subroutine nwxcm_x_b86b_d2(param,tol_rho,ipol,nq,wght, 185 +rho,rgamma,fnc,Amat,Amat2,Cmat,Cmat2) 186c $Id: $ 187#ifdef NWXC_QUAD_PREC 188 implicit real(kind=selected_real_kind(30))(a-h,o-z),integer(i-n) 189 integer, parameter :: rk=selected_real_kind(30) 190#else 191 implicit real(kind=selected_real_kind(15))(a-h,o-z),integer(i-n) 192 integer, parameter :: rk=selected_real_kind(15) 193#endif 194#include "nwxc_param.fh" 195 double precision param(*) !< [Input] Parameters of functional 196 double precision tol_rho !< [Input] The lower limit on the density 197 integer ipol !< [Input] The number of spin channels 198 integer nq !< [Input] The number of points 199 double precision wght !< [Input] The weight of the functional 200 double precision rho(nq,*) !< [Input] The density 201 double precision rgamma(nq,*) !< [Input] The norm of the density 202 !< gradients 203 double precision fnc(nq) !< [Output] The value of the functional 204c 205c Sampling Matrices for the XC Kernel 206c 207 double precision Amat(nq,*) !< [Output] The derivative wrt rho 208 double precision Cmat(nq,*) !< [Output] The derivative wrt rgamma 209c 210c Sampling Matrices for the XC Kernel 211c 212 double precision Amat2(nq,*) !< [Output] The 2nd derivative wrt rho 213 double precision Cmat2(nq,*) !< [Output] The 2nd derivative wrt rgamma 214 !< and possibly rho 215 integer iq 216 double precision tmp 217 double precision rhoa,rhob 218 double precision gammaaa,gammaab,gammabb 219 double precision taua,taub 220 double precision nwxcm_heaviside 221 external nwxcm_heaviside 222CDIR$ NOVECTOR 223 do iq = 1, nq 224 if (ipol.eq.1) then 225 rhoa = 0.5d0*rho(iq,R_T) 226 gammaaa = 0.25d0*rgamma(iq,G_TT) 227 if (rhoa.gt.tol_rho) then 228 t1 = 7.000000000000001d-3*gammaaa/rhoa**2.6666666666666666d+ 229 1 0+1.0d+0 230 t2 = 1/t1**8.0d-1 231 t3 = 1/rhoa**1.3333333333333333d+0 232 t4 = gammaaa**2 233 t5 = 1/t1**1.8d+0 234 t6 = 1/rhoa**5 235 t7 = 1/rhoa**2.3333333333333334d+0 236 t8 = 1/rhoa**4 237 t9 = 1/t1**2.7999999999999997d+0 238 fnc(iq) = fnc(iq)-7.5d-3*gammaaa*t2*t3*wght 239 Amat(iq,D1_RA) = (4.9999999999999994d-3*gammaaa*t2*t7-5.5999 240 1 99999999999d-5*t4*t5*t6)*wght+Amat(iq,D1_RA) 241 Cmat(iq,D1_GAA) = (2.1000000000000002d-5*gammaaa*t5*t8-3.75d 242 1 -3*t2*t3)*wght+Cmat(iq,D1_GAA) 243 Cmat(iq,D1_GAB) = Cmat(iq,D1_GAB) 244 Amat2(iq,D2_RA_RA) = (-1.8816d-6*gammaaa**3*t9/rhoa**8.66666 245 1 6666666666d+0+3.5466666666666663d-4*t4*t5/rhoa**6-1.16666 246 2 66666666665d-2*gammaaa*t2/rhoa**3.3333333333333337d+0)*wg 247 3 ht+Amat2(iq,D2_RA_RA) 248 Amat2(iq,D2_RA_RB) = Amat2(iq,D2_RA_RB) 249 Cmat2(iq,D2_RA_GAA) = (7.056d-7*t4*t9/rhoa**7.66666666666666 250 1 7d+0+4.9999999999999994d-3*t2*t7-1.3999999999999999d-4*ga 251 2 mmaaa*t5*t6)*wght+Cmat2(iq,D2_RA_GAA) 252 Cmat2(iq,D2_RA_GAB) = Cmat2(iq,D2_RA_GAB) 253 Cmat2(iq,D2_RA_GBB) = Cmat2(iq,D2_RA_GBB) 254 Cmat2(iq,D2_GAA_GAA) = (4.2000000000000004d-5*t5*t8-2.646d-7 255 1 *gammaaa*t9/rhoa**6.666666666666667d+0)*wght+Cmat2(iq,D2_ 256 2 GAA_GAA) 257 Cmat2(iq,D2_GAA_GAB) = Cmat2(iq,D2_GAA_GAB) 258 Cmat2(iq,D2_GAA_GBB) = Cmat2(iq,D2_GAA_GBB) 259 Cmat2(iq,D2_GAB_GAB) = Cmat2(iq,D2_GAB_GAB) 260 endif ! rhoa.gt.tol_rho 261 else ! ipol.eq.1 262 rhoa = rho(iq,R_A) 263 rhob = rho(iq,R_B) 264 gammaaa = rgamma(iq,G_AA) 265 gammaab = rgamma(iq,G_AB) 266 gammabb = rgamma(iq,G_BB) 267 if (rhoa.gt.tol_rho.and.rhob.gt.tol_rho) then 268 t1 = 7.000000000000001d-3*gammaaa/rhoa**2.6666666666666666d+ 269 1 0+1.0d+0 270 t2 = 1/t1**8.0d-1 271 t3 = 1/rhoa**1.3333333333333333d+0 272 t4 = 7.000000000000001d-3*gammabb/rhob**2.6666666666666666d+ 273 1 0+1.0d+0 274 t5 = 1/t4**8.0d-1 275 t6 = 1/rhob**1.3333333333333333d+0 276 t7 = gammaaa**2 277 t8 = 1/t1**1.8d+0 278 t9 = 1/rhoa**5 279 t10 = 1/rhoa**2.3333333333333334d+0 280 t11 = gammabb**2 281 t12 = 1/t4**1.8d+0 282 t13 = 1/rhob**5 283 t14 = 1/rhob**2.3333333333333334d+0 284 t15 = 1/rhoa**4 285 t16 = 1/rhob**4 286 t17 = 1/t1**2.7999999999999997d+0 287 t18 = 1/t4**2.7999999999999997d+0 288 fnc(iq) = (-3.75d-3*gammabb*t5*t6-3.75d-3*gammaaa*t2*t3)*wgh 289 1 t+fnc(iq) 290 Amat(iq,D1_RA) = (4.9999999999999994d-3*gammaaa*t10*t2-5.599 291 1 999999999999d-5*t7*t8*t9)*wght+Amat(iq,D1_RA) 292 Amat(iq,D1_RB) = (4.9999999999999994d-3*gammabb*t14*t5-5.599 293 1 999999999999d-5*t11*t12*t13)*wght+Amat(iq,D1_RB) 294 Cmat(iq,D1_GAA) = (2.1000000000000002d-5*gammaaa*t15*t8-3.75 295 1 d-3*t2*t3)*wght+Cmat(iq,D1_GAA) 296 Cmat(iq,D1_GAB) = Cmat(iq,D1_GAB) 297 Cmat(iq,D1_GBB) = (2.1000000000000002d-5*gammabb*t12*t16-3.7 298 1 5d-3*t5*t6)*wght+Cmat(iq,D1_GBB) 299 Amat2(iq,D2_RA_RA) = (3.5466666666666663d-4*t7*t8/rhoa**6-1. 300 1 1666666666666665d-2*gammaaa*t2/rhoa**3.3333333333333337d+ 301 2 0-1.8816d-6*gammaaa**3*t17/rhoa**8.666666666666666d+0)*wg 302 3 ht+Amat2(iq,D2_RA_RA) 303 Amat2(iq,D2_RA_RB) = Amat2(iq,D2_RA_RB) 304 Amat2(iq,D2_RB_RB) = (-1.1666666666666665d-2*gammabb*t5/rhob 305 1 **3.3333333333333337d+0-1.8816d-6*gammabb**3*t18/rhob**8. 306 2 666666666666666d+0+3.5466666666666663d-4*t11*t12/rhob**6) 307 3 *wght+Amat2(iq,D2_RB_RB) 308 Cmat2(iq,D2_RA_GAA) = (-1.3999999999999999d-4*gammaaa*t8*t9+ 309 1 7.056d-7*t17*t7/rhoa**7.666666666666667d+0+4.999999999999 310 2 9994d-3*t10*t2)*wght+Cmat2(iq,D2_RA_GAA) 311 Cmat2(iq,D2_RA_GAB) = Cmat2(iq,D2_RA_GAB) 312 Cmat2(iq,D2_RA_GBB) = Cmat2(iq,D2_RA_GBB) 313 Cmat2(iq,D2_RB_GAA) = Cmat2(iq,D2_RB_GAA) 314 Cmat2(iq,D2_RB_GAB) = Cmat2(iq,D2_RB_GAB) 315 Cmat2(iq,D2_RB_GBB) = (4.9999999999999994d-3*t14*t5+7.056d-7 316 1 *t11*t18/rhob**7.666666666666667d+0-1.3999999999999999d-4 317 2 *gammabb*t12*t13)*wght+Cmat2(iq,D2_RB_GBB) 318 Cmat2(iq,D2_GAA_GAA) = (4.2000000000000004d-5*t15*t8-2.646d- 319 1 7*gammaaa*t17/rhoa**6.666666666666667d+0)*wght+Cmat2(iq,D 320 2 2_GAA_GAA) 321 Cmat2(iq,D2_GAA_GAB) = Cmat2(iq,D2_GAA_GAB) 322 Cmat2(iq,D2_GAA_GBB) = Cmat2(iq,D2_GAA_GBB) 323 Cmat2(iq,D2_GAB_GAB) = Cmat2(iq,D2_GAB_GAB) 324 Cmat2(iq,D2_GAB_GBB) = Cmat2(iq,D2_GAB_GBB) 325 Cmat2(iq,D2_GBB_GBB) = (4.2000000000000004d-5*t12*t16-2.646d 326 1 -7*gammabb*t18/rhob**6.666666666666667d+0)*wght+Cmat2(iq, 327 2 D2_GBB_GBB) 328 elseif (rhoa.gt.tol_rho.and.rhob.le.tol_rho) then 329 t1 = 7.000000000000001d-3*gammaaa/rhoa**2.6666666666666666d+ 330 1 0+1.0d+0 331 t2 = 1/t1**8.0d-1 332 t3 = 1/rhoa**1.3333333333333333d+0 333 t4 = gammaaa**2 334 t5 = 1/t1**1.8d+0 335 t6 = 1/rhoa**5 336 t7 = 1/rhoa**2.3333333333333334d+0 337 t8 = 1/rhoa**4 338 t9 = 1/t1**2.7999999999999997d+0 339 fnc(iq) = fnc(iq)-3.75d-3*gammaaa*t2*t3*wght 340 Amat(iq,D1_RA) = 4.9999999999999994d-3*gammaaa*t2*t7*wght-5. 341 1 599999999999999d-5*t4*t5*t6*wght+Amat(iq,D1_RA) 342 Cmat(iq,D1_GAA) = 2.1000000000000002d-5*gammaaa*t5*t8*wght-3 343 1 .75d-3*t2*t3*wght+Cmat(iq,D1_GAA) 344 Amat2(iq,D2_RA_RA) = -1.8816d-6*gammaaa**3*t9*wght/rhoa**8.6 345 1 66666666666666d+0+3.5466666666666663d-4*t4*t5*wght/rhoa** 346 2 6-1.1666666666666665d-2*gammaaa*t2*wght/rhoa**3.333333333 347 3 3333337d+0+Amat2(iq,D2_RA_RA) 348 Cmat2(iq,D2_RA_GAA) = 7.056d-7*t4*t9*wght/rhoa**7.6666666666 349 1 66667d+0+4.9999999999999994d-3*t2*t7*wght-1.3999999999999 350 2 999d-4*gammaaa*t5*t6*wght+Cmat2(iq,D2_RA_GAA) 351 Cmat2(iq,D2_GAA_GAA) = -2.646d-7*gammaaa*t9*wght/rhoa**6.666 352 1 666666666667d+0+4.2000000000000004d-5*t5*t8*wght+Cmat2(iq 353 2 ,D2_GAA_GAA) 354 elseif (rhoa.le.tol_rho.and.rhob.gt.tol_rho) then 355 t1 = 7.000000000000001d-3*gammabb/rhob**2.6666666666666666d+ 356 1 0+1.0d+0 357 t2 = 1/t1**8.0d-1 358 t3 = 1/rhob**1.3333333333333333d+0 359 t4 = gammabb**2 360 t5 = 1/t1**1.8d+0 361 t6 = 1/rhob**5 362 t7 = 1/rhob**2.3333333333333334d+0 363 t8 = 1/rhob**4 364 t9 = 1/t1**2.7999999999999997d+0 365 fnc(iq) = fnc(iq)-3.75d-3*gammabb*t2*t3*wght 366 Amat(iq,D1_RB) = 4.9999999999999994d-3*gammabb*t2*t7*wght-5. 367 1 599999999999999d-5*t4*t5*t6*wght+Amat(iq,D1_RB) 368 Cmat(iq,D1_GBB) = 2.1000000000000002d-5*gammabb*t5*t8*wght-3 369 1 .75d-3*t2*t3*wght+Cmat(iq,D1_GBB) 370 Amat2(iq,D2_RB_RB) = -1.8816d-6*gammabb**3*t9*wght/rhob**8.6 371 1 66666666666666d+0+3.5466666666666663d-4*t4*t5*wght/rhob** 372 2 6-1.1666666666666665d-2*gammabb*t2*wght/rhob**3.333333333 373 3 3333337d+0+Amat2(iq,D2_RB_RB) 374 Cmat2(iq,D2_RB_GBB) = 7.056d-7*t4*t9*wght/rhob**7.6666666666 375 1 66667d+0+4.9999999999999994d-3*t2*t7*wght-1.3999999999999 376 2 999d-4*gammabb*t5*t6*wght+Cmat2(iq,D2_RB_GBB) 377 Cmat2(iq,D2_GBB_GBB) = -2.646d-7*gammabb*t9*wght/rhob**6.666 378 1 666666666667d+0+4.2000000000000004d-5*t5*t8*wght+Cmat2(iq 379 2 ,D2_GBB_GBB) 380 endif ! rhoa.gt.tol_rho.and.rhob.gt.tol_rho 381 endif ! ipol.eq.1 382 enddo ! iq 383 end 384C> 385C> \brief Evaluate the nwxcm_x_b86b functional [1] 386C> 387C> \f{eqnarray*}{ 388C> f &=& -{{0.00375\,\sigma_{\beta\beta}} 389C> \over{\rho_\beta^{{{4}\over{3}}}\,\left({{0.007\, 390C> \sigma_{\beta\beta}}\over{\rho_\beta^{{{8}\over{3}}}}} 391C> +1.0\right)^{0.8}}}-{{0.00375\,\sigma_{\alpha\alpha}} 392C> \over{\rho_\alpha^{{{4}\over{3}}}\,\left({{0.007\, 393C> \sigma_{\alpha\alpha}}\over{\rho_\alpha^{{{8}\over{3}}}}} 394C> +1.0\right)^{0.8}}}\\\\ 395C> g &=& 0\\\\ 396C> G &=& -{{0.00375\,\sigma_{ss}}\over{\rho_s^{{{4}\over{3}}} 397C> \,\left({{0.007\,\sigma_{ss}}\over{\rho_s^{{{8}\over{3}}}}} 398C> +1.0\right)^{0.8}}}\\\\ 399C> \f} 400C> 401C> Code generated with Maxima 5.34.0 [2] 402C> driven by autoxc [3]. 403C> 404C> ### References ### 405C> 406C> [1] AD Becke, J.Chem.Phys. 85, 7184 (1986) , DOI: 407C> <a href="https://doi.org/10.1063/1.451353 "> 408C> 10.1063/1.451353 </a> 409C> 410C> [2] Maxima, a computer algebra system, 411C> <a href="http://maxima.sourceforge.net/"> 412C> http://maxima.sourceforge.net/</a> 413C> 414C> [3] autoxc, revision 27097 2015-05-08 415C> 416 subroutine nwxcm_x_b86b_d3(param,tol_rho,ipol,nq,wght, 417 +rho,rgamma,fnc,Amat,Amat2,Amat3, 418 +Cmat,Cmat2,Cmat3) 419c $Id: $ 420#ifdef NWXC_QUAD_PREC 421 implicit real(kind=selected_real_kind(30))(a-h,o-z),integer(i-n) 422 integer, parameter :: rk=selected_real_kind(30) 423#else 424 implicit real(kind=selected_real_kind(15))(a-h,o-z),integer(i-n) 425 integer, parameter :: rk=selected_real_kind(15) 426#endif 427#include "nwxc_param.fh" 428 double precision param(*) !< [Input] Parameters of functional 429 double precision tol_rho !< [Input] The lower limit on the density 430 integer ipol !< [Input] The number of spin channels 431 integer nq !< [Input] The number of points 432 double precision wght !< [Input] The weight of the functional 433 double precision rho(nq,*) !< [Input] The density 434 double precision rgamma(nq,*) !< [Input] The norm of the density 435 !< gradients 436 double precision fnc(nq) !< [Output] The value of the functional 437c 438c Sampling Matrices for the XC Kernel 439c 440 double precision Amat(nq,*) !< [Output] The derivative wrt rho 441 double precision Cmat(nq,*) !< [Output] The derivative wrt rgamma 442c 443c Sampling Matrices for the XC Kernel 444c 445 double precision Amat2(nq,*) !< [Output] The 2nd derivative wrt rho 446 double precision Cmat2(nq,*) !< [Output] The 2nd derivative wrt rgamma 447 !< and possibly rho 448c 449c Sampling Matrices for the XC Kernel 450c 451 double precision Amat3(nq,*) !< [Output] The 3rd derivative wrt rho 452 double precision Cmat3(nq,*) !< [Output] The 3rd derivative wrt rgamma 453 !< and possibly rho 454 integer iq 455 double precision tmp 456 double precision rhoa,rhob 457 double precision gammaaa,gammaab,gammabb 458 double precision taua,taub 459 double precision nwxcm_heaviside 460 external nwxcm_heaviside 461CDIR$ NOVECTOR 462 do iq = 1, nq 463 if (ipol.eq.1) then 464 rhoa = 0.5d0*rho(iq,R_T) 465 gammaaa = 0.25d0*rgamma(iq,G_TT) 466 if (rhoa.gt.tol_rho) then 467 t1 = 7.000000000000001d-3*gammaaa/rhoa**2.6666666666666666d+ 468 1 0+1.0d+0 469 t2 = 1/t1**8.0d-1 470 t3 = 1/rhoa**1.3333333333333333d+0 471 t4 = gammaaa**2 472 t5 = 1/t1**1.8d+0 473 t6 = 1/rhoa**5 474 t7 = 1/rhoa**2.3333333333333334d+0 475 t8 = 1/rhoa**4 476 t9 = gammaaa**3 477 t10 = 1/t1**2.7999999999999997d+0 478 t11 = 1/rhoa**8.666666666666666d+0 479 t12 = 1/rhoa**6 480 t13 = 1/rhoa**3.3333333333333337d+0 481 t14 = 1/rhoa**7.666666666666667d+0 482 t15 = 1/rhoa**6.666666666666667d+0 483 t16 = 1/t1**3.8d+0 484 fnc(iq) = fnc(iq)-7.5d-3*gammaaa*t2*t3*wght 485 Amat(iq,D1_RA) = (4.9999999999999994d-3*gammaaa*t2*t7-5.5999 486 1 99999999999d-5*t4*t5*t6)*wght+Amat(iq,D1_RA) 487 Cmat(iq,D1_GAA) = (2.1000000000000002d-5*gammaaa*t5*t8-3.75d 488 1 -3*t2*t3)*wght+Cmat(iq,D1_GAA) 489 Cmat(iq,D1_GAB) = Cmat(iq,D1_GAB) 490 Amat2(iq,D2_RA_RA) = (-1.8816d-6*t10*t11*t9+3.54666666666666 491 1 63d-4*t12*t4*t5-1.1666666666666665d-2*gammaaa*t13*t2)*wgh 492 2 t+Amat2(iq,D2_RA_RA) 493 Amat2(iq,D2_RA_RB) = Amat2(iq,D2_RA_RB) 494 Cmat2(iq,D2_RA_GAA) = (4.9999999999999994d-3*t2*t7-1.3999999 495 1 999999999d-4*gammaaa*t5*t6+7.056d-7*t10*t14*t4)*wght+Cmat 496 2 2(iq,D2_RA_GAA) 497 Cmat2(iq,D2_RA_GAB) = Cmat2(iq,D2_RA_GAB) 498 Cmat2(iq,D2_RA_GBB) = Cmat2(iq,D2_RA_GBB) 499 Cmat2(iq,D2_GAA_GAA) = (4.2000000000000004d-5*t5*t8-2.646d-7 500 1 *gammaaa*t10*t15)*wght+Cmat2(iq,D2_GAA_GAA) 501 Cmat2(iq,D2_GAA_GAB) = Cmat2(iq,D2_GAA_GAB) 502 Cmat2(iq,D2_GAA_GBB) = Cmat2(iq,D2_GAA_GBB) 503 Cmat2(iq,D2_GAB_GAB) = Cmat2(iq,D2_GAB_GAB) 504 Amat3(iq,D3_RA_RA_RA) = (2.8224d-5*t10*t9/rhoa**9.6666666666 505 1 66666d+0-2.302222222222222d-3*t4*t5/rhoa**7+3.88888888888 506 2 8889d-2*gammaaa*t2/rhoa**4.333333333333333d+0-9.834495999 507 3 999997d-8*gammaaa**4*t16/rhoa**1.2333333333333334d+1)*wgh 508 4 t+Amat3(iq,D3_RA_RA_RA) 509 Amat3(iq,D3_RA_RA_RB) = Amat3(iq,D3_RA_RA_RB) 510 Cmat3(iq,D3_RA_RA_GAA) = (3.687936d-8*t16*t9/rhoa**1.1333333 511 1 333333334d+1+7.746666666666666d-4*gammaaa*t12*t5-1.011359 512 2 9999999999d-5*t10*t11*t4-1.1666666666666665d-2*t13*t2)*wg 513 3 ht+Cmat3(iq,D3_RA_RA_GAA) 514 Cmat3(iq,D3_RA_RA_GAB) = Cmat3(iq,D3_RA_RA_GAB) 515 Cmat3(iq,D3_RA_RA_GBB) = Cmat3(iq,D3_RA_RA_GBB) 516 Cmat3(iq,D3_RA_RB_GAA) = Cmat3(iq,D3_RA_RB_GAA) 517 Cmat3(iq,D3_RA_RB_GAB) = Cmat3(iq,D3_RA_RB_GAB) 518 Cmat3(iq,D3_RA_GAA_GAA) = (-1.6799999999999998d-4*t5*t6-1.38 519 1 2976d-8*t16*t4/rhoa**1.0333333333333333d+1+3.175199999999 520 2 9997d-6*gammaaa*t10*t14)*wght+Cmat3(iq,D3_RA_GAA_GAA) 521 Cmat3(iq,D3_RA_GAA_GAB) = Cmat3(iq,D3_RA_GAA_GAB) 522 Cmat3(iq,D3_RA_GAA_GBB) = Cmat3(iq,D3_RA_GAA_GBB) 523 Cmat3(iq,D3_RA_GAB_GAB) = Cmat3(iq,D3_RA_GAB_GAB) 524 Cmat3(iq,D3_RA_GAB_GBB) = Cmat3(iq,D3_RA_GAB_GBB) 525 Cmat3(iq,D3_RA_GBB_GBB) = Cmat3(iq,D3_RA_GBB_GBB) 526 Cmat3(iq,D3_GAA_GAA_GAA) = (5.186160000000001d-9*gammaaa*t16 527 1 /rhoa**9.333333333333333d+0-7.938000000000002d-7*t10*t15) 528 2 *wght+Cmat3(iq,D3_GAA_GAA_GAA) 529 Cmat3(iq,D3_GAA_GAA_GAB) = Cmat3(iq,D3_GAA_GAA_GAB) 530 Cmat3(iq,D3_GAA_GAA_GBB) = Cmat3(iq,D3_GAA_GAA_GBB) 531 Cmat3(iq,D3_GAA_GAB_GAB) = Cmat3(iq,D3_GAA_GAB_GAB) 532 Cmat3(iq,D3_GAA_GAB_GBB) = Cmat3(iq,D3_GAA_GAB_GBB) 533 Cmat3(iq,D3_GAB_GAB_GAB) = Cmat3(iq,D3_GAB_GAB_GAB) 534 endif ! rhoa.gt.tol_rho 535 else ! ipol.eq.1 536 rhoa = rho(iq,R_A) 537 rhob = rho(iq,R_B) 538 gammaaa = rgamma(iq,G_AA) 539 gammaab = rgamma(iq,G_AB) 540 gammabb = rgamma(iq,G_BB) 541 if (rhoa.gt.tol_rho.and.rhob.gt.tol_rho) then 542 t1 = 7.000000000000001d-3*gammaaa/rhoa**2.6666666666666666d+ 543 1 0+1.0d+0 544 t2 = 1/t1**8.0d-1 545 t3 = 1/rhoa**1.3333333333333333d+0 546 t4 = 7.000000000000001d-3*gammabb/rhob**2.6666666666666666d+ 547 1 0+1.0d+0 548 t5 = 1/t4**8.0d-1 549 t6 = 1/rhob**1.3333333333333333d+0 550 t7 = gammaaa**2 551 t8 = 1/t1**1.8d+0 552 t9 = 1/rhoa**5 553 t10 = 1/rhoa**2.3333333333333334d+0 554 t11 = gammabb**2 555 t12 = 1/t4**1.8d+0 556 t13 = 1/rhob**5 557 t14 = 1/rhob**2.3333333333333334d+0 558 t15 = 1/rhoa**4 559 t16 = 1/rhob**4 560 t17 = gammaaa**3 561 t18 = 1/t1**2.7999999999999997d+0 562 t19 = 1/rhoa**8.666666666666666d+0 563 t20 = 1/rhoa**6 564 t21 = 1/rhoa**3.3333333333333337d+0 565 t22 = gammabb**3 566 t23 = 1/t4**2.7999999999999997d+0 567 t24 = 1/rhob**8.666666666666666d+0 568 t25 = 1/rhob**6 569 t26 = 1/rhob**3.3333333333333337d+0 570 t27 = 1/rhoa**7.666666666666667d+0 571 t28 = 1/rhob**7.666666666666667d+0 572 t29 = 1/rhoa**6.666666666666667d+0 573 t30 = 1/rhob**6.666666666666667d+0 574 t31 = 1/t1**3.8d+0 575 t32 = 1/t4**3.8d+0 576 fnc(iq) = (-3.75d-3*gammabb*t5*t6-3.75d-3*gammaaa*t2*t3)*wgh 577 1 t+fnc(iq) 578 Amat(iq,D1_RA) = (4.9999999999999994d-3*gammaaa*t10*t2-5.599 579 1 999999999999d-5*t7*t8*t9)*wght+Amat(iq,D1_RA) 580 Amat(iq,D1_RB) = (4.9999999999999994d-3*gammabb*t14*t5-5.599 581 1 999999999999d-5*t11*t12*t13)*wght+Amat(iq,D1_RB) 582 Cmat(iq,D1_GAA) = (2.1000000000000002d-5*gammaaa*t15*t8-3.75 583 1 d-3*t2*t3)*wght+Cmat(iq,D1_GAA) 584 Cmat(iq,D1_GAB) = Cmat(iq,D1_GAB) 585 Cmat(iq,D1_GBB) = (2.1000000000000002d-5*gammabb*t12*t16-3.7 586 1 5d-3*t5*t6)*wght+Cmat(iq,D1_GBB) 587 Amat2(iq,D2_RA_RA) = (3.5466666666666663d-4*t20*t7*t8-1.1666 588 1 666666666665d-2*gammaaa*t2*t21-1.8816d-6*t17*t18*t19)*wgh 589 2 t+Amat2(iq,D2_RA_RA) 590 Amat2(iq,D2_RA_RB) = Amat2(iq,D2_RA_RB) 591 Amat2(iq,D2_RB_RB) = (-1.1666666666666665d-2*gammabb*t26*t5+ 592 1 3.5466666666666663d-4*t11*t12*t25-1.8816d-6*t22*t23*t24)* 593 2 wght+Amat2(iq,D2_RB_RB) 594 Cmat2(iq,D2_RA_GAA) = (-1.3999999999999999d-4*gammaaa*t8*t9+ 595 1 7.056d-7*t18*t27*t7+4.9999999999999994d-3*t10*t2)*wght+Cm 596 2 at2(iq,D2_RA_GAA) 597 Cmat2(iq,D2_RA_GAB) = Cmat2(iq,D2_RA_GAB) 598 Cmat2(iq,D2_RA_GBB) = Cmat2(iq,D2_RA_GBB) 599 Cmat2(iq,D2_RB_GAA) = Cmat2(iq,D2_RB_GAA) 600 Cmat2(iq,D2_RB_GAB) = Cmat2(iq,D2_RB_GAB) 601 Cmat2(iq,D2_RB_GBB) = (4.9999999999999994d-3*t14*t5+7.056d-7 602 1 *t11*t23*t28-1.3999999999999999d-4*gammabb*t12*t13)*wght+ 603 2 Cmat2(iq,D2_RB_GBB) 604 Cmat2(iq,D2_GAA_GAA) = (4.2000000000000004d-5*t15*t8-2.646d- 605 1 7*gammaaa*t18*t29)*wght+Cmat2(iq,D2_GAA_GAA) 606 Cmat2(iq,D2_GAA_GAB) = Cmat2(iq,D2_GAA_GAB) 607 Cmat2(iq,D2_GAA_GBB) = Cmat2(iq,D2_GAA_GBB) 608 Cmat2(iq,D2_GAB_GAB) = Cmat2(iq,D2_GAB_GAB) 609 Cmat2(iq,D2_GAB_GBB) = Cmat2(iq,D2_GAB_GBB) 610 Cmat2(iq,D2_GBB_GBB) = (4.2000000000000004d-5*t12*t16-2.646d 611 1 -7*gammabb*t23*t30)*wght+Cmat2(iq,D2_GBB_GBB) 612 Amat3(iq,D3_RA_RA_RA) = (-2.302222222222222d-3*t7*t8/rhoa**7 613 1 -9.834495999999997d-8*gammaaa**4*t31/rhoa**1.233333333333 614 2 3334d+1+3.888888888888889d-2*gammaaa*t2/rhoa**4.333333333 615 3 333333d+0+2.8224d-5*t17*t18/rhoa**9.666666666666666d+0)*w 616 4 ght+Amat3(iq,D3_RA_RA_RA) 617 Amat3(iq,D3_RA_RA_RB) = Amat3(iq,D3_RA_RA_RB) 618 Amat3(iq,D3_RA_RB_RB) = Amat3(iq,D3_RA_RB_RB) 619 Amat3(iq,D3_RB_RB_RB) = (3.888888888888889d-2*gammabb*t5/rho 620 1 b**4.333333333333333d+0-9.834495999999997d-8*gammabb**4*t 621 2 32/rhob**1.2333333333333334d+1+2.8224d-5*t22*t23/rhob**9. 622 3 666666666666666d+0-2.302222222222222d-3*t11*t12/rhob**7)* 623 4 wght+Amat3(iq,D3_RB_RB_RB) 624 Cmat3(iq,D3_RA_RA_GAA) = (7.746666666666666d-4*gammaaa*t20*t 625 1 8-1.0113599999999999d-5*t18*t19*t7+3.687936d-8*t17*t31/rh 626 2 oa**1.1333333333333334d+1-1.1666666666666665d-2*t2*t21)*w 627 3 ght+Cmat3(iq,D3_RA_RA_GAA) 628 Cmat3(iq,D3_RA_RA_GAB) = Cmat3(iq,D3_RA_RA_GAB) 629 Cmat3(iq,D3_RA_RA_GBB) = Cmat3(iq,D3_RA_RA_GBB) 630 Cmat3(iq,D3_RA_RB_GAA) = Cmat3(iq,D3_RA_RB_GAA) 631 Cmat3(iq,D3_RA_RB_GAB) = Cmat3(iq,D3_RA_RB_GAB) 632 Cmat3(iq,D3_RA_RB_GBB) = Cmat3(iq,D3_RA_RB_GBB) 633 Cmat3(iq,D3_RB_RB_GAA) = Cmat3(iq,D3_RB_RB_GAA) 634 Cmat3(iq,D3_RB_RB_GAB) = Cmat3(iq,D3_RB_RB_GAB) 635 Cmat3(iq,D3_RB_RB_GBB) = (-1.1666666666666665d-2*t26*t5+3.68 636 1 7936d-8*t22*t32/rhob**1.1333333333333334d+1+7.74666666666 637 2 6666d-4*gammabb*t12*t25-1.0113599999999999d-5*t11*t23*t24 638 3 )*wght+Cmat3(iq,D3_RB_RB_GBB) 639 Cmat3(iq,D3_RA_GAA_GAA) = (-1.6799999999999998d-4*t8*t9-1.38 640 1 2976d-8*t31*t7/rhoa**1.0333333333333333d+1+3.175199999999 641 2 9997d-6*gammaaa*t18*t27)*wght+Cmat3(iq,D3_RA_GAA_GAA) 642 Cmat3(iq,D3_RA_GAA_GAB) = Cmat3(iq,D3_RA_GAA_GAB) 643 Cmat3(iq,D3_RA_GAA_GBB) = Cmat3(iq,D3_RA_GAA_GBB) 644 Cmat3(iq,D3_RA_GAB_GAB) = Cmat3(iq,D3_RA_GAB_GAB) 645 Cmat3(iq,D3_RA_GAB_GBB) = Cmat3(iq,D3_RA_GAB_GBB) 646 Cmat3(iq,D3_RA_GBB_GBB) = Cmat3(iq,D3_RA_GBB_GBB) 647 Cmat3(iq,D3_RB_GAA_GAA) = Cmat3(iq,D3_RB_GAA_GAA) 648 Cmat3(iq,D3_RB_GAA_GAB) = Cmat3(iq,D3_RB_GAA_GAB) 649 Cmat3(iq,D3_RB_GAA_GBB) = Cmat3(iq,D3_RB_GAA_GBB) 650 Cmat3(iq,D3_RB_GAB_GAB) = Cmat3(iq,D3_RB_GAB_GAB) 651 Cmat3(iq,D3_RB_GAB_GBB) = Cmat3(iq,D3_RB_GAB_GBB) 652 Cmat3(iq,D3_RB_GBB_GBB) = (-1.382976d-8*t11*t32/rhob**1.0333 653 1 333333333333d+1+3.1751999999999997d-6*gammabb*t23*t28-1.6 654 2 799999999999998d-4*t12*t13)*wght+Cmat3(iq,D3_RB_GBB_GBB) 655 Cmat3(iq,D3_GAA_GAA_GAA) = (5.186160000000001d-9*gammaaa*t31 656 1 /rhoa**9.333333333333333d+0-7.938000000000002d-7*t18*t29) 657 2 *wght+Cmat3(iq,D3_GAA_GAA_GAA) 658 Cmat3(iq,D3_GAA_GAA_GAB) = Cmat3(iq,D3_GAA_GAA_GAB) 659 Cmat3(iq,D3_GAA_GAA_GBB) = Cmat3(iq,D3_GAA_GAA_GBB) 660 Cmat3(iq,D3_GAA_GAB_GAB) = Cmat3(iq,D3_GAA_GAB_GAB) 661 Cmat3(iq,D3_GAA_GAB_GBB) = Cmat3(iq,D3_GAA_GAB_GBB) 662 Cmat3(iq,D3_GAA_GBB_GBB) = Cmat3(iq,D3_GAA_GBB_GBB) 663 Cmat3(iq,D3_GAB_GAB_GAB) = Cmat3(iq,D3_GAB_GAB_GAB) 664 Cmat3(iq,D3_GAB_GAB_GBB) = Cmat3(iq,D3_GAB_GAB_GBB) 665 Cmat3(iq,D3_GAB_GBB_GBB) = Cmat3(iq,D3_GAB_GBB_GBB) 666 Cmat3(iq,D3_GBB_GBB_GBB) = (5.186160000000001d-9*gammabb*t32 667 1 /rhob**9.333333333333333d+0-7.938000000000002d-7*t23*t30) 668 2 *wght+Cmat3(iq,D3_GBB_GBB_GBB) 669 elseif (rhoa.gt.tol_rho.and.rhob.le.tol_rho) then 670 t1 = 7.000000000000001d-3*gammaaa/rhoa**2.6666666666666666d+ 671 1 0+1.0d+0 672 t2 = 1/t1**8.0d-1 673 t3 = 1/rhoa**1.3333333333333333d+0 674 t4 = gammaaa**2 675 t5 = 1/t1**1.8d+0 676 t6 = 1/rhoa**5 677 t7 = 1/rhoa**2.3333333333333334d+0 678 t8 = 1/rhoa**4 679 t9 = gammaaa**3 680 t10 = 1/t1**2.7999999999999997d+0 681 t11 = 1/rhoa**8.666666666666666d+0 682 t12 = 1/rhoa**6 683 t13 = 1/rhoa**3.3333333333333337d+0 684 t14 = 1/rhoa**7.666666666666667d+0 685 t15 = 1/rhoa**6.666666666666667d+0 686 t16 = 1/t1**3.8d+0 687 fnc(iq) = fnc(iq)-3.75d-3*gammaaa*t2*t3*wght 688 Amat(iq,D1_RA) = 4.9999999999999994d-3*gammaaa*t2*t7*wght-5. 689 1 599999999999999d-5*t4*t5*t6*wght+Amat(iq,D1_RA) 690 Cmat(iq,D1_GAA) = 2.1000000000000002d-5*gammaaa*t5*t8*wght-3 691 1 .75d-3*t2*t3*wght+Cmat(iq,D1_GAA) 692 Amat2(iq,D2_RA_RA) = -1.8816d-6*t10*t11*t9*wght+3.5466666666 693 1 666663d-4*t12*t4*t5*wght-1.1666666666666665d-2*gammaaa*t1 694 2 3*t2*wght+Amat2(iq,D2_RA_RA) 695 Cmat2(iq,D2_RA_GAA) = 4.9999999999999994d-3*t2*t7*wght-1.399 696 1 9999999999999d-4*gammaaa*t5*t6*wght+7.056d-7*t10*t14*t4*w 697 2 ght+Cmat2(iq,D2_RA_GAA) 698 Cmat2(iq,D2_GAA_GAA) = 4.2000000000000004d-5*t5*t8*wght-2.64 699 1 6d-7*gammaaa*t10*t15*wght+Cmat2(iq,D2_GAA_GAA) 700 Amat3(iq,D3_RA_RA_RA) = 2.8224d-5*t10*t9*wght/rhoa**9.666666 701 1 666666666d+0-2.302222222222222d-3*t4*t5*wght/rhoa**7+3.88 702 2 8888888888889d-2*gammaaa*t2*wght/rhoa**4.333333333333333d 703 3 +0-9.834495999999997d-8*gammaaa**4*t16*wght/rhoa**1.23333 704 4 33333333334d+1+Amat3(iq,D3_RA_RA_RA) 705 Cmat3(iq,D3_RA_RA_GAA) = 3.687936d-8*t16*t9*wght/rhoa**1.133 706 1 3333333333334d+1+7.746666666666666d-4*gammaaa*t12*t5*wght 707 2 -1.0113599999999999d-5*t10*t11*t4*wght-1.1666666666666665 708 3 d-2*t13*t2*wght+Cmat3(iq,D3_RA_RA_GAA) 709 Cmat3(iq,D3_RA_GAA_GAA) = -1.6799999999999998d-4*t5*t6*wght- 710 1 1.382976d-8*t16*t4*wght/rhoa**1.0333333333333333d+1+3.175 711 2 1999999999997d-6*gammaaa*t10*t14*wght+Cmat3(iq,D3_RA_GAA_ 712 3 GAA) 713 Cmat3(iq,D3_GAA_GAA_GAA) = 5.186160000000001d-9*gammaaa*t16* 714 1 wght/rhoa**9.333333333333333d+0-7.938000000000002d-7*t10* 715 2 t15*wght+Cmat3(iq,D3_GAA_GAA_GAA) 716 elseif (rhoa.le.tol_rho.and.rhob.gt.tol_rho) then 717 t1 = 7.000000000000001d-3*gammabb/rhob**2.6666666666666666d+ 718 1 0+1.0d+0 719 t2 = 1/t1**8.0d-1 720 t3 = 1/rhob**1.3333333333333333d+0 721 t4 = gammabb**2 722 t5 = 1/t1**1.8d+0 723 t6 = 1/rhob**5 724 t7 = 1/rhob**2.3333333333333334d+0 725 t8 = 1/rhob**4 726 t9 = gammabb**3 727 t10 = 1/t1**2.7999999999999997d+0 728 t11 = 1/rhob**8.666666666666666d+0 729 t12 = 1/rhob**6 730 t13 = 1/rhob**3.3333333333333337d+0 731 t14 = 1/rhob**7.666666666666667d+0 732 t15 = 1/rhob**6.666666666666667d+0 733 t16 = 1/t1**3.8d+0 734 fnc(iq) = fnc(iq)-3.75d-3*gammabb*t2*t3*wght 735 Amat(iq,D1_RB) = 4.9999999999999994d-3*gammabb*t2*t7*wght-5. 736 1 599999999999999d-5*t4*t5*t6*wght+Amat(iq,D1_RB) 737 Cmat(iq,D1_GBB) = 2.1000000000000002d-5*gammabb*t5*t8*wght-3 738 1 .75d-3*t2*t3*wght+Cmat(iq,D1_GBB) 739 Amat2(iq,D2_RB_RB) = -1.8816d-6*t10*t11*t9*wght+3.5466666666 740 1 666663d-4*t12*t4*t5*wght-1.1666666666666665d-2*gammabb*t1 741 2 3*t2*wght+Amat2(iq,D2_RB_RB) 742 Cmat2(iq,D2_RB_GBB) = 4.9999999999999994d-3*t2*t7*wght-1.399 743 1 9999999999999d-4*gammabb*t5*t6*wght+7.056d-7*t10*t14*t4*w 744 2 ght+Cmat2(iq,D2_RB_GBB) 745 Cmat2(iq,D2_GBB_GBB) = 4.2000000000000004d-5*t5*t8*wght-2.64 746 1 6d-7*gammabb*t10*t15*wght+Cmat2(iq,D2_GBB_GBB) 747 Amat3(iq,D3_RB_RB_RB) = 2.8224d-5*t10*t9*wght/rhob**9.666666 748 1 666666666d+0-2.302222222222222d-3*t4*t5*wght/rhob**7+3.88 749 2 8888888888889d-2*gammabb*t2*wght/rhob**4.333333333333333d 750 3 +0-9.834495999999997d-8*gammabb**4*t16*wght/rhob**1.23333 751 4 33333333334d+1+Amat3(iq,D3_RB_RB_RB) 752 Cmat3(iq,D3_RB_RB_GBB) = 3.687936d-8*t16*t9*wght/rhob**1.133 753 1 3333333333334d+1+7.746666666666666d-4*gammabb*t12*t5*wght 754 2 -1.0113599999999999d-5*t10*t11*t4*wght-1.1666666666666665 755 3 d-2*t13*t2*wght+Cmat3(iq,D3_RB_RB_GBB) 756 Cmat3(iq,D3_RB_GBB_GBB) = -1.6799999999999998d-4*t5*t6*wght- 757 1 1.382976d-8*t16*t4*wght/rhob**1.0333333333333333d+1+3.175 758 2 1999999999997d-6*gammabb*t10*t14*wght+Cmat3(iq,D3_RB_GBB_ 759 3 GBB) 760 Cmat3(iq,D3_GBB_GBB_GBB) = 5.186160000000001d-9*gammabb*t16* 761 1 wght/rhob**9.333333333333333d+0-7.938000000000002d-7*t10* 762 2 t15*wght+Cmat3(iq,D3_GBB_GBB_GBB) 763 endif ! rhoa.gt.tol_rho.and.rhob.gt.tol_rho 764 endif ! ipol.eq.1 765 enddo ! iq 766 end 767C> @} 768