1 /* 2 * Copyright (c) 2020 Robert Shaw 3 * This file was generated as a part of Libecpint. 4 * 5 * Permission is hereby granted, free of charge, to any person obtaining 6 * a copy of this software and associated documentation files (the 7 * "Software"), to deal in the Software without restriction, including 8 * without limitation the rights to use, copy, modify, merge, publish, 9 * distribute, sublicense, and/or sell copies of the Software, and to 10 * permit persons to whom the Software is furnished to do so, subject to 11 * the following conditions: 12 * 13 * The above copyright notice and this permission notice shall be 14 * included in all copies or substantial portions of the Software. 15 * 16 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, 17 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF 18 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND 19 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE 20 * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION 21 * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION 22 * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. 23 */ 24 25 #include "radial.hpp" 26 #include "mathutil.hpp" 27 #include "Faddeeva.hpp" 28 #include <iostream> 29 30 namespace libecpint { 31 compute_base_integrals(const int N_min,const int N_max,const double p,const double o_root_p,const double P1,const double P2,const double P1_2,const double P2_2,const double X1,const double X2,const double oP1,const double oP2,double * values)32 void RadialIntegral::compute_base_integrals( 33 const int N_min, const int N_max, const double p, const double o_root_p, const double P1, 34 const double P2, const double P1_2, const double P2_2, const double X1, const double X2, 35 const double oP1, const double oP2, double* values) { 36 37 // Recursively construct the base integrals in order F2, G3, F4, G5, etc... as described in Shaw2017 38 39 int imax = N_max / 2; 40 int imin = (N_min + 1) / 2; 41 int gmax = (N_max - 1) / 2; 42 int gmin = N_min / 2; 43 44 double P1_2k = 1.0; 45 double P2_2k = 1.0; 46 47 for (int k = 2; k < imin; k++) { 48 P1_2k *= P1_2; 49 P2_2k *= P2_2; 50 } 51 52 double ck, dk, ek, val; 53 double C0 = o_root_p * ROOT_PI; 54 for (int n = imin; n <= imax; n++) { 55 ck = C0; 56 dk = P1_2k * X1; 57 ek = P2_2k * X2; 58 val = ck * (dk - ek); 59 60 for (int k = n - 1; k > 1; k--) { 61 ck *= 2*k*(2*k - 1)*(n-k-0.5) / ((2*n - 2*k) * (2*n - 2*k - 1) * p); 62 dk *= oP1; 63 ek *= oP2; 64 val += ck * (dk - ek); 65 } 66 67 if (n > 1) { 68 ck *= 2*(n-1.5) / ((2*n - 2) * (2*n - 3) * p); 69 val += ck * (X1 - X2); 70 } 71 72 values[2*n - N_min] = val; 73 74 P1_2k *= P1_2; 75 P2_2k *= P2_2; 76 } 77 78 P1_2k = P1; 79 P2_2k = P2; 80 for (int k = 1; k < gmin; k++) { 81 P1_2k *= P1_2; 82 P2_2k *= P2_2; 83 } 84 85 86 for (int n = gmin; n <= gmax; n++) { 87 ck = C0; 88 dk = P1_2k * X1; 89 ek = P2_2k * X2; 90 val = ck * (dk - ek); 91 92 for (int k = n-1; k >0; k--) { 93 ck *= 2*k*(2*k+1)*(n-k-0.5) / ((2*n-2*k) * (2*n - 1 - 2*k) * p); 94 dk *= oP1; 95 ek *= oP2; 96 val += ck * (dk - ek); 97 } 98 99 values[2*n + 1 - N_min] = val; 100 101 P1_2k *= P1_2; 102 P2_2k *= P2_2; 103 } 104 105 } 106 integrate_small(const int N,const int l1,const int l2,const double n,const double a,const double b,const double A,const double B)107 std::pair<double, bool> RadialIntegral::integrate_small( 108 const int N, const int l1, const int l2, const double n, 109 const double a, const double b, const double A, const double B) { 110 int gridSize = smallGrid.getN(); 111 std::vector<double> &gridPoints = smallGrid.getX(); 112 113 double Ftab[gridSize]; 114 std::vector<double> besselValues1, besselValues2; 115 116 double z, zA, zB; 117 double aA = 2.0 * a * A; 118 double bB = 2.0 * b * B; 119 for (int i = 0; i < gridSize; i++) { 120 z = gridPoints[i]; 121 zA = z - A; 122 zB = z - B; 123 124 // TODO: Efficiencies could be found here by calculating Bessel function for only l1/l2, not all l up to l1/l2 125 bessie.calculate(aA * z, l1, besselValues1); 126 bessie.calculate(bB * z, l2, besselValues2); 127 128 Ftab[i] = pow(z, N) * exp(-n * z * z - a * zA * zA - b * zB * zB) * besselValues1[l1] * besselValues2[l2]; 129 } 130 131 std::function<double(double, double*, int)> intgd = RadialIntegral::integrand; 132 133 // There should be no instances where this fails, so no backup plan to large grid, but return check just in case 134 bool success = smallGrid.integrate(intgd, Ftab, 1e-12); 135 std::pair<double, bool> rval = {smallGrid.getI(), success}; 136 return rval; 137 } 138 type2(const std::vector<Triple> & triples,const int nbase,const int lam,const ECP & U,const GaussianShell & shellA,const GaussianShell & shellB,const double A,const double B,ThreeIndex<double> & radials)139 void RadialIntegral::type2( 140 const std::vector<Triple>& triples, const int nbase, const int lam, 141 const ECP &U, const GaussianShell &shellA, const GaussianShell &shellB, 142 const double A, const double B, ThreeIndex<double> &radials) 143 { 144 int npA = shellA.nprimitive(); 145 int npB = shellB.nprimitive(); 146 147 // Loop over primitives in ECP, only considering correct ang. momentum 148 for(const auto& u : U.gaussians) { 149 if (u.l == lam) { 150 151 // Loop over primitives in orbital basis shellß 152 for(int na = 0; na < npA; na++) { 153 double a = shellA.exp(na); 154 double da = shellA.coef(na); 155 156 for (int nb = 0; nb < npB; nb++) { 157 double b = shellB.exp(nb); 158 double db = shellB.coef(nb); 159 160 // Construct values that will be reused across all radial integrals 161 double p = u.a + a + b; 162 double x = a * A; 163 double y = b * B; 164 165 double P1 = (x + y) / p; 166 double P2 = (y - x) / p; 167 double P1_2 = P1 * P1; 168 double P2_2 = P2 * P2; 169 double oP1 = 1.0 / P1_2; 170 double oP2 = std::abs(P2) < 1e-7 ? 0.0 : 1.0 / P2_2; 171 double root_p = sqrt(p); 172 double o_root_p = 1.0 / root_p; 173 double aAbB = a*A*A + b*B*B; 174 double Kab = 1.0 / (16.0 * x * y); 175 double X1 = exp(p * P1_2 - aAbB) * Kab; 176 double X2 = exp(p * P2_2 - aAbB) * Kab; 177 178 double x2 = x * x; 179 double y2 = y * y; 180 double p2 = p * p; 181 182 double result = 0.0; 183 184 // G1A, G1B may not be required, but it seems to be quicker to calculate than to check if needed 185 double daw1 = X1 * Faddeeva::Dawson(root_p * P1); 186 double daw2 = X2 * Faddeeva::Dawson(root_p * P2); 187 double G1B = 2.0 * ROOT_PI * (daw1 - daw2); 188 double G1A = 2.0 * ROOT_PI * (daw1 + daw2); 189 double H2 = ROOT_PI * ( X1 + X2 ) * o_root_p; 190 191 // Compute base integrals 192 double *values = new double[nbase+2]; 193 compute_base_integrals(2, 3+nbase, p, o_root_p, P1, P2, P1_2, P2_2, X1, X2, oP1, oP2, values); 194 195 // Loop over all radial integrals required, divert to generated code 196 for (const Triple& triple : triples ) { 197 int i = std::get<1>(triple); 198 int j = std::get<2>(triple); 199 int k = std::get<0>(triple) + u.n + 2; 200 201 int ijk = i*10000 + j*100 + k; 202 double result = 0.0; 203 if (a > 0.1 && b > 0.1) { 204 switch(ijk) { 205 case 2 : { 206 result = ( 1 ) * values[0]; 207 break; 208 } 209 210 case 4 : { 211 result += ( 1 ) * values[ 2 ]; 212 break; 213 } 214 215 case 6 : { 216 result += ( 1 ) * values[ 4 ]; 217 break; 218 } 219 220 case 8 : { 221 result += ( 1 ) * values[ 6 ]; 222 break; 223 } 224 225 case 10 : { 226 result += ( 1 ) * values[ 8 ]; 227 break; 228 } 229 230 case 12 : { 231 result += ( 1 ) * values[ 10 ]; 232 break; 233 } 234 235 case 101 : { 236 result = ( p/y ) * values[0]; 237 result += ( -x/y ) * G1A; 238 break; 239 } 240 241 case 103 : { 242 result = ( -1/(2*y) ) * values[0]; 243 result += ( 1 ) * values[ 1 ]; 244 break; 245 } 246 247 case 105 : { 248 result += ( -1/(2*y) ) * values[ 2 ]; 249 result += ( 1 ) * values[ 3 ]; 250 break; 251 } 252 253 case 107 : { 254 result += ( -1/(2*y) ) * values[ 4 ]; 255 result += ( 1 ) * values[ 5 ]; 256 break; 257 } 258 259 case 109 : { 260 result += ( -1/(2*y) ) * values[ 6 ]; 261 result += ( 1 ) * values[ 7 ]; 262 break; 263 } 264 265 case 111 : { 266 result += ( -1/(2*y) ) * values[ 8 ]; 267 result += ( 1 ) * values[ 9 ]; 268 break; 269 } 270 271 case 10102 : { 272 result = ( -(p/2 + y2)/(x*y) ) * values[0]; 273 result += ( p/x ) * values[ 1 ]; 274 break; 275 } 276 277 case 10104 : { 278 result = ( 1/(2*x*y) ) * values[0]; 279 result += ( -(p/2 + y2)/(x*y) ) * values[ 2 ]; 280 result += ( -1/x ) * values[ 1 ]; 281 result += ( p/x ) * values[ 3 ]; 282 break; 283 } 284 285 case 10106 : { 286 result += ( 1/(x*y) ) * values[ 2 ]; 287 result += ( -(p/2 + y2)/(x*y) ) * values[ 4 ]; 288 result += ( -2/x ) * values[ 3 ]; 289 result += ( p/x ) * values[ 5 ]; 290 break; 291 } 292 293 case 10108 : { 294 result += ( 3/(2*x*y) ) * values[ 4 ]; 295 result += ( -(p/2 + y2)/(x*y) ) * values[ 6 ]; 296 result += ( -3/x ) * values[ 5 ]; 297 result += ( p/x ) * values[ 7 ]; 298 break; 299 } 300 301 case 10110 : { 302 result += ( -1/(4*x*y) ) * values[ 6 ]; 303 result += ( -(p/2 + y2)/(x*y) ) * values[ 8 ]; 304 result += ( 1/(2*x) ) * values[ 7 ]; 305 result += ( p/x ) * values[ 9 ]; 306 break; 307 } 308 309 case 202 : { 310 result = ( -3*p/(2*y2) + 1 ) * values[0]; 311 result += ( 3*x/(2*y2) ) * G1A; 312 break; 313 } 314 315 case 204 : { 316 result = ( 3/(4*y2) ) * values[0]; 317 result += ( 1 ) * values[ 2 ]; 318 result += ( -3/(2*y) ) * values[ 1 ]; 319 break; 320 } 321 322 case 206 : { 323 result += ( 3/(4*y2) ) * values[ 2 ]; 324 result += ( 1 ) * values[ 4 ]; 325 result += ( -3/(2*y) ) * values[ 3 ]; 326 break; 327 } 328 329 case 208 : { 330 result += ( 3/(4*y2) ) * values[ 4 ]; 331 result += ( 1 ) * values[ 6 ]; 332 result += ( -3/(2*y) ) * values[ 5 ]; 333 break; 334 } 335 336 case 210 : { 337 result += ( 3/(4*y2) ) * values[ 6 ]; 338 result += ( 1 ) * values[ 8 ]; 339 result += ( -3/(2*y) ) * values[ 7 ]; 340 break; 341 } 342 343 case 10201 : { 344 result = ( -p*(p + 2*x2)/(2*x*y2) ) * values[0]; 345 result += ( x2/y2 ) * G1A; 346 result += ( p/y ) * H2; 347 break; 348 } 349 350 case 10203 : { 351 result = ( (3*p + 2*y2)/(4*x*y2) ) * values[0]; 352 result += ( p/x ) * values[ 2 ]; 353 result += ( -(3*p/2 + y2)/(x*y) ) * values[ 1 ]; 354 break; 355 } 356 357 case 10205 : { 358 result = ( -3/(4*x*y2) ) * values[0]; 359 result += ( (3*p - 2*y2)/(4*x*y2) ) * values[ 2 ]; 360 result += ( p/x ) * values[ 4 ]; 361 result += ( 3/(2*x*y) ) * values[ 1 ]; 362 result += ( -(3*p/2 + y2)/(x*y) ) * values[ 3 ]; 363 break; 364 } 365 366 case 10207 : { 367 result += ( -3/(2*x*y2) ) * values[ 2 ]; 368 result += ( 3*(p - 2*y2)/(4*x*y2) ) * values[ 4 ]; 369 result += ( p/x ) * values[ 6 ]; 370 result += ( 3/(x*y) ) * values[ 3 ]; 371 result += ( -(3*p/2 + y2)/(x*y) ) * values[ 5 ]; 372 break; 373 } 374 375 case 10209 : { 376 result += ( -9/(4*x*y2) ) * values[ 4 ]; 377 result += ( (3*p - 10*y2)/(4*x*y2) ) * values[ 6 ]; 378 result += ( p/x ) * values[ 8 ]; 379 result += ( 9/(2*x*y) ) * values[ 5 ]; 380 result += ( -(3*p/2 + y2)/(x*y) ) * values[ 7 ]; 381 break; 382 } 383 384 case 20202 : { 385 result = ( (3*p2 + 4*y2*(p + y2))/(4*x2*y2) ) * values[0]; 386 result += ( p2/x2 ) * values[ 2 ]; 387 result += ( -p*(3*p + 4*y2)/(2*x2*y) ) * values[ 1 ]; 388 break; 389 } 390 391 case 20204 : { 392 result = ( -(3*p/2 + y2)/(x2*y2) ) * values[0]; 393 result += ( (3*p2 + 4*y2*(-p + y2))/(4*x2*y2) ) * values[ 2 ]; 394 result += ( p2/x2 ) * values[ 4 ]; 395 result += ( (3*p + 2*y2)/(x2*y) ) * values[ 1 ]; 396 result += ( -p*(3*p + 4*y2)/(2*x2*y) ) * values[ 3 ]; 397 break; 398 } 399 400 case 20206 : { 401 result = ( 3/(2*x2*y2) ) * values[0]; 402 result += ( -3*p/(x2*y2) ) * values[ 2 ]; 403 result += ( (3*p2 + 4*y2*(-3*p + y2))/(4*x2*y2) ) * values[ 4 ]; 404 result += ( p2/x2 ) * values[ 6 ]; 405 result += ( -3/(x2*y) ) * values[ 1 ]; 406 result += ( 2*(3*p + 2*y2)/(x2*y) ) * values[ 3 ]; 407 result += ( -p*(3*p + 4*y2)/(2*x2*y) ) * values[ 5 ]; 408 break; 409 } 410 411 case 20208 : { 412 result += ( 9/(2*x2*y2) ) * values[ 2 ]; 413 result += ( 3*(-3*p + 2*y2)/(2*x2*y2) ) * values[ 4 ]; 414 result += ( (3*p2 + 4*y2*(-5*p + y2))/(4*x2*y2) ) * values[ 6 ]; 415 result += ( p2/x2 ) * values[ 8 ]; 416 result += ( -9/(x2*y) ) * values[ 3 ]; 417 result += ( 3*(3*p + 2*y2)/(x2*y) ) * values[ 5 ]; 418 result += ( -p*(3*p + 4*y2)/(2*x2*y) ) * values[ 7 ]; 419 break; 420 } 421 422 case 301 : { 423 result = ( p*(-5*p + 5*x2 + 2*y2)/(2*(y2*y)) ) * values[0]; 424 result += ( x*(15*p - 10*x2 + 6*y2)/(4*(y2*y)) ) * G1A; 425 result += ( -5*p*x/(2*y2) ) * H2; 426 break; 427 } 428 429 case 303 : { 430 result = ( 15*p/(4*(y2*y)) - 3/y ) * values[0]; 431 result += ( 1 ) * values[ 1 ]; 432 result += ( -15*x/(4*(y2*y)) ) * G1A; 433 break; 434 } 435 436 case 305 : { 437 result = ( -15/(8*(y2*y)) ) * values[0]; 438 result += ( -3/y ) * values[ 2 ]; 439 result += ( 15/(4*y2) ) * values[ 1 ]; 440 result += ( 1 ) * values[ 3 ]; 441 break; 442 } 443 444 case 307 : { 445 result += ( -15/(8*(y2*y)) ) * values[ 2 ]; 446 result += ( -3/y ) * values[ 4 ]; 447 result += ( 15/(4*y2) ) * values[ 3 ]; 448 result += ( 1 ) * values[ 5 ]; 449 break; 450 } 451 452 case 309 : { 453 result += ( -15/(8*(y2*y)) ) * values[ 4 ]; 454 result += ( -3/y ) * values[ 6 ]; 455 result += ( 15/(4*y2) ) * values[ 5 ]; 456 result += ( 1 ) * values[ 7 ]; 457 break; 458 } 459 460 case 10302 : { 461 result = ( (5*p2 + 10*p*x2 - 2*p*y2 - 4*(y2*y2))/(4*x*(y2*y)) ) * values[0]; 462 result += ( p/x ) * values[ 1 ]; 463 result += ( -5*x2/(2*(y2*y)) ) * G1A; 464 result += ( -5*p/(2*y2) ) * H2; 465 break; 466 } 467 468 case 10304 : { 469 result = ( -(15*p + 6*y2)/(8*x*(y2*y)) ) * values[0]; 470 result += ( -(3*p + y2)/(x*y) ) * values[ 2 ]; 471 result += ( 3*(5*p + 2*y2)/(4*x*y2) ) * values[ 1 ]; 472 result += ( p/x ) * values[ 3 ]; 473 break; 474 } 475 476 case 10306 : { 477 result = ( 15/(8*x*(y2*y)) ) * values[0]; 478 result += ( 3*(-5*p + 6*y2)/(8*x*(y2*y)) ) * values[ 2 ]; 479 result += ( -(3*p + y2)/(x*y) ) * values[ 4 ]; 480 result += ( -15/(4*x*y2) ) * values[ 1 ]; 481 result += ( (15*p + 2*y2)/(4*x*y2) ) * values[ 3 ]; 482 result += ( p/x ) * values[ 5 ]; 483 break; 484 } 485 486 case 10308 : { 487 result += ( 15/(4*x*(y2*y)) ) * values[ 2 ]; 488 result += ( 3*(-5*p + 14*y2)/(8*x*(y2*y)) ) * values[ 4 ]; 489 result += ( -(3*p + y2)/(x*y) ) * values[ 6 ]; 490 result += ( -15/(2*x*y2) ) * values[ 3 ]; 491 result += ( (15*p - 2*y2)/(4*x*y2) ) * values[ 5 ]; 492 result += ( p/x ) * values[ 7 ]; 493 break; 494 } 495 496 case 20301 : { 497 result = ( p*(3*p2 + 2*p*x2 + 4*(x2*x2) + 4*x2*y2 - 4*(y2*y2))/(4*x2*(y2*y)) ) * values[0]; 498 result += ( p2/x2 ) * values[ 1 ]; 499 result += ( -(x2*x)/(y2*y) ) * G1A; 500 result += ( -p*(3*p + 2*x2 + 2*y2)/(2*x*y2) ) * H2; 501 break; 502 } 503 504 case 20303 : { 505 result = ( -(15*p2 + 12*p*y2 + 4*(y2*y2))/(8*x2*(y2*y)) ) * values[0]; 506 result += ( -p*(3*p + 2*y2)/(x2*y) ) * values[ 2 ]; 507 result += ( (15*p2 + 4*y2*(3*p + y2))/(4*x2*y2) ) * values[ 1 ]; 508 result += ( p2/x2 ) * values[ 3 ]; 509 break; 510 } 511 512 case 20305 : { 513 result = ( 3*(5*p + 2*y2)/(4*x2*(y2*y)) ) * values[0]; 514 result += ( 3*(-5*p2 + 12*p*y2 + 4*(y2*y2))/(8*x2*(y2*y)) ) * values[ 2 ]; 515 result += ( -p*(3*p + 2*y2)/(x2*y) ) * values[ 4 ]; 516 result += ( -(15*p + 6*y2)/(2*x2*y2) ) * values[ 1 ]; 517 result += ( (15*p2 + 4*y2*(p + y2))/(4*x2*y2) ) * values[ 3 ]; 518 result += ( p2/x2 ) * values[ 5 ]; 519 break; 520 } 521 522 case 20307 : { 523 result = ( -15/(4*x2*(y2*y)) ) * values[0]; 524 result += ( 3*(5*p - 2*y2)/(2*x2*(y2*y)) ) * values[ 2 ]; 525 result += ( (-15*p2 + 84*p*y2 + 28*(y2*y2))/(8*x2*(y2*y)) ) * values[ 4 ]; 526 result += ( -p*(3*p + 2*y2)/(x2*y) ) * values[ 6 ]; 527 result += ( 15/(2*x2*y2) ) * values[ 1 ]; 528 result += ( -(15*p + 4*y2)/(x2*y2) ) * values[ 3 ]; 529 result += ( (15*p2 + 4*y2*(-p + y2))/(4*x2*y2) ) * values[ 5 ]; 530 result += ( p2/x2 ) * values[ 7 ]; 531 break; 532 } 533 534 case 30302 : { 535 result = ( -(15*(p2*p) + 18*p2*y2 + 12*p*(y2*y2) + 8*(y2*y2*y2))/(8*(x2*x)*(y2*y)) ) * values[0]; 536 result += ( -3*p2*(p + y2)/((x2*x)*y) ) * values[ 2 ]; 537 result += ( 3*p*(5*p2 + 6*p*y2 + 4*(y2*y2))/(4*(x2*x)*y2) ) * values[ 1 ]; 538 result += ( (p2*p)/(x2*x) ) * values[ 3 ]; 539 break; 540 } 541 542 case 30304 : { 543 result = ( 3*(15*p2 + 12*p*y2 + 4*(y2*y2))/(8*(x2*x)*(y2*y)) ) * values[0]; 544 result += ( (-15*(p2*p) + 54*p2*y2 + 4*(y2*y2)*(9*p - 2*y2))/(8*(x2*x)*(y2*y)) ) * values[ 2 ]; 545 result += ( -3*p2*(p + y2)/((x2*x)*y) ) * values[ 4 ]; 546 result += ( -(45*p2 + 12*y2*(3*p + y2))/(4*(x2*x)*y2) ) * values[ 1 ]; 547 result += ( 3*p*(5*p2 + 2*y2*(p + 2*y2))/(4*(x2*x)*y2) ) * values[ 3 ]; 548 result += ( (p2*p)/(x2*x) ) * values[ 5 ]; 549 break; 550 } 551 552 case 30306 : { 553 result = ( -(45*p + 18*y2)/(4*(x2*x)*(y2*y)) ) * values[0]; 554 result += ( 3*(15*p2 - 12*p*y2 - 4*(y2*y2))/(4*(x2*x)*(y2*y)) ) * values[ 2 ]; 555 result += ( (-15*(p2*p) + 126*p2*y2 + 4*(y2*y2)*(21*p - 2*y2))/(8*(x2*x)*(y2*y)) ) * values[ 4 ]; 556 result += ( -3*p2*(p + y2)/((x2*x)*y) ) * values[ 6 ]; 557 result += ( 9*(5*p + 2*y2)/(2*(x2*x)*y2) ) * values[ 1 ]; 558 result += ( -(45*p2 + 12*y2*(2*p + y2))/(2*(x2*x)*y2) ) * values[ 3 ]; 559 result += ( 3*p*(5*p2 + 2*y2*(-p + 2*y2))/(4*(x2*x)*y2) ) * values[ 5 ]; 560 result += ( (p2*p)/(x2*x) ) * values[ 7 ]; 561 break; 562 } 563 564 case 402 : { 565 result = ( (-5*p*(-7*p + 7*x2 + 4*y2)/4 + (y2*y2))/(y2*y2) ) * values[0]; 566 result += ( 5*x*(-21*p + 14*x2 - 6*y2)/(8*(y2*y2)) ) * G1A; 567 result += ( 35*p*x/(4*(y2*y)) ) * H2; 568 break; 569 } 570 571 case 404 : { 572 result = ( 15*(-7*p + 6*y2)/(8*(y2*y2)) ) * values[0]; 573 result += ( 1 ) * values[ 2 ]; 574 result += ( -5/y ) * values[ 1 ]; 575 result += ( 105*x/(8*(y2*y2)) ) * G1A; 576 break; 577 } 578 579 case 406 : { 580 result = ( 105/(16*(y2*y2)) ) * values[0]; 581 result += ( 45/(4*y2) ) * values[ 2 ]; 582 result += ( 1 ) * values[ 4 ]; 583 result += ( -105/(8*(y2*y)) ) * values[ 1 ]; 584 result += ( -5/y ) * values[ 3 ]; 585 break; 586 } 587 588 case 408 : { 589 result += ( 105/(16*(y2*y2)) ) * values[ 2 ]; 590 result += ( 45/(4*y2) ) * values[ 4 ]; 591 result += ( 1 ) * values[ 6 ]; 592 result += ( -105/(8*(y2*y)) ) * values[ 3 ]; 593 result += ( -5/y ) * values[ 5 ]; 594 break; 595 } 596 597 case 10401 : { 598 result = ( -p*(-7*p2 - 28*p*x2 + 2*p*y2 + 14*(x2*x2) + 4*x2*y2)/(4*x*(y2*y2)) ) * values[0]; 599 result += ( x2*(-35*p + 14*x2 - 10*y2)/(4*(y2*y2)) ) * G1A; 600 result += ( p*(-7*p + 7*x2 + 2*y2)/(2*(y2*y)) ) * H2; 601 break; 602 } 603 604 case 10403 : { 605 result = ( -(35*p2 + 70*p*x2 - 20*p*y2 - 32*(y2*y2))/(8*x*(y2*y2)) ) * values[0]; 606 result += ( p/x ) * values[ 2 ]; 607 result += ( -(5*p + y2)/(x*y) ) * values[ 1 ]; 608 result += ( 35*x2/(4*(y2*y2)) ) * G1A; 609 result += ( 35*p/(4*(y2*y)) ) * H2; 610 break; 611 } 612 613 case 10405 : { 614 result = ( 15*(7*p + 2*y2)/(16*x*(y2*y2)) ) * values[0]; 615 result += ( 45*p/(4*x*y2) + 3/x ) * values[ 2 ]; 616 result += ( p/x ) * values[ 4 ]; 617 result += ( -(105*p + 30*y2)/(8*x*(y2*y)) ) * values[ 1 ]; 618 result += ( -(5*p + y2)/(x*y) ) * values[ 3 ]; 619 break; 620 } 621 622 case 10407 : { 623 result = ( -105/(16*x*(y2*y2)) ) * values[0]; 624 result += ( 15*(7*p - 10*y2)/(16*x*(y2*y2)) ) * values[ 2 ]; 625 result += ( 45*p/(4*x*y2) + 2/x ) * values[ 4 ]; 626 result += ( p/x ) * values[ 6 ]; 627 result += ( 105/(8*x*(y2*y)) ) * values[ 1 ]; 628 result += ( 5*(-21*p + 2*y2)/(8*x*(y2*y)) ) * values[ 3 ]; 629 result += ( -(5*p + y2)/(x*y) ) * values[ 5 ]; 630 break; 631 } 632 633 case 20402 : { 634 result = ( -(21*(p2*p) + 14*p2*x2 - 6*p2*y2 + 28*p*(x2*x2) + 28*p*x2*y2 - 36*p*(y2*y2) - 8*(y2*y2*y2))/(8*x2*(y2*y2)) ) * values[0]; 635 result += ( p2/x2 ) * values[ 2 ]; 636 result += ( -p*(5*p + 2*y2)/(x2*y) ) * values[ 1 ]; 637 result += ( 7*(x2*x)/(2*(y2*y2)) ) * G1A; 638 result += ( 7*p*(3*p + 2*x2 + 2*y2)/(4*x*(y2*y)) ) * H2; 639 break; 640 } 641 642 case 20404 : { 643 result = ( 3*(35*p2 + 20*p*y2 + 4*(y2*y2))/(16*x2*(y2*y2)) ) * values[0]; 644 result += ( (45*p2 + 4*y2*(6*p + y2))/(4*x2*y2) ) * values[ 2 ]; 645 result += ( p2/x2 ) * values[ 4 ]; 646 result += ( -(105*p2 + 60*p*y2 + 12*(y2*y2))/(8*x2*(y2*y)) ) * values[ 1 ]; 647 result += ( -p*(5*p + 2*y2)/(x2*y) ) * values[ 3 ]; 648 break; 649 } 650 651 case 20406 : { 652 result = ( -(105*p + 30*y2)/(8*x2*(y2*y2)) ) * values[0]; 653 result += ( 3*(35*p2 - 100*p*y2 - 28*(y2*y2))/(16*x2*(y2*y2)) ) * values[ 2 ]; 654 result += ( (45*p2 + 4*y2*(4*p + y2))/(4*x2*y2) ) * values[ 4 ]; 655 result += ( p2/x2 ) * values[ 6 ]; 656 result += ( 15*(7*p + 2*y2)/(4*x2*(y2*y)) ) * values[ 1 ]; 657 result += ( (-105*p2 + 20*p*y2 + 4*(y2*y2))/(8*x2*(y2*y)) ) * values[ 3 ]; 658 result += ( -p*(5*p + 2*y2)/(x2*y) ) * values[ 5 ]; 659 break; 660 } 661 662 case 30401 : { 663 result = ( -p*(15*(p2*p) + 6*p2*x2 + 4*p*(x2*x2) + 24*p*x2*y2 - 36*p*(y2*y2) + 8*(x2*x2*x2) + 8*(x2*x2)*y2 + 8*x2*(y2*y2) - 16*(y2*y2*y2))/(8*(x2*x)*(y2*y2)) ) * values[0]; 664 result += ( (p2*p)/(x2*x) ) * values[ 2 ]; 665 result += ( -p2*(5*p + 3*y2)/((x2*x)*y) ) * values[ 1 ]; 666 result += ( (x2*x2)/(y2*y2) ) * G1A; 667 result += ( p*(15*p2 + 6*p*x2 + 20*p*y2 + 4*(x2*x2) + 4*x2*y2 + 4*(y2*y2))/(4*x2*(y2*y)) ) * H2; 668 break; 669 } 670 671 case 30403 : { 672 result = ( (105*(p2*p) + 90*p2*y2 + 36*p*(y2*y2) + 8*(y2*y2*y2))/(16*(x2*x)*(y2*y2)) ) * values[0]; 673 result += ( 3*p*(15*p2 + 4*y2*(3*p + y2))/(4*(x2*x)*y2) ) * values[ 2 ]; 674 result += ( (p2*p)/(x2*x) ) * values[ 4 ]; 675 result += ( -(105*(p2*p) + 90*p2*y2 + 36*p*(y2*y2) + 8*(y2*y2*y2))/(8*(x2*x)*(y2*y)) ) * values[ 1 ]; 676 result += ( -p2*(5*p + 3*y2)/((x2*x)*y) ) * values[ 3 ]; 677 break; 678 } 679 680 case 30405 : { 681 result = ( -(315*p2 + 180*p*y2 + 36*(y2*y2))/(16*(x2*x)*(y2*y2)) ) * values[0]; 682 result += ( -(-105*(p2*p) + 450*p2*y2 + 252*p*(y2*y2) + 40*(y2*y2*y2))/(16*(x2*x)*(y2*y2)) ) * values[ 2 ]; 683 result += ( 3*p*(15*p2 + 4*y2*(2*p + y2))/(4*(x2*x)*y2) ) * values[ 4 ]; 684 result += ( (p2*p)/(x2*x) ) * values[ 6 ]; 685 result += ( 9*(35*p2 + 20*p*y2 + 4*(y2*y2))/(8*(x2*x)*(y2*y)) ) * values[ 1 ]; 686 result += ( (-105*(p2*p) + 30*p2*y2 + 4*(y2*y2)*(3*p - 2*y2))/(8*(x2*x)*(y2*y)) ) * values[ 3 ]; 687 result += ( -p2*(5*p + 3*y2)/((x2*x)*y) ) * values[ 5 ]; 688 break; 689 } 690 691 case 40402 : { 692 result = ( (105*(p2*p2) + 120*(p2*p)*y2 + 72*p2*(y2*y2) + 32*p*(y2*y2*y2) + 16*(y2*y2*y2*y2))/(16*(x2*x2)*(y2*y2)) ) * values[0]; 693 result += ( 3*p2*(15*p2 + 8*y2*(2*p + y2))/(4*(x2*x2)*y2) ) * values[ 2 ]; 694 result += ( (p2*p2)/(x2*x2) ) * values[ 4 ]; 695 result += ( -p*(105*(p2*p) + 120*p2*y2 + 72*p*(y2*y2) + 32*(y2*y2*y2))/(8*(x2*x2)*(y2*y)) ) * values[ 1 ]; 696 result += ( -(p2*p)*(5*p + 4*y2)/((x2*x2)*y) ) * values[ 3 ]; 697 break; 698 } 699 700 case 40404 : { 701 result = ( -(105*(p2*p) + 90*p2*y2 + 36*p*(y2*y2) + 8*(y2*y2*y2))/(4*(x2*x2)*(y2*y2)) ) * values[0]; 702 result += ( (105*(p2*p2) - 600*(p2*p)*y2 - 504*p2*(y2*y2) - 160*p*(y2*y2*y2) + 16*(y2*y2*y2*y2))/(16*(x2*x2)*(y2*y2)) ) * values[ 2 ]; 703 result += ( p2*(45*p2 + 32*p*y2 + 24*(y2*y2))/(4*(x2*x2)*y2) ) * values[ 4 ]; 704 result += ( (p2*p2)/(x2*x2) ) * values[ 6 ]; 705 result += ( (105*(p2*p) + 90*p2*y2 + 36*p*(y2*y2) + 8*(y2*y2*y2))/(2*(x2*x2)*(y2*y)) ) * values[ 1 ]; 706 result += ( p*(-105*(p2*p) + 40*p2*y2 + (y2*y2)*(24*p - 32*y2))/(8*(x2*x2)*(y2*y)) ) * values[ 3 ]; 707 result += ( -(p2*p)*(5*p + 4*y2)/((x2*x2)*y) ) * values[ 5 ]; 708 break; 709 } 710 711 default: { 712 std::pair<double, bool> quadval = integrate_small(k, i, j, u.a, a, b, A, B); 713 result = quadval.first; 714 if (!quadval.second) std::cout << "Quadrature failed" << std::endl; 715 } 716 } 717 } else { 718 std::pair<double, bool> quadval = integrate_small(k, i, j, u.a, a, b, A, B); 719 result = quadval.first; 720 if (!quadval.second) std::cout << "Quadrature failed" << std::endl; 721 } 722 723 radials(k-2-u.n, i, j) += da * db * u.d * result; 724 } 725 726 delete[] values; 727 } 728 } 729 } 730 } 731 } 732 } 733