1C> \ingroup nwxc 2C> @{ 3C> 4C> \file nwxcm_x_gill.F 5C> The nwxcm_x_gill functional 6C> 7C> @} 8C> 9C> \ingroup nwxc_priv 10C> @{ 11C> 12C> \brief Evaluate the nwxcm_x_gill functional [1] 13C> 14C> \f{eqnarray*}{ 15C> f &=& -{{0.0072992700729927\,\sigma_{\beta\beta}^{{{3} 16C> \over{4}}}}\over{\rho_\beta^{{{2}\over{3}}}}} 17C> -{{0.0072992700729927\,\sigma_{\alpha\alpha}^{{{3} 18C> \over{4}}}}\over{\rho_\alpha^{{{2}\over{3}}}}} 19C> -0.9305257363491002\,\rho_\beta^{{{4}\over{3}}} 20C> -0.9305257363491002\,\rho_\alpha^{{{4}\over{3}}}\\\\ 21C> g &=& 0\\\\ 22C> G &=& -{{0.0072992700729927\,\sigma_{ss}^{{{3}\over{4}}}} 23C> \over{\rho_s^{{{2}\over{3}}}}}-0.9305257363491002 24C> \,\rho_s^{{{4}\over{3}}}\\\\ 25C> \f} 26C> 27C> Code generated with Maxima 5.34.0 [2] 28C> driven by autoxc [3]. 29C> 30C> ### References ### 31C> 32C> [1] PMW Gill, Mol.Phys. 89, 433 (1996) , DOI: 33C> <a href="https://doi.org/10.1080/002689796173813 "> 34C> 10.1080/002689796173813 </a> 35C> 36C> [2] Maxima, a computer algebra system, 37C> <a href="http://maxima.sourceforge.net/"> 38C> http://maxima.sourceforge.net/</a> 39C> 40C> [3] autoxc, revision 27097 2015-05-08 41C> 42 subroutine nwxcm_x_gill(param,tol_rho,ipol,nq,wght, 43 +rho,rgamma,fnc,Amat,Cmat) 44c $Id: $ 45#ifdef NWXC_QUAD_PREC 46 implicit real(kind=selected_real_kind(30))(a-h,o-z),integer(i-n) 47 integer, parameter :: rk=selected_real_kind(30) 48#else 49 implicit real(kind=selected_real_kind(15))(a-h,o-z),integer(i-n) 50 integer, parameter :: rk=selected_real_kind(15) 51#endif 52#include "nwxc_param.fh" 53 double precision param(*) !< [Input] Parameters of functional 54 double precision tol_rho !< [Input] The lower limit on the density 55 integer ipol !< [Input] The number of spin channels 56 integer nq !< [Input] The number of points 57 double precision wght !< [Input] The weight of the functional 58 double precision rho(nq,*) !< [Input] The density 59 double precision rgamma(nq,*) !< [Input] The norm of the density 60 !< gradients 61 double precision fnc(nq) !< [Output] The value of the functional 62c 63c Sampling Matrices for the XC Kernel 64c 65 double precision Amat(nq,*) !< [Output] The derivative wrt rho 66 double precision Cmat(nq,*) !< [Output] The derivative wrt rgamma 67 integer iq 68 double precision tmp 69 double precision rhoa,rhob 70 double precision gammaaa,gammaab,gammabb 71 double precision taua,taub 72 double precision nwxcm_heaviside 73 external nwxcm_heaviside 74CDIR$ NOVECTOR 75 do iq = 1, nq 76 if (ipol.eq.1) then 77 rhoa = 0.5d0*rho(iq,R_T) 78 gammaaa = 0.25d0*rgamma(iq,G_TT) 79 if (rhoa.gt.tol_rho) then 80 t1 = gammaaa**7.5d-1 81 t2 = 1/rhoa**6.666666666666666d-1 82 fnc(iq) = (-1.45985401459854d-2*t1*t2-1.8610514726982003d+0* 83 1 rhoa**1.3333333333333333d+0)*wght+fnc(iq) 84 Amat(iq,D1_RA) = (4.8661800486617995d-3*t1/rhoa**1.666666666 85 1 6666669d+0-1.2407009817988002d+0*rhoa**3.333333333333333d 86 2 -1)*wght+Amat(iq,D1_RA) 87 Cmat(iq,D1_GAA) = Cmat(iq,D1_GAA)-5.474452554744524d-3*t2*wg 88 1 ht/gammaaa**2.5d-1 89 Cmat(iq,D1_GAB) = Cmat(iq,D1_GAB) 90 endif ! rhoa.gt.tol_rho 91 else ! ipol.eq.1 92 rhoa = rho(iq,R_A) 93 rhob = rho(iq,R_B) 94 gammaaa = rgamma(iq,G_AA) 95 gammaab = rgamma(iq,G_AB) 96 gammabb = rgamma(iq,G_BB) 97 if (rhoa.gt.tol_rho.and.rhob.gt.tol_rho) then 98 t1 = gammaaa**7.5d-1 99 t2 = 1/rhoa**6.666666666666666d-1 100 t3 = gammabb**7.5d-1 101 t4 = 1/rhob**6.666666666666666d-1 102 fnc(iq) = (-7.2992700729927d-3*t3*t4-7.2992700729927d-3*t1*t 103 1 2-9.305257363491002d-1*rhob**1.3333333333333333d+0-9.3052 104 2 57363491002d-1*rhoa**1.3333333333333333d+0)*wght+fnc(iq) 105 Amat(iq,D1_RA) = (4.8661800486617995d-3*t1/rhoa**1.666666666 106 1 6666669d+0-1.2407009817988002d+0*rhoa**3.333333333333333d 107 2 -1)*wght+Amat(iq,D1_RA) 108 Amat(iq,D1_RB) = (4.8661800486617995d-3*t3/rhob**1.666666666 109 1 6666669d+0-1.2407009817988002d+0*rhob**3.333333333333333d 110 2 -1)*wght+Amat(iq,D1_RB) 111 Cmat(iq,D1_GAA) = Cmat(iq,D1_GAA)-5.474452554744524d-3*t2*wg 112 1 ht/gammaaa**2.5d-1 113 Cmat(iq,D1_GAB) = Cmat(iq,D1_GAB) 114 Cmat(iq,D1_GBB) = Cmat(iq,D1_GBB)-5.474452554744524d-3*t4*wg 115 1 ht/gammabb**2.5d-1 116 elseif (rhoa.gt.tol_rho.and.rhob.le.tol_rho) then 117 t1 = gammaaa**7.5d-1 118 t2 = 1/rhoa**6.666666666666666d-1 119 fnc(iq) = (-7.2992700729927d-3*t1*t2-9.305257363491002d-1*rh 120 1 oa**1.3333333333333333d+0)*wght+fnc(iq) 121 Amat(iq,D1_RA) = (4.8661800486617995d-3*t1/rhoa**1.666666666 122 1 6666669d+0-1.2407009817988002d+0*rhoa**3.333333333333333d 123 2 -1)*wght+Amat(iq,D1_RA) 124 Cmat(iq,D1_GAA) = Cmat(iq,D1_GAA)-5.474452554744524d-3*t2*wg 125 1 ht/gammaaa**2.5d-1 126 elseif (rhoa.le.tol_rho.and.rhob.gt.tol_rho) then 127 t1 = gammabb**7.5d-1 128 t2 = 1/rhob**6.666666666666666d-1 129 fnc(iq) = (-7.2992700729927d-3*t1*t2-9.305257363491002d-1*rh 130 1 ob**1.3333333333333333d+0)*wght+fnc(iq) 131 Amat(iq,D1_RB) = (4.8661800486617995d-3*t1/rhob**1.666666666 132 1 6666669d+0-1.2407009817988002d+0*rhob**3.333333333333333d 133 2 -1)*wght+Amat(iq,D1_RB) 134 Cmat(iq,D1_GBB) = Cmat(iq,D1_GBB)-5.474452554744524d-3*t2*wg 135 1 ht/gammabb**2.5d-1 136 endif ! rhoa.gt.tol_rho.and.rhob.gt.tol_rho 137 endif ! ipol.eq.1 138 enddo ! iq 139 end 140C> 141C> \brief Evaluate the nwxcm_x_gill functional [1] 142C> 143C> \f{eqnarray*}{ 144C> f &=& -{{0.0072992700729927\,\sigma_{\beta\beta}^{{{3} 145C> \over{4}}}}\over{\rho_\beta^{{{2}\over{3}}}}} 146C> -{{0.0072992700729927\,\sigma_{\alpha\alpha}^{{{3} 147C> \over{4}}}}\over{\rho_\alpha^{{{2}\over{3}}}}} 148C> -0.9305257363491002\,\rho_\beta^{{{4}\over{3}}} 149C> -0.9305257363491002\,\rho_\alpha^{{{4}\over{3}}}\\\\ 150C> g &=& 0\\\\ 151C> G &=& -{{0.0072992700729927\,\sigma_{ss}^{{{3}\over{4}}}} 152C> \over{\rho_s^{{{2}\over{3}}}}}-0.9305257363491002 153C> \,\rho_s^{{{4}\over{3}}}\\\\ 154C> \f} 155C> 156C> Code generated with Maxima 5.34.0 [2] 157C> driven by autoxc [3]. 158C> 159C> ### References ### 160C> 161C> [1] PMW Gill, Mol.Phys. 89, 433 (1996) , DOI: 162C> <a href="https://doi.org/10.1080/002689796173813 "> 163C> 10.1080/002689796173813 </a> 164C> 165C> [2] Maxima, a computer algebra system, 166C> <a href="http://maxima.sourceforge.net/"> 167C> http://maxima.sourceforge.net/</a> 168C> 169C> [3] autoxc, revision 27097 2015-05-08 170C> 171 subroutine nwxcm_x_gill_d2(param,tol_rho,ipol,nq,wght, 172 +rho,rgamma,fnc,Amat,Amat2,Cmat,Cmat2) 173c $Id: $ 174#ifdef NWXC_QUAD_PREC 175 implicit real(kind=selected_real_kind(30))(a-h,o-z),integer(i-n) 176 integer, parameter :: rk=selected_real_kind(30) 177#else 178 implicit real(kind=selected_real_kind(15))(a-h,o-z),integer(i-n) 179 integer, parameter :: rk=selected_real_kind(15) 180#endif 181#include "nwxc_param.fh" 182 double precision param(*) !< [Input] Parameters of functional 183 double precision tol_rho !< [Input] The lower limit on the density 184 integer ipol !< [Input] The number of spin channels 185 integer nq !< [Input] The number of points 186 double precision wght !< [Input] The weight of the functional 187 double precision rho(nq,*) !< [Input] The density 188 double precision rgamma(nq,*) !< [Input] The norm of the density 189 !< gradients 190 double precision fnc(nq) !< [Output] The value of the functional 191c 192c Sampling Matrices for the XC Kernel 193c 194 double precision Amat(nq,*) !< [Output] The derivative wrt rho 195 double precision Cmat(nq,*) !< [Output] The derivative wrt rgamma 196c 197c Sampling Matrices for the XC Kernel 198c 199 double precision Amat2(nq,*) !< [Output] The 2nd derivative wrt rho 200 double precision Cmat2(nq,*) !< [Output] The 2nd derivative wrt rgamma 201 !< and possibly rho 202 integer iq 203 double precision tmp 204 double precision rhoa,rhob 205 double precision gammaaa,gammaab,gammabb 206 double precision taua,taub 207 double precision nwxcm_heaviside 208 external nwxcm_heaviside 209CDIR$ NOVECTOR 210 do iq = 1, nq 211 if (ipol.eq.1) then 212 rhoa = 0.5d0*rho(iq,R_T) 213 gammaaa = 0.25d0*rgamma(iq,G_TT) 214 if (rhoa.gt.tol_rho) then 215 t1 = gammaaa**7.5d-1 216 t2 = 1/rhoa**6.666666666666666d-1 217 t3 = 1/rhoa**1.6666666666666669d+0 218 t4 = 1/gammaaa**2.5d-1 219 fnc(iq) = (-1.45985401459854d-2*t1*t2-1.8610514726982003d+0* 220 1 rhoa**1.3333333333333333d+0)*wght+fnc(iq) 221 Amat(iq,D1_RA) = (4.8661800486617995d-3*t1*t3-1.240700981798 222 1 8002d+0*rhoa**3.333333333333333d-1)*wght+Amat(iq,D1_RA) 223 Cmat(iq,D1_GAA) = Cmat(iq,D1_GAA)-5.474452554744524d-3*t2*t4 224 1 *wght 225 Cmat(iq,D1_GAB) = Cmat(iq,D1_GAB) 226 Amat2(iq,D2_RA_RA) = (-4.135669939329334d-1*t2-8.11030008110 227 1 3d-3*t1/rhoa**2.6666666666666666d+0)*wght+Amat2(iq,D2_RA_ 228 2 RA) 229 Amat2(iq,D2_RA_RB) = Amat2(iq,D2_RA_RB) 230 Cmat2(iq,D2_RA_GAA) = 3.64963503649635d-3*t3*t4*wght+Cmat2(i 231 1 q,D2_RA_GAA) 232 Cmat2(iq,D2_RA_GAB) = Cmat2(iq,D2_RA_GAB) 233 Cmat2(iq,D2_RA_GBB) = Cmat2(iq,D2_RA_GBB) 234 Cmat2(iq,D2_GAA_GAA) = 1.368613138686131d-3*t2*wght/gammaaa* 235 1 *1.25d+0+Cmat2(iq,D2_GAA_GAA) 236 Cmat2(iq,D2_GAA_GAB) = Cmat2(iq,D2_GAA_GAB) 237 Cmat2(iq,D2_GAA_GBB) = Cmat2(iq,D2_GAA_GBB) 238 Cmat2(iq,D2_GAB_GAB) = Cmat2(iq,D2_GAB_GAB) 239 endif ! rhoa.gt.tol_rho 240 else ! ipol.eq.1 241 rhoa = rho(iq,R_A) 242 rhob = rho(iq,R_B) 243 gammaaa = rgamma(iq,G_AA) 244 gammaab = rgamma(iq,G_AB) 245 gammabb = rgamma(iq,G_BB) 246 if (rhoa.gt.tol_rho.and.rhob.gt.tol_rho) then 247 t1 = gammaaa**7.5d-1 248 t2 = 1/rhoa**6.666666666666666d-1 249 t3 = gammabb**7.5d-1 250 t4 = 1/rhob**6.666666666666666d-1 251 t5 = 1/rhoa**1.6666666666666669d+0 252 t6 = 1/rhob**1.6666666666666669d+0 253 t7 = 1/gammaaa**2.5d-1 254 t8 = 1/gammabb**2.5d-1 255 fnc(iq) = (-7.2992700729927d-3*t3*t4-7.2992700729927d-3*t1*t 256 1 2-9.305257363491002d-1*rhob**1.3333333333333333d+0-9.3052 257 2 57363491002d-1*rhoa**1.3333333333333333d+0)*wght+fnc(iq) 258 Amat(iq,D1_RA) = (4.8661800486617995d-3*t1*t5-1.240700981798 259 1 8002d+0*rhoa**3.333333333333333d-1)*wght+Amat(iq,D1_RA) 260 Amat(iq,D1_RB) = (4.8661800486617995d-3*t3*t6-1.240700981798 261 1 8002d+0*rhob**3.333333333333333d-1)*wght+Amat(iq,D1_RB) 262 Cmat(iq,D1_GAA) = Cmat(iq,D1_GAA)-5.474452554744524d-3*t2*t7 263 1 *wght 264 Cmat(iq,D1_GAB) = Cmat(iq,D1_GAB) 265 Cmat(iq,D1_GBB) = Cmat(iq,D1_GBB)-5.474452554744524d-3*t4*t8 266 1 *wght 267 Amat2(iq,D2_RA_RA) = (-4.135669939329334d-1*t2-8.11030008110 268 1 3d-3*t1/rhoa**2.6666666666666666d+0)*wght+Amat2(iq,D2_RA_ 269 2 RA) 270 Amat2(iq,D2_RA_RB) = Amat2(iq,D2_RA_RB) 271 Amat2(iq,D2_RB_RB) = (-4.135669939329334d-1*t4-8.11030008110 272 1 3d-3*t3/rhob**2.6666666666666666d+0)*wght+Amat2(iq,D2_RB_ 273 2 RB) 274 Cmat2(iq,D2_RA_GAA) = 3.64963503649635d-3*t5*t7*wght+Cmat2(i 275 1 q,D2_RA_GAA) 276 Cmat2(iq,D2_RA_GAB) = Cmat2(iq,D2_RA_GAB) 277 Cmat2(iq,D2_RA_GBB) = Cmat2(iq,D2_RA_GBB) 278 Cmat2(iq,D2_RB_GAA) = Cmat2(iq,D2_RB_GAA) 279 Cmat2(iq,D2_RB_GAB) = Cmat2(iq,D2_RB_GAB) 280 Cmat2(iq,D2_RB_GBB) = 3.64963503649635d-3*t6*t8*wght+Cmat2(i 281 1 q,D2_RB_GBB) 282 Cmat2(iq,D2_GAA_GAA) = 1.368613138686131d-3*t2*wght/gammaaa* 283 1 *1.25d+0+Cmat2(iq,D2_GAA_GAA) 284 Cmat2(iq,D2_GAA_GAB) = Cmat2(iq,D2_GAA_GAB) 285 Cmat2(iq,D2_GAA_GBB) = Cmat2(iq,D2_GAA_GBB) 286 Cmat2(iq,D2_GAB_GAB) = Cmat2(iq,D2_GAB_GAB) 287 Cmat2(iq,D2_GAB_GBB) = Cmat2(iq,D2_GAB_GBB) 288 Cmat2(iq,D2_GBB_GBB) = 1.368613138686131d-3*t4*wght/gammabb* 289 1 *1.25d+0+Cmat2(iq,D2_GBB_GBB) 290 elseif (rhoa.gt.tol_rho.and.rhob.le.tol_rho) then 291 t1 = gammaaa**7.5d-1 292 t2 = 1/rhoa**6.666666666666666d-1 293 t3 = 1/rhoa**1.6666666666666669d+0 294 t4 = 1/gammaaa**2.5d-1 295 fnc(iq) = (-7.2992700729927d-3*t1*t2-9.305257363491002d-1*rh 296 1 oa**1.3333333333333333d+0)*wght+fnc(iq) 297 Amat(iq,D1_RA) = (4.8661800486617995d-3*t1*t3-1.240700981798 298 1 8002d+0*rhoa**3.333333333333333d-1)*wght+Amat(iq,D1_RA) 299 Cmat(iq,D1_GAA) = Cmat(iq,D1_GAA)-5.474452554744524d-3*t2*t4 300 1 *wght 301 Amat2(iq,D2_RA_RA) = (-4.135669939329334d-1*t2-8.11030008110 302 1 3d-3*t1/rhoa**2.6666666666666666d+0)*wght+Amat2(iq,D2_RA_ 303 2 RA) 304 Cmat2(iq,D2_RA_GAA) = 3.64963503649635d-3*t3*t4*wght+Cmat2(i 305 1 q,D2_RA_GAA) 306 Cmat2(iq,D2_GAA_GAA) = 1.368613138686131d-3*t2*wght/gammaaa* 307 1 *1.25d+0+Cmat2(iq,D2_GAA_GAA) 308 elseif (rhoa.le.tol_rho.and.rhob.gt.tol_rho) then 309 t1 = gammabb**7.5d-1 310 t2 = 1/rhob**6.666666666666666d-1 311 t3 = 1/rhob**1.6666666666666669d+0 312 t4 = 1/gammabb**2.5d-1 313 fnc(iq) = (-7.2992700729927d-3*t1*t2-9.305257363491002d-1*rh 314 1 ob**1.3333333333333333d+0)*wght+fnc(iq) 315 Amat(iq,D1_RB) = (4.8661800486617995d-3*t1*t3-1.240700981798 316 1 8002d+0*rhob**3.333333333333333d-1)*wght+Amat(iq,D1_RB) 317 Cmat(iq,D1_GBB) = Cmat(iq,D1_GBB)-5.474452554744524d-3*t2*t4 318 1 *wght 319 Amat2(iq,D2_RB_RB) = (-4.135669939329334d-1*t2-8.11030008110 320 1 3d-3*t1/rhob**2.6666666666666666d+0)*wght+Amat2(iq,D2_RB_ 321 2 RB) 322 Cmat2(iq,D2_RB_GBB) = 3.64963503649635d-3*t3*t4*wght+Cmat2(i 323 1 q,D2_RB_GBB) 324 Cmat2(iq,D2_GBB_GBB) = 1.368613138686131d-3*t2*wght/gammabb* 325 1 *1.25d+0+Cmat2(iq,D2_GBB_GBB) 326 endif ! rhoa.gt.tol_rho.and.rhob.gt.tol_rho 327 endif ! ipol.eq.1 328 enddo ! iq 329 end 330C> 331C> \brief Evaluate the nwxcm_x_gill functional [1] 332C> 333C> \f{eqnarray*}{ 334C> f &=& -{{0.0072992700729927\,\sigma_{\beta\beta}^{{{3} 335C> \over{4}}}}\over{\rho_\beta^{{{2}\over{3}}}}} 336C> -{{0.0072992700729927\,\sigma_{\alpha\alpha}^{{{3} 337C> \over{4}}}}\over{\rho_\alpha^{{{2}\over{3}}}}} 338C> -0.9305257363491002\,\rho_\beta^{{{4}\over{3}}} 339C> -0.9305257363491002\,\rho_\alpha^{{{4}\over{3}}}\\\\ 340C> g &=& 0\\\\ 341C> G &=& -{{0.0072992700729927\,\sigma_{ss}^{{{3}\over{4}}}} 342C> \over{\rho_s^{{{2}\over{3}}}}}-0.9305257363491002 343C> \,\rho_s^{{{4}\over{3}}}\\\\ 344C> \f} 345C> 346C> Code generated with Maxima 5.34.0 [2] 347C> driven by autoxc [3]. 348C> 349C> ### References ### 350C> 351C> [1] PMW Gill, Mol.Phys. 89, 433 (1996) , DOI: 352C> <a href="https://doi.org/10.1080/002689796173813 "> 353C> 10.1080/002689796173813 </a> 354C> 355C> [2] Maxima, a computer algebra system, 356C> <a href="http://maxima.sourceforge.net/"> 357C> http://maxima.sourceforge.net/</a> 358C> 359C> [3] autoxc, revision 27097 2015-05-08 360C> 361 subroutine nwxcm_x_gill_d3(param,tol_rho,ipol,nq,wght, 362 +rho,rgamma,fnc,Amat,Amat2,Amat3, 363 +Cmat,Cmat2,Cmat3) 364c $Id: $ 365#ifdef NWXC_QUAD_PREC 366 implicit real(kind=selected_real_kind(30))(a-h,o-z),integer(i-n) 367 integer, parameter :: rk=selected_real_kind(30) 368#else 369 implicit real(kind=selected_real_kind(15))(a-h,o-z),integer(i-n) 370 integer, parameter :: rk=selected_real_kind(15) 371#endif 372#include "nwxc_param.fh" 373 double precision param(*) !< [Input] Parameters of functional 374 double precision tol_rho !< [Input] The lower limit on the density 375 integer ipol !< [Input] The number of spin channels 376 integer nq !< [Input] The number of points 377 double precision wght !< [Input] The weight of the functional 378 double precision rho(nq,*) !< [Input] The density 379 double precision rgamma(nq,*) !< [Input] The norm of the density 380 !< gradients 381 double precision fnc(nq) !< [Output] The value of the functional 382c 383c Sampling Matrices for the XC Kernel 384c 385 double precision Amat(nq,*) !< [Output] The derivative wrt rho 386 double precision Cmat(nq,*) !< [Output] The derivative wrt rgamma 387c 388c Sampling Matrices for the XC Kernel 389c 390 double precision Amat2(nq,*) !< [Output] The 2nd derivative wrt rho 391 double precision Cmat2(nq,*) !< [Output] The 2nd derivative wrt rgamma 392 !< and possibly rho 393c 394c Sampling Matrices for the XC Kernel 395c 396 double precision Amat3(nq,*) !< [Output] The 3rd derivative wrt rho 397 double precision Cmat3(nq,*) !< [Output] The 3rd derivative wrt rgamma 398 !< and possibly rho 399 integer iq 400 double precision tmp 401 double precision rhoa,rhob 402 double precision gammaaa,gammaab,gammabb 403 double precision taua,taub 404 double precision nwxcm_heaviside 405 external nwxcm_heaviside 406CDIR$ NOVECTOR 407 do iq = 1, nq 408 if (ipol.eq.1) then 409 rhoa = 0.5d0*rho(iq,R_T) 410 gammaaa = 0.25d0*rgamma(iq,G_TT) 411 if (rhoa.gt.tol_rho) then 412 t1 = gammaaa**7.5d-1 413 t2 = 1/rhoa**6.666666666666666d-1 414 t3 = 1/rhoa**1.6666666666666669d+0 415 t4 = 1/gammaaa**2.5d-1 416 t5 = 1/rhoa**2.6666666666666666d+0 417 t6 = 1/gammaaa**1.25d+0 418 fnc(iq) = (-1.45985401459854d-2*t1*t2-1.8610514726982003d+0* 419 1 rhoa**1.3333333333333333d+0)*wght+fnc(iq) 420 Amat(iq,D1_RA) = (4.8661800486617995d-3*t1*t3-1.240700981798 421 1 8002d+0*rhoa**3.333333333333333d-1)*wght+Amat(iq,D1_RA) 422 Cmat(iq,D1_GAA) = Cmat(iq,D1_GAA)-5.474452554744524d-3*t2*t4 423 1 *wght 424 Cmat(iq,D1_GAB) = Cmat(iq,D1_GAB) 425 Amat2(iq,D2_RA_RA) = (-8.110300081103d-3*t1*t5-4.13566993932 426 1 9334d-1*t2)*wght+Amat2(iq,D2_RA_RA) 427 Amat2(iq,D2_RA_RB) = Amat2(iq,D2_RA_RB) 428 Cmat2(iq,D2_RA_GAA) = 3.64963503649635d-3*t3*t4*wght+Cmat2(i 429 1 q,D2_RA_GAA) 430 Cmat2(iq,D2_RA_GAB) = Cmat2(iq,D2_RA_GAB) 431 Cmat2(iq,D2_RA_GBB) = Cmat2(iq,D2_RA_GBB) 432 Cmat2(iq,D2_GAA_GAA) = 1.368613138686131d-3*t2*t6*wght+Cmat2 433 1 (iq,D2_GAA_GAA) 434 Cmat2(iq,D2_GAA_GAB) = Cmat2(iq,D2_GAA_GAB) 435 Cmat2(iq,D2_GAA_GBB) = Cmat2(iq,D2_GAA_GBB) 436 Cmat2(iq,D2_GAB_GAB) = Cmat2(iq,D2_GAB_GAB) 437 Amat3(iq,D3_RA_RA_RA) = (2.7571132928862224d-1*t3+2.16274668 438 1 82941332d-2*t1/rhoa**3.6666666666666664d+0)*wght+Amat3(iq 439 2 ,D3_RA_RA_RA) 440 Amat3(iq,D3_RA_RA_RB) = Amat3(iq,D3_RA_RA_RB) 441 Cmat3(iq,D3_RA_RA_GAA) = Cmat3(iq,D3_RA_RA_GAA)-6.0827250608 442 1 2725d-3*t4*t5*wght 443 Cmat3(iq,D3_RA_RA_GAB) = Cmat3(iq,D3_RA_RA_GAB) 444 Cmat3(iq,D3_RA_RA_GBB) = Cmat3(iq,D3_RA_RA_GBB) 445 Cmat3(iq,D3_RA_RB_GAA) = Cmat3(iq,D3_RA_RB_GAA) 446 Cmat3(iq,D3_RA_RB_GAB) = Cmat3(iq,D3_RA_RB_GAB) 447 Cmat3(iq,D3_RA_GAA_GAA) = Cmat3(iq,D3_RA_GAA_GAA)-9.12408759 448 1 1240875d-4*t3*t6*wght 449 Cmat3(iq,D3_RA_GAA_GAB) = Cmat3(iq,D3_RA_GAA_GAB) 450 Cmat3(iq,D3_RA_GAA_GBB) = Cmat3(iq,D3_RA_GAA_GBB) 451 Cmat3(iq,D3_RA_GAB_GAB) = Cmat3(iq,D3_RA_GAB_GAB) 452 Cmat3(iq,D3_RA_GAB_GBB) = Cmat3(iq,D3_RA_GAB_GBB) 453 Cmat3(iq,D3_RA_GBB_GBB) = Cmat3(iq,D3_RA_GBB_GBB) 454 Cmat3(iq,D3_GAA_GAA_GAA) = Cmat3(iq,D3_GAA_GAA_GAA)-1.710766 455 1 4233576642d-3*t2*wght/gammaaa**2.25d+0 456 Cmat3(iq,D3_GAA_GAA_GAB) = Cmat3(iq,D3_GAA_GAA_GAB) 457 Cmat3(iq,D3_GAA_GAA_GBB) = Cmat3(iq,D3_GAA_GAA_GBB) 458 Cmat3(iq,D3_GAA_GAB_GAB) = Cmat3(iq,D3_GAA_GAB_GAB) 459 Cmat3(iq,D3_GAA_GAB_GBB) = Cmat3(iq,D3_GAA_GAB_GBB) 460 Cmat3(iq,D3_GAB_GAB_GAB) = Cmat3(iq,D3_GAB_GAB_GAB) 461 endif ! rhoa.gt.tol_rho 462 else ! ipol.eq.1 463 rhoa = rho(iq,R_A) 464 rhob = rho(iq,R_B) 465 gammaaa = rgamma(iq,G_AA) 466 gammaab = rgamma(iq,G_AB) 467 gammabb = rgamma(iq,G_BB) 468 if (rhoa.gt.tol_rho.and.rhob.gt.tol_rho) then 469 t1 = gammaaa**7.5d-1 470 t2 = 1/rhoa**6.666666666666666d-1 471 t3 = gammabb**7.5d-1 472 t4 = 1/rhob**6.666666666666666d-1 473 t5 = 1/rhoa**1.6666666666666669d+0 474 t6 = 1/rhob**1.6666666666666669d+0 475 t7 = 1/gammaaa**2.5d-1 476 t8 = 1/gammabb**2.5d-1 477 t9 = 1/rhoa**2.6666666666666666d+0 478 t10 = 1/rhob**2.6666666666666666d+0 479 t11 = 1/gammaaa**1.25d+0 480 t12 = 1/gammabb**1.25d+0 481 fnc(iq) = (-7.2992700729927d-3*t3*t4-7.2992700729927d-3*t1*t 482 1 2-9.305257363491002d-1*rhob**1.3333333333333333d+0-9.3052 483 2 57363491002d-1*rhoa**1.3333333333333333d+0)*wght+fnc(iq) 484 Amat(iq,D1_RA) = (4.8661800486617995d-3*t1*t5-1.240700981798 485 1 8002d+0*rhoa**3.333333333333333d-1)*wght+Amat(iq,D1_RA) 486 Amat(iq,D1_RB) = (4.8661800486617995d-3*t3*t6-1.240700981798 487 1 8002d+0*rhob**3.333333333333333d-1)*wght+Amat(iq,D1_RB) 488 Cmat(iq,D1_GAA) = Cmat(iq,D1_GAA)-5.474452554744524d-3*t2*t7 489 1 *wght 490 Cmat(iq,D1_GAB) = Cmat(iq,D1_GAB) 491 Cmat(iq,D1_GBB) = Cmat(iq,D1_GBB)-5.474452554744524d-3*t4*t8 492 1 *wght 493 Amat2(iq,D2_RA_RA) = (-8.110300081103d-3*t1*t9-4.13566993932 494 1 9334d-1*t2)*wght+Amat2(iq,D2_RA_RA) 495 Amat2(iq,D2_RA_RB) = Amat2(iq,D2_RA_RB) 496 Amat2(iq,D2_RB_RB) = (-4.135669939329334d-1*t4-8.11030008110 497 1 3d-3*t10*t3)*wght+Amat2(iq,D2_RB_RB) 498 Cmat2(iq,D2_RA_GAA) = 3.64963503649635d-3*t5*t7*wght+Cmat2(i 499 1 q,D2_RA_GAA) 500 Cmat2(iq,D2_RA_GAB) = Cmat2(iq,D2_RA_GAB) 501 Cmat2(iq,D2_RA_GBB) = Cmat2(iq,D2_RA_GBB) 502 Cmat2(iq,D2_RB_GAA) = Cmat2(iq,D2_RB_GAA) 503 Cmat2(iq,D2_RB_GAB) = Cmat2(iq,D2_RB_GAB) 504 Cmat2(iq,D2_RB_GBB) = 3.64963503649635d-3*t6*t8*wght+Cmat2(i 505 1 q,D2_RB_GBB) 506 Cmat2(iq,D2_GAA_GAA) = 1.368613138686131d-3*t11*t2*wght+Cmat 507 1 2(iq,D2_GAA_GAA) 508 Cmat2(iq,D2_GAA_GAB) = Cmat2(iq,D2_GAA_GAB) 509 Cmat2(iq,D2_GAA_GBB) = Cmat2(iq,D2_GAA_GBB) 510 Cmat2(iq,D2_GAB_GAB) = Cmat2(iq,D2_GAB_GAB) 511 Cmat2(iq,D2_GAB_GBB) = Cmat2(iq,D2_GAB_GBB) 512 Cmat2(iq,D2_GBB_GBB) = 1.368613138686131d-3*t12*t4*wght+Cmat 513 1 2(iq,D2_GBB_GBB) 514 Amat3(iq,D3_RA_RA_RA) = (2.7571132928862224d-1*t5+2.16274668 515 1 82941332d-2*t1/rhoa**3.6666666666666664d+0)*wght+Amat3(iq 516 2 ,D3_RA_RA_RA) 517 Amat3(iq,D3_RA_RA_RB) = Amat3(iq,D3_RA_RA_RB) 518 Amat3(iq,D3_RA_RB_RB) = Amat3(iq,D3_RA_RB_RB) 519 Amat3(iq,D3_RB_RB_RB) = (2.7571132928862224d-1*t6+2.16274668 520 1 82941332d-2*t3/rhob**3.6666666666666664d+0)*wght+Amat3(iq 521 2 ,D3_RB_RB_RB) 522 Cmat3(iq,D3_RA_RA_GAA) = Cmat3(iq,D3_RA_RA_GAA)-6.0827250608 523 1 2725d-3*t7*t9*wght 524 Cmat3(iq,D3_RA_RA_GAB) = Cmat3(iq,D3_RA_RA_GAB) 525 Cmat3(iq,D3_RA_RA_GBB) = Cmat3(iq,D3_RA_RA_GBB) 526 Cmat3(iq,D3_RA_RB_GAA) = Cmat3(iq,D3_RA_RB_GAA) 527 Cmat3(iq,D3_RA_RB_GAB) = Cmat3(iq,D3_RA_RB_GAB) 528 Cmat3(iq,D3_RA_RB_GBB) = Cmat3(iq,D3_RA_RB_GBB) 529 Cmat3(iq,D3_RB_RB_GAA) = Cmat3(iq,D3_RB_RB_GAA) 530 Cmat3(iq,D3_RB_RB_GAB) = Cmat3(iq,D3_RB_RB_GAB) 531 Cmat3(iq,D3_RB_RB_GBB) = Cmat3(iq,D3_RB_RB_GBB)-6.0827250608 532 1 2725d-3*t10*t8*wght 533 Cmat3(iq,D3_RA_GAA_GAA) = Cmat3(iq,D3_RA_GAA_GAA)-9.12408759 534 1 1240875d-4*t11*t5*wght 535 Cmat3(iq,D3_RA_GAA_GAB) = Cmat3(iq,D3_RA_GAA_GAB) 536 Cmat3(iq,D3_RA_GAA_GBB) = Cmat3(iq,D3_RA_GAA_GBB) 537 Cmat3(iq,D3_RA_GAB_GAB) = Cmat3(iq,D3_RA_GAB_GAB) 538 Cmat3(iq,D3_RA_GAB_GBB) = Cmat3(iq,D3_RA_GAB_GBB) 539 Cmat3(iq,D3_RA_GBB_GBB) = Cmat3(iq,D3_RA_GBB_GBB) 540 Cmat3(iq,D3_RB_GAA_GAA) = Cmat3(iq,D3_RB_GAA_GAA) 541 Cmat3(iq,D3_RB_GAA_GAB) = Cmat3(iq,D3_RB_GAA_GAB) 542 Cmat3(iq,D3_RB_GAA_GBB) = Cmat3(iq,D3_RB_GAA_GBB) 543 Cmat3(iq,D3_RB_GAB_GAB) = Cmat3(iq,D3_RB_GAB_GAB) 544 Cmat3(iq,D3_RB_GAB_GBB) = Cmat3(iq,D3_RB_GAB_GBB) 545 Cmat3(iq,D3_RB_GBB_GBB) = Cmat3(iq,D3_RB_GBB_GBB)-9.12408759 546 1 1240875d-4*t12*t6*wght 547 Cmat3(iq,D3_GAA_GAA_GAA) = Cmat3(iq,D3_GAA_GAA_GAA)-1.710766 548 1 4233576642d-3*t2*wght/gammaaa**2.25d+0 549 Cmat3(iq,D3_GAA_GAA_GAB) = Cmat3(iq,D3_GAA_GAA_GAB) 550 Cmat3(iq,D3_GAA_GAA_GBB) = Cmat3(iq,D3_GAA_GAA_GBB) 551 Cmat3(iq,D3_GAA_GAB_GAB) = Cmat3(iq,D3_GAA_GAB_GAB) 552 Cmat3(iq,D3_GAA_GAB_GBB) = Cmat3(iq,D3_GAA_GAB_GBB) 553 Cmat3(iq,D3_GAA_GBB_GBB) = Cmat3(iq,D3_GAA_GBB_GBB) 554 Cmat3(iq,D3_GAB_GAB_GAB) = Cmat3(iq,D3_GAB_GAB_GAB) 555 Cmat3(iq,D3_GAB_GAB_GBB) = Cmat3(iq,D3_GAB_GAB_GBB) 556 Cmat3(iq,D3_GAB_GBB_GBB) = Cmat3(iq,D3_GAB_GBB_GBB) 557 Cmat3(iq,D3_GBB_GBB_GBB) = Cmat3(iq,D3_GBB_GBB_GBB)-1.710766 558 1 4233576642d-3*t4*wght/gammabb**2.25d+0 559 elseif (rhoa.gt.tol_rho.and.rhob.le.tol_rho) then 560 t1 = gammaaa**7.5d-1 561 t2 = 1/rhoa**6.666666666666666d-1 562 t3 = 1/rhoa**1.6666666666666669d+0 563 t4 = 1/gammaaa**2.5d-1 564 t5 = 1/rhoa**2.6666666666666666d+0 565 t6 = 1/gammaaa**1.25d+0 566 fnc(iq) = (-7.2992700729927d-3*t1*t2-9.305257363491002d-1*rh 567 1 oa**1.3333333333333333d+0)*wght+fnc(iq) 568 Amat(iq,D1_RA) = (4.8661800486617995d-3*t1*t3-1.240700981798 569 1 8002d+0*rhoa**3.333333333333333d-1)*wght+Amat(iq,D1_RA) 570 Cmat(iq,D1_GAA) = Cmat(iq,D1_GAA)-5.474452554744524d-3*t2*t4 571 1 *wght 572 Amat2(iq,D2_RA_RA) = (-8.110300081103d-3*t1*t5-4.13566993932 573 1 9334d-1*t2)*wght+Amat2(iq,D2_RA_RA) 574 Cmat2(iq,D2_RA_GAA) = 3.64963503649635d-3*t3*t4*wght+Cmat2(i 575 1 q,D2_RA_GAA) 576 Cmat2(iq,D2_GAA_GAA) = 1.368613138686131d-3*t2*t6*wght+Cmat2 577 1 (iq,D2_GAA_GAA) 578 Amat3(iq,D3_RA_RA_RA) = (2.7571132928862224d-1*t3+2.16274668 579 1 82941332d-2*t1/rhoa**3.6666666666666664d+0)*wght+Amat3(iq 580 2 ,D3_RA_RA_RA) 581 Cmat3(iq,D3_RA_RA_GAA) = Cmat3(iq,D3_RA_RA_GAA)-6.0827250608 582 1 2725d-3*t4*t5*wght 583 Cmat3(iq,D3_RA_GAA_GAA) = Cmat3(iq,D3_RA_GAA_GAA)-9.12408759 584 1 1240875d-4*t3*t6*wght 585 Cmat3(iq,D3_GAA_GAA_GAA) = Cmat3(iq,D3_GAA_GAA_GAA)-1.710766 586 1 4233576642d-3*t2*wght/gammaaa**2.25d+0 587 elseif (rhoa.le.tol_rho.and.rhob.gt.tol_rho) then 588 t1 = gammabb**7.5d-1 589 t2 = 1/rhob**6.666666666666666d-1 590 t3 = 1/rhob**1.6666666666666669d+0 591 t4 = 1/gammabb**2.5d-1 592 t5 = 1/rhob**2.6666666666666666d+0 593 t6 = 1/gammabb**1.25d+0 594 fnc(iq) = (-7.2992700729927d-3*t1*t2-9.305257363491002d-1*rh 595 1 ob**1.3333333333333333d+0)*wght+fnc(iq) 596 Amat(iq,D1_RB) = (4.8661800486617995d-3*t1*t3-1.240700981798 597 1 8002d+0*rhob**3.333333333333333d-1)*wght+Amat(iq,D1_RB) 598 Cmat(iq,D1_GBB) = Cmat(iq,D1_GBB)-5.474452554744524d-3*t2*t4 599 1 *wght 600 Amat2(iq,D2_RB_RB) = (-8.110300081103d-3*t1*t5-4.13566993932 601 1 9334d-1*t2)*wght+Amat2(iq,D2_RB_RB) 602 Cmat2(iq,D2_RB_GBB) = 3.64963503649635d-3*t3*t4*wght+Cmat2(i 603 1 q,D2_RB_GBB) 604 Cmat2(iq,D2_GBB_GBB) = 1.368613138686131d-3*t2*t6*wght+Cmat2 605 1 (iq,D2_GBB_GBB) 606 Amat3(iq,D3_RB_RB_RB) = (2.7571132928862224d-1*t3+2.16274668 607 1 82941332d-2*t1/rhob**3.6666666666666664d+0)*wght+Amat3(iq 608 2 ,D3_RB_RB_RB) 609 Cmat3(iq,D3_RB_RB_GBB) = Cmat3(iq,D3_RB_RB_GBB)-6.0827250608 610 1 2725d-3*t4*t5*wght 611 Cmat3(iq,D3_RB_GBB_GBB) = Cmat3(iq,D3_RB_GBB_GBB)-9.12408759 612 1 1240875d-4*t3*t6*wght 613 Cmat3(iq,D3_GBB_GBB_GBB) = Cmat3(iq,D3_GBB_GBB_GBB)-1.710766 614 1 4233576642d-3*t2*wght/gammabb**2.25d+0 615 endif ! rhoa.gt.tol_rho.and.rhob.gt.tol_rho 616 endif ! ipol.eq.1 617 enddo ! iq 618 end 619C> @} 620