1 /* 2 * Copyright (C) 2004-2021 Edward F. Valeev 3 * 4 * This file is part of Libint. 5 * 6 * Libint is free software: you can redistribute it and/or modify 7 * it under the terms of the GNU Lesser General Public License as published by 8 * the Free Software Foundation, either version 3 of the License, or 9 * (at your option) any later version. 10 * 11 * Libint is distributed in the hope that it will be useful, 12 * but WITHOUT ANY WARRANTY; without even the implied warranty of 13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 * GNU Lesser General Public License for more details. 15 * 16 * You should have received a copy of the GNU Lesser General Public License 17 * along with Libint. If not, see <http://www.gnu.org/licenses/>. 18 * 19 */ 20 21 #include <cmath> 22 #include <vector> 23 #include <cstdlib> 24 #include <random> 25 #include <libint2.h> 26 #include <libint2/boys.h> 27 #include <time.h> 28 29 extern libint2::FmEval_Chebyshev7<double> fmeval_chebyshev; 30 extern libint2::FmEval_Taylor<double,6> fmeval_taylor; 31 32 typedef unsigned int uint; 33 template <unsigned int N> 34 struct RandomShellSet { 35 public: RandomShellSetRandomShellSet36 RandomShellSet(uint* am, 37 uint veclen, 38 uint contrdepth) { 39 40 std::copy(am, am+N, l); 41 42 std::random_device rd; 43 std::mt19937 rng(rd()); // produces randomness out of thin air 44 std::uniform_real_distribution<> rdist(0.1, 3.0); // distribution that maps to 0.1 .. 3.0 45 auto die = [&rng, &rdist]() -> double { return rdist(rng); }; // glues randomness with mapping 46 47 for(uint c=0; c<N; ++c) { 48 49 R[c].resize(3); std::generate(R[c].begin(), R[c].end(), die); 50 51 exp[c].resize(veclen); 52 coef[c].resize(veclen); 53 for(uint v=0; v<veclen; ++v) { 54 exp[c][v].resize(contrdepth); std::generate(exp[c][v].begin(), exp[c][v].end(), die); 55 coef[c][v].resize(contrdepth); std::generate(coef[c][v].begin(), coef[c][v].end(), die); 56 } 57 } 58 59 } 60 61 uint l[N]; // angular momenta 62 std::vector<double> R[N]; // origins 63 std::vector< std::vector<double> > exp[N]; // exponents 64 std::vector< std::vector<double> > coef[N]; // coefficients 65 }; 66 67 template<typename LibintEval> 68 void prep_libint2(LibintEval* erievals, 69 const RandomShellSet<4u>& rsqset, 70 int norm_flag, 71 int deriv_order = 0) { 72 const uint veclen = rsqset.exp[0].size(); 73 const uint contrdepth = rsqset.exp[0][0].size(); 74 75 const double* A = &(rsqset.R[0][0]); 76 const double* B = &(rsqset.R[1][0]); 77 const double* C = &(rsqset.R[2][0]); 78 const double* D = &(rsqset.R[3][0]); 79 80 const uint* am = rsqset.l; 81 const unsigned int am01 = am[0] + am[1]; 82 const unsigned int am23 = am[2] + am[3]; 83 const unsigned int amtot = am01 + am23 + deriv_order; 84 // static seems to be important for performance on OS X 10.7 with Apple clang 3.1 (318.0.58) 85 // 8/24/2012 also seems to affect other compilers (macports g++ 4.7.1) on OS X 10.7 and 10.8 86 static double F[LIBINT_MAX_AM*4 + 6]; 87 88 uint p0123 = 0; 89 for (uint p0 = 0; p0 < contrdepth; p0++) { 90 for (uint p1 = 0; p1 < contrdepth; p1++) { 91 for (uint p2 = 0; p2 < contrdepth; p2++) { 92 for (uint p3 = 0; p3 < contrdepth; p3++, p0123++) { 93 94 LibintEval* erieval = &erievals[p0123]; 95 erieval->veclen = veclen; 96 #if LIBINT2_FLOP_COUNT 97 erieval->nflops = erievals[0].nflops; 98 #endif 99 100 for (uint v = 0; v < veclen; v++) { 101 102 const double alpha0 = rsqset.exp[0][v][p0]; 103 const double alpha1 = rsqset.exp[1][v][p1]; 104 const double alpha2 = rsqset.exp[2][v][p2]; 105 const double alpha3 = rsqset.exp[3][v][p3]; 106 107 const double c0 = rsqset.coef[0][v][p0]; 108 const double c1 = rsqset.coef[1][v][p1]; 109 const double c2 = rsqset.coef[2][v][p2]; 110 const double c3 = rsqset.coef[3][v][p3]; 111 112 const double gammap = alpha0 + alpha1; 113 const double oogammap = 1.0 / gammap; 114 const double rhop = alpha0 * alpha1 * oogammap; 115 const double Px = (alpha0 * A[0] + alpha1 * B[0]) * oogammap; 116 const double Py = (alpha0 * A[1] + alpha1 * B[1]) * oogammap; 117 const double Pz = (alpha0 * A[2] + alpha1 * B[2]) * oogammap; 118 const double AB_x = A[0] - B[0]; 119 const double AB_y = A[1] - B[1]; 120 const double AB_z = A[2] - B[2]; 121 const double AB2 = AB_x * AB_x + AB_y * AB_y + AB_z * AB_z; 122 123 const double PAx = Px - A[0]; 124 const double PAy = Py - A[1]; 125 const double PAz = Pz - A[2]; 126 const double PBx = Px - B[0]; 127 const double PBy = Py - B[1]; 128 const double PBz = Pz - B[2]; 129 130 #if LIBINT2_DEFINED(eri,PA_x) 131 erieval->PA_x[v] = PAx; 132 #endif 133 #if LIBINT2_DEFINED(eri,PA_y) 134 erieval->PA_y[v] = PAy; 135 #endif 136 #if LIBINT2_DEFINED(eri,PA_z) 137 erieval->PA_z[v] = PAz; 138 #endif 139 #if LIBINT2_DEFINED(eri,PB_x) 140 erieval->PB_x[v] = PBx; 141 #endif 142 #if LIBINT2_DEFINED(eri,PB_y) 143 erieval->PB_y[v] = PBy; 144 #endif 145 #if LIBINT2_DEFINED(eri,PB_z) 146 erieval->PB_z[v] = PBz; 147 #endif 148 149 #if LIBINT2_DEFINED(eri,AB_x) 150 erieval->AB_x[v] = AB_x; 151 #endif 152 #if LIBINT2_DEFINED(eri,AB_y) 153 erieval->AB_y[v] = AB_y; 154 #endif 155 #if LIBINT2_DEFINED(eri,AB_z) 156 erieval->AB_z[v] = AB_z; 157 #endif 158 #if LIBINT2_DEFINED(eri,BA_x) 159 erieval->BA_x[v] = -AB_x; 160 #endif 161 #if LIBINT2_DEFINED(eri,BA_y) 162 erieval->BA_y[v] = -AB_y; 163 #endif 164 #if LIBINT2_DEFINED(eri,BA_z) 165 erieval->BA_z[v] = -AB_z; 166 #endif 167 #if LIBINT2_DEFINED(eri,oo2z) 168 erieval->oo2z[v] = 0.5*oogammap; 169 #endif 170 171 const double gammaq = alpha2 + alpha3; 172 const double oogammaq = 1.0 / gammaq; 173 const double rhoq = alpha2 * alpha3 * oogammaq; 174 const double one_o_gammap_plus_gammaq = 1.0 / (gammap + gammaq); 175 const double gammapq = gammap * gammaq * one_o_gammap_plus_gammaq; 176 const double gammap_o_gammapgammaq = gammapq * oogammaq; 177 const double gammaq_o_gammapgammaq = gammapq * oogammap; 178 const double Qx = (alpha2 * C[0] + alpha3 * D[0]) * oogammaq; 179 const double Qy = (alpha2 * C[1] + alpha3 * D[1]) * oogammaq; 180 const double Qz = (alpha2 * C[2] + alpha3 * D[2]) * oogammaq; 181 const double CD_x = C[0] - D[0]; 182 const double CD_y = C[1] - D[1]; 183 const double CD_z = C[2] - D[2]; 184 const double CD2 = CD_x * CD_x + CD_y * CD_y + CD_z * CD_z; 185 186 const double PQx = Px - Qx; 187 const double PQy = Py - Qy; 188 const double PQz = Pz - Qz; 189 const double PQ2 = PQx * PQx + PQy * PQy + PQz * PQz; 190 191 const double QCx = Qx - C[0]; 192 const double QCy = Qy - C[1]; 193 const double QCz = Qz - C[2]; 194 const double QDx = Qx - D[0]; 195 const double QDy = Qy - D[1]; 196 const double QDz = Qz - D[2]; 197 198 #if LIBINT2_DEFINED(eri,QC_x) 199 erieval->QC_x[v] = QCx; 200 #endif 201 #if LIBINT2_DEFINED(eri,QC_y) 202 erieval->QC_y[v] = QCy; 203 #endif 204 #if LIBINT2_DEFINED(eri,QC_z) 205 erieval->QC_z[v] = QCz; 206 #endif 207 #if LIBINT2_DEFINED(eri,QD_x) 208 erieval->QD_x[v] = QDx; 209 #endif 210 #if LIBINT2_DEFINED(eri,QD_y) 211 erieval->QD_y[v] = QDy; 212 #endif 213 #if LIBINT2_DEFINED(eri,QD_z) 214 erieval->QD_z[v] = QDz; 215 #endif 216 217 #if LIBINT2_DEFINED(eri,CD_x) 218 erieval->CD_x[v] = CD_x; 219 #endif 220 #if LIBINT2_DEFINED(eri,CD_y) 221 erieval->CD_y[v] = CD_y; 222 #endif 223 #if LIBINT2_DEFINED(eri,CD_z) 224 erieval->CD_z[v] = CD_z; 225 #endif 226 #if LIBINT2_DEFINED(eri,DC_x) 227 erieval->DC_x[v] = -CD_x; 228 #endif 229 #if LIBINT2_DEFINED(eri,DC_y) 230 erieval->DC_y[v] = -CD_y; 231 #endif 232 #if LIBINT2_DEFINED(eri,DC_z) 233 erieval->DC_z[v] = -CD_z; 234 #endif 235 #if LIBINT2_DEFINED(eri,oo2e) 236 erieval->oo2e[v] = 0.5*oogammaq; 237 #endif 238 239 // Prefactors for interelectron transfer relation 240 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_0_x) 241 erieval->TwoPRepITR_pfac0_0_0_x[v] = - (alpha1*AB_x + alpha3*CD_x)*oogammap; 242 #endif 243 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_0_y) 244 erieval->TwoPRepITR_pfac0_0_0_y[v] = - (alpha1*AB_y + alpha3*CD_y)*oogammap; 245 #endif 246 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_0_z) 247 erieval->TwoPRepITR_pfac0_0_0_z[v] = - (alpha1*AB_z + alpha3*CD_z)*oogammap; 248 #endif 249 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_0_x) 250 erieval->TwoPRepITR_pfac0_1_0_x[v] = - (alpha1*AB_x + alpha3*CD_x)*oogammaq; 251 #endif 252 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_0_y) 253 erieval->TwoPRepITR_pfac0_1_0_y[v] = - (alpha1*AB_y + alpha3*CD_y)*oogammaq; 254 #endif 255 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_0_z) 256 erieval->TwoPRepITR_pfac0_1_0_z[v] = - (alpha1*AB_z + alpha3*CD_z)*oogammaq; 257 #endif 258 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_1_x) 259 erieval->TwoPRepITR_pfac0_0_1_x[v] = (alpha0*AB_x + alpha2*CD_x)*oogammap; 260 #endif 261 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_1_y) 262 erieval->TwoPRepITR_pfac0_0_1_y[v] = (alpha0*AB_y + alpha2*CD_y)*oogammap; 263 #endif 264 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_1_z) 265 erieval->TwoPRepITR_pfac0_0_1_z[v] = (alpha0*AB_z + alpha2*CD_z)*oogammap; 266 #endif 267 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_1_x) 268 erieval->TwoPRepITR_pfac0_1_1_x[v] = (alpha0*AB_x + alpha2*CD_x)*oogammaq; 269 #endif 270 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_1_y) 271 erieval->TwoPRepITR_pfac0_1_1_y[v] = (alpha0*AB_y + alpha2*CD_y)*oogammaq; 272 #endif 273 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_1_z) 274 erieval->TwoPRepITR_pfac0_1_1_z[v] = (alpha0*AB_z + alpha2*CD_z)*oogammaq; 275 #endif 276 #if LIBINT2_DEFINED(eri,eoz) 277 erieval->eoz[v] = gammaq*oogammap; 278 #endif 279 #if LIBINT2_DEFINED(eri,zoe) 280 erieval->zoe[v] = gammap*oogammaq; 281 #endif 282 283 const double Wx = (gammap_o_gammapgammaq * Px + gammaq_o_gammapgammaq * Qx); 284 const double Wy = (gammap_o_gammapgammaq * Py + gammaq_o_gammapgammaq * Qy); 285 const double Wz = (gammap_o_gammapgammaq * Pz + gammaq_o_gammapgammaq * Qz); 286 287 #if LIBINT2_DEFINED(eri,WP_x) 288 erieval->WP_x[v] = Wx - Px; 289 #endif 290 #if LIBINT2_DEFINED(eri,WP_y) 291 erieval->WP_y[v] = Wy - Py; 292 #endif 293 #if LIBINT2_DEFINED(eri,WP_z) 294 erieval->WP_z[v] = Wz - Pz; 295 #endif 296 #if LIBINT2_DEFINED(eri,WQ_x) 297 erieval->WQ_x[v] = Wx - Qx; 298 #endif 299 #if LIBINT2_DEFINED(eri,WQ_y) 300 erieval->WQ_y[v] = Wy - Qy; 301 #endif 302 #if LIBINT2_DEFINED(eri,WQ_z) 303 erieval->WQ_z[v] = Wz - Qz; 304 #endif 305 #if LIBINT2_DEFINED(eri,oo2ze) 306 erieval->oo2ze[v] = 0.5/(gammap+gammaq); 307 #endif 308 #if LIBINT2_DEFINED(eri,roz) 309 erieval->roz[v] = gammapq*oogammap; 310 #endif 311 #if LIBINT2_DEFINED(eri,roe) 312 erieval->roe[v] = gammapq*oogammaq; 313 #endif 314 315 if (deriv_order > 0) { 316 // prefactors for derivative ERI relations 317 #if LIBINT2_DEFINED(eri,alpha1_rho_over_zeta2) 318 erieval->alpha1_rho_over_zeta2[v] = alpha0 * gammapq / (gammap * gammap); 319 #endif 320 #if LIBINT2_DEFINED(eri,alpha2_rho_over_zeta2) 321 erieval->alpha2_rho_over_zeta2[v] = alpha1 * gammapq / (gammap * gammap); 322 #endif 323 #if LIBINT2_DEFINED(eri,alpha3_rho_over_eta2) 324 erieval->alpha3_rho_over_eta2[v] = alpha2 * gammapq / (gammaq * gammaq); 325 #endif 326 #if LIBINT2_DEFINED(eri,alpha4_rho_over_eta2) 327 erieval->alpha4_rho_over_eta2[v] = alpha3 * gammapq / (gammaq * gammaq); 328 #endif 329 #if LIBINT2_DEFINED(eri,alpha1_over_zetapluseta) 330 erieval->alpha1_over_zetapluseta[v] = alpha0 / (gammap + gammaq); 331 #endif 332 #if LIBINT2_DEFINED(eri,alpha2_over_zetapluseta) 333 erieval->alpha2_over_zetapluseta[v] = alpha1 / (gammap + gammaq); 334 #endif 335 #if LIBINT2_DEFINED(eri,alpha3_over_zetapluseta) 336 erieval->alpha3_over_zetapluseta[v] = alpha2 / (gammap + gammaq); 337 #endif 338 #if LIBINT2_DEFINED(eri,alpha4_over_zetapluseta) 339 erieval->alpha4_over_zetapluseta[v] = alpha3 / (gammap + gammaq); 340 #endif 341 #if LIBINT2_DEFINED(eri,rho12_over_alpha1) 342 erieval->rho12_over_alpha1[v] = rhop / alpha0; 343 #endif 344 #if LIBINT2_DEFINED(eri,rho12_over_alpha2) 345 erieval->rho12_over_alpha2[v] = rhop / alpha1; 346 #endif 347 #if LIBINT2_DEFINED(eri,rho34_over_alpha3) 348 erieval->rho34_over_alpha3[v] = rhoq / alpha2; 349 #endif 350 #if LIBINT2_DEFINED(eri,rho34_over_alpha4) 351 erieval->rho34_over_alpha4[v] = rhoq / alpha3; 352 #endif 353 #if LIBINT2_DEFINED(eri,two_alpha0_bra) 354 erieval->two_alpha0_bra[v] = 2.0 * alpha0; 355 #endif 356 #if LIBINT2_DEFINED(eri,two_alpha0_ket) 357 erieval->two_alpha0_ket[v] = 2.0 * alpha1; 358 #endif 359 #if LIBINT2_DEFINED(eri,two_alpha1_bra) 360 erieval->two_alpha1_bra[v] = 2.0 * alpha2; 361 #endif 362 #if LIBINT2_DEFINED(eri,two_alpha1_ket) 363 erieval->two_alpha1_ket[v] = 2.0 * alpha3; 364 #endif 365 } 366 367 const double K1 = exp(- rhop * AB2) * oogammap; 368 const double K2 = exp(- rhoq * CD2) * oogammaq; 369 const double two_times_M_PI_to_25 = 34.986836655249725693; 370 double pfac = two_times_M_PI_to_25 * K1 * K2 * sqrt(one_o_gammap_plus_gammaq); 371 pfac *= c0 * c1 * c2 * c3; 372 373 // if veclen=1, ignore the F scratch, use the erieval directly 374 if (veclen == 1 && std::is_same<double,LIBINT2_REALTYPE>::value) { 375 { // this is only used for double realtype, so nothing nefarious here, just a workaround the type system 376 double* ssss_ptr = reinterpret_cast<double*>(erieval->LIBINT_T_SS_EREP_SS(0)); 377 fmeval_chebyshev.eval(ssss_ptr,PQ2*gammapq,amtot); 378 } 379 LIBINT2_REALTYPE* ssss_ptr = erieval->LIBINT_T_SS_EREP_SS(0); 380 for(unsigned int l=0; l<=amtot; ++l, ++ssss_ptr) 381 *ssss_ptr = *ssss_ptr * pfac; 382 } 383 else { 384 fmeval_chebyshev.eval(F,PQ2*gammapq,amtot); 385 //fmeval_taylor.eval(F,PQ2*gammapq,amtot); 386 { 387 LIBINT2_REALTYPE* ssss_ptr = erieval->LIBINT_T_SS_EREP_SS(0) + v; 388 for(unsigned int l=0; l<=amtot; ++l, ssss_ptr+=veclen) 389 *ssss_ptr = pfac*F[l]; 390 } 391 } 392 393 } // end of v loop 394 395 } 396 } 397 } 398 } // end of primitive loops 399 } 400 401 template<typename LibintEval> 402 void prep_libint2(LibintEval* erievals, 403 const RandomShellSet<3u>& rsqset, 404 int norm_flag, 405 int deriv_order = 0) { 406 const uint veclen = rsqset.exp[0].size(); 407 const uint contrdepth = rsqset.exp[0][0].size(); 408 409 const unsigned int dummy_center = (LIBINT_SHELL_SET == LIBINT_SHELL_SET_STANDARD) ? 1 : 0; 410 const double* A = &(rsqset.R[0][0]); 411 const double* B = &(rsqset.R[0][0]); 412 const double* C = &(rsqset.R[1][0]); 413 const double* D = &(rsqset.R[2][0]); 414 415 const uint* am = rsqset.l; 416 const unsigned int amtot = am[0] + am[1] + am[2] + deriv_order; 417 static double F[LIBINT_MAX_AM*3 + 6]; 418 419 uint p012 = 0; 420 for (uint p0 = 0; p0 < contrdepth; p0++) { 421 for (uint p1 = 0; p1 < contrdepth; p1++) { 422 for (uint p2 = 0; p2 < contrdepth; p2++, p012++) { 423 424 LibintEval* erieval = &erievals[p012]; 425 erieval->veclen = veclen; 426 #if LIBINT2_FLOP_COUNT 427 erieval->nflops = erievals[0].nflops; 428 #endif 429 430 for (uint v = 0; v < veclen; v++) { 431 432 const double alpha0 = (dummy_center == 0) ? 0.0 : rsqset.exp[0][v][p0]; 433 const double alpha1 = (dummy_center == 1) ? 0.0 : rsqset.exp[0][v][p0]; 434 const double alpha2 = rsqset.exp[1][v][p1]; 435 const double alpha3 = rsqset.exp[2][v][p2]; 436 437 const double c0 = (dummy_center == 0) ? 1.0 : rsqset.coef[0][v][p0]; 438 const double c1 = (dummy_center == 1) ? 1.0 : rsqset.coef[0][v][p0]; 439 const double c2 = rsqset.coef[1][v][p1]; 440 const double c3 = rsqset.coef[2][v][p2]; 441 442 const double gammap = alpha0 + alpha1; 443 const double oogammap = 1.0 / gammap; 444 const double rhop = alpha0 * alpha1 * oogammap; 445 const double Px = (alpha0 * A[0] + alpha1 * B[0]) * oogammap; 446 const double Py = (alpha0 * A[1] + alpha1 * B[1]) * oogammap; 447 const double Pz = (alpha0 * A[2] + alpha1 * B[2]) * oogammap; 448 const double PAx = Px - A[0]; 449 const double PAy = Py - A[1]; 450 const double PAz = Pz - A[2]; 451 const double PBx = Px - B[0]; 452 const double PBy = Py - B[1]; 453 const double PBz = Pz - B[2]; 454 const double AB_x = A[0] - B[0]; 455 const double AB_y = A[1] - B[1]; 456 const double AB_z = A[2] - B[2]; 457 const double AB2 = AB_x * AB_x + AB_y * AB_y + AB_z * AB_z; 458 459 #if LIBINT2_DEFINED(eri,PA_x) 460 erieval->PA_x[v] = PAx; 461 #endif 462 #if LIBINT2_DEFINED(eri,PA_y) 463 erieval->PA_y[v] = PAy; 464 #endif 465 #if LIBINT2_DEFINED(eri,PA_z) 466 erieval->PA_z[v] = PAz; 467 #endif 468 #if LIBINT2_DEFINED(eri,PB_x) 469 erieval->PB_x[v] = PBx; 470 #endif 471 #if LIBINT2_DEFINED(eri,PB_y) 472 erieval->PB_y[v] = PBy; 473 #endif 474 #if LIBINT2_DEFINED(eri,PB_z) 475 erieval->PB_z[v] = PBz; 476 #endif 477 478 #if LIBINT2_DEFINED(eri,AB_x) 479 erieval->AB_x[v] = AB_x; 480 #endif 481 #if LIBINT2_DEFINED(eri,AB_y) 482 erieval->AB_y[v] = AB_y; 483 #endif 484 #if LIBINT2_DEFINED(eri,AB_z) 485 erieval->AB_z[v] = AB_z; 486 #endif 487 #if LIBINT2_DEFINED(eri,BA_x) 488 erieval->BA_x[v] = -AB_x; 489 #endif 490 #if LIBINT2_DEFINED(eri,BA_y) 491 erieval->BA_y[v] = -AB_y; 492 #endif 493 #if LIBINT2_DEFINED(eri,BA_z) 494 erieval->BA_z[v] = -AB_z; 495 #endif 496 #if LIBINT2_DEFINED(eri,oo2z) 497 erieval->oo2z[v] = 0.5*oogammap; 498 #endif 499 500 const double gammaq = alpha2 + alpha3; 501 const double oogammaq = 1.0 / gammaq; 502 const double rhoq = alpha2 * alpha3 * oogammaq; 503 const double gammapq = gammap * gammaq / (gammap + gammaq); 504 const double gammap_o_gammapgammaq = gammapq * oogammaq; 505 const double gammaq_o_gammapgammaq = gammapq * oogammap; 506 const double Qx = (alpha2 * C[0] + alpha3 * D[0]) * oogammaq; 507 const double Qy = (alpha2 * C[1] + alpha3 * D[1]) * oogammaq; 508 const double Qz = (alpha2 * C[2] + alpha3 * D[2]) * oogammaq; 509 const double QCx = Qx - C[0]; 510 const double QCy = Qy - C[1]; 511 const double QCz = Qz - C[2]; 512 const double QDx = Qx - D[0]; 513 const double QDy = Qy - D[1]; 514 const double QDz = Qz - D[2]; 515 const double CD_x = C[0] - D[0]; 516 const double CD_y = C[1] - D[1]; 517 const double CD_z = C[2] - D[2]; 518 const double CD2 = CD_x * CD_x + CD_y * CD_y + CD_z * CD_z; 519 520 #if LIBINT2_DEFINED(eri,QC_x) 521 erieval->QC_x[v] = QCx; 522 #endif 523 #if LIBINT2_DEFINED(eri,QC_y) 524 erieval->QC_y[v] = QCy; 525 #endif 526 #if LIBINT2_DEFINED(eri,QC_z) 527 erieval->QC_z[v] = QCz; 528 #endif 529 #if LIBINT2_DEFINED(eri,QD_x) 530 erieval->QD_x[v] = QDx; 531 #endif 532 #if LIBINT2_DEFINED(eri,QD_y) 533 erieval->QD_y[v] = QDy; 534 #endif 535 #if LIBINT2_DEFINED(eri,QD_z) 536 erieval->QD_z[v] = QDz; 537 #endif 538 539 #if LIBINT2_DEFINED(eri,CD_x) 540 erieval->CD_x[v] = CD_x; 541 #endif 542 #if LIBINT2_DEFINED(eri,CD_y) 543 erieval->CD_y[v] = CD_y; 544 #endif 545 #if LIBINT2_DEFINED(eri,CD_z) 546 erieval->CD_z[v] = CD_z; 547 #endif 548 #if LIBINT2_DEFINED(eri,DC_x) 549 erieval->DC_x[v] = -CD_x; 550 #endif 551 #if LIBINT2_DEFINED(eri,DC_y) 552 erieval->DC_y[v] = -CD_y; 553 #endif 554 #if LIBINT2_DEFINED(eri,DC_z) 555 erieval->DC_z[v] = -CD_z; 556 #endif 557 #if LIBINT2_DEFINED(eri,oo2e) 558 erieval->oo2e[v] = 0.5*oogammaq; 559 #endif 560 561 // Prefactors for interelectron transfer relation 562 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_0_x) 563 erieval->TwoPRepITR_pfac0_0_0_x[v] = - (alpha1*AB_x + alpha3*CD_x)*oogammap; 564 #endif 565 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_0_y) 566 erieval->TwoPRepITR_pfac0_0_0_y[v] = - (alpha1*AB_y + alpha3*CD_y)*oogammap; 567 #endif 568 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_0_z) 569 erieval->TwoPRepITR_pfac0_0_0_z[v] = - (alpha1*AB_z + alpha3*CD_z)*oogammap; 570 #endif 571 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_0_x) 572 erieval->TwoPRepITR_pfac0_1_0_x[v] = - (alpha1*AB_x + alpha3*CD_x)*oogammaq; 573 #endif 574 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_0_y) 575 erieval->TwoPRepITR_pfac0_1_0_y[v] = - (alpha1*AB_y + alpha3*CD_y)*oogammaq; 576 #endif 577 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_0_z) 578 erieval->TwoPRepITR_pfac0_1_0_z[v] = - (alpha1*AB_z + alpha3*CD_z)*oogammaq; 579 #endif 580 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_1_x) 581 erieval->TwoPRepITR_pfac0_0_1_x[v] = (alpha0*AB_x + alpha2*CD_x)*oogammap; 582 #endif 583 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_1_y) 584 erieval->TwoPRepITR_pfac0_0_1_y[v] = (alpha0*AB_y + alpha2*CD_y)*oogammap; 585 #endif 586 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_1_z) 587 erieval->TwoPRepITR_pfac0_0_1_z[v] = (alpha0*AB_z + alpha2*CD_z)*oogammap; 588 #endif 589 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_1_x) 590 erieval->TwoPRepITR_pfac0_1_1_x[v] = (alpha0*AB_x + alpha2*CD_x)*oogammaq; 591 #endif 592 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_1_y) 593 erieval->TwoPRepITR_pfac0_1_1_y[v] = (alpha0*AB_y + alpha2*CD_y)*oogammaq; 594 #endif 595 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_1_z) 596 erieval->TwoPRepITR_pfac0_1_1_z[v] = (alpha0*AB_z + alpha2*CD_z)*oogammaq; 597 #endif 598 #if LIBINT2_DEFINED(eri,eoz) 599 erieval->eoz[v] = gammaq*oogammap; 600 #endif 601 #if LIBINT2_DEFINED(eri,zoe) 602 erieval->zoe[v] = gammap*oogammaq; 603 #endif 604 605 const double PQx = Px - Qx; 606 const double PQy = Py - Qy; 607 const double PQz = Pz - Qz; 608 const double PQ2 = PQx * PQx + PQy * PQy + PQz * PQz; 609 const double Wx = (gammap_o_gammapgammaq * Px + gammaq_o_gammapgammaq * Qx); 610 const double Wy = (gammap_o_gammapgammaq * Py + gammaq_o_gammapgammaq * Qy); 611 const double Wz = (gammap_o_gammapgammaq * Pz + gammaq_o_gammapgammaq * Qz); 612 613 #if LIBINT2_DEFINED(eri,WP_x) 614 erieval->WP_x[v] = Wx - Px; 615 #endif 616 #if LIBINT2_DEFINED(eri,WP_y) 617 erieval->WP_y[v] = Wy - Py; 618 #endif 619 #if LIBINT2_DEFINED(eri,WP_z) 620 erieval->WP_z[v] = Wz - Pz; 621 #endif 622 #if LIBINT2_DEFINED(eri,WQ_x) 623 erieval->WQ_x[v] = Wx - Qx; 624 #endif 625 #if LIBINT2_DEFINED(eri,WQ_y) 626 erieval->WQ_y[v] = Wy - Qy; 627 #endif 628 #if LIBINT2_DEFINED(eri,WQ_z) 629 erieval->WQ_z[v] = Wz - Qz; 630 #endif 631 #if LIBINT2_DEFINED(eri,oo2ze) 632 erieval->oo2ze[v] = 0.5/(gammap+gammaq); 633 #endif 634 #if LIBINT2_DEFINED(eri,roz) 635 erieval->roz[v] = gammapq*oogammap; 636 #endif 637 #if LIBINT2_DEFINED(eri,roe) 638 erieval->roe[v] = gammapq*oogammaq; 639 #endif 640 641 // prefactors for derivative ERI relations 642 if (deriv_order > 0) { 643 #if LIBINT2_DEFINED(eri,alpha1_rho_over_zeta2) 644 erieval->alpha1_rho_over_zeta2[v] = alpha0 * gammapq / (gammap * gammap); 645 #endif 646 #if LIBINT2_DEFINED(eri,alpha2_rho_over_zeta2) 647 erieval->alpha2_rho_over_zeta2[v] = alpha1 * gammapq / (gammap * gammap); 648 #endif 649 #if LIBINT2_DEFINED(eri,alpha3_rho_over_eta2) 650 erieval->alpha3_rho_over_eta2[v] = alpha2 * gammapq / (gammaq * gammaq); 651 #endif 652 #if LIBINT2_DEFINED(eri,alpha4_rho_over_eta2) 653 erieval->alpha4_rho_over_eta2[v] = alpha3 * gammapq / (gammaq * gammaq); 654 #endif 655 #if LIBINT2_DEFINED(eri,alpha1_over_zetapluseta) 656 erieval->alpha1_over_zetapluseta[v] = alpha0 / (gammap + gammaq); 657 #endif 658 #if LIBINT2_DEFINED(eri,alpha2_over_zetapluseta) 659 erieval->alpha2_over_zetapluseta[v] = alpha1 / (gammap + gammaq); 660 #endif 661 #if LIBINT2_DEFINED(eri,alpha3_over_zetapluseta) 662 erieval->alpha3_over_zetapluseta[v] = alpha2 / (gammap + gammaq); 663 #endif 664 #if LIBINT2_DEFINED(eri,alpha4_over_zetapluseta) 665 erieval->alpha4_over_zetapluseta[v] = alpha3 / (gammap + gammaq); 666 #endif 667 #if LIBINT2_DEFINED(eri,rho12_over_alpha1) 668 erieval->rho12_over_alpha1[v] = rhop / alpha0; 669 #endif 670 #if LIBINT2_DEFINED(eri,rho12_over_alpha2) 671 erieval->rho12_over_alpha2[v] = rhop / alpha1; 672 #endif 673 #if LIBINT2_DEFINED(eri,rho34_over_alpha3) 674 erieval->rho34_over_alpha3[v] = rhoq / alpha2; 675 #endif 676 #if LIBINT2_DEFINED(eri,rho34_over_alpha4) 677 erieval->rho34_over_alpha4[v] = rhoq / alpha3; 678 #endif 679 #if LIBINT2_DEFINED(eri,two_alpha0_bra) 680 erieval->two_alpha0_bra[v] = 2.0 * alpha0; 681 #endif 682 #if LIBINT2_DEFINED(eri,two_alpha0_ket) 683 erieval->two_alpha0_ket[v] = 2.0 * alpha1; 684 #endif 685 #if LIBINT2_DEFINED(eri,two_alpha1_bra) 686 erieval->two_alpha1_bra[v] = 2.0 * alpha2; 687 #endif 688 #if LIBINT2_DEFINED(eri,two_alpha1_ket) 689 erieval->two_alpha1_ket[v] = 2.0 * alpha3; 690 #endif 691 } 692 693 const double K1 = exp(- rhop * AB2); 694 const double K2 = exp(- rhoq * CD2); 695 const double two_times_M_PI_to_25 = 34.986836655249725693; 696 double pfac = two_times_M_PI_to_25 * K1 * K2 / (gammap * gammaq * sqrt(gammap 697 + gammaq)); 698 pfac *= c0 * c1 * c2 * c3; 699 700 libint2::FmEval_Reference2<double>::eval(F,PQ2*gammapq,amtot); 701 //fmeval_chebyshev.eval(F,PQ2*gammapq,amtot); 702 //fmeval_taylor.eval(F,PQ2*gammapq,amtot); 703 704 // using dangerous macros from libint2.h 705 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(0)) 706 erieval->LIBINT_T_SS_EREP_SS(0)[v] = pfac*F[0]; 707 #endif 708 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(1)) 709 erieval->LIBINT_T_SS_EREP_SS(1)[v] = pfac*F[1]; 710 #endif 711 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(2)) 712 erieval->LIBINT_T_SS_EREP_SS(2)[v] = pfac*F[2]; 713 #endif 714 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(3)) 715 erieval->LIBINT_T_SS_EREP_SS(3)[v] = pfac*F[3]; 716 #endif 717 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(4)) 718 erieval->LIBINT_T_SS_EREP_SS(4)[v] = pfac*F[4]; 719 #endif 720 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(5)) 721 erieval->LIBINT_T_SS_EREP_SS(5)[v] = pfac*F[5]; 722 #endif 723 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(6)) 724 erieval->LIBINT_T_SS_EREP_SS(6)[v] = pfac*F[6]; 725 #endif 726 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(7)) 727 erieval->LIBINT_T_SS_EREP_SS(7)[v] = pfac*F[7]; 728 #endif 729 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(8)) 730 erieval->LIBINT_T_SS_EREP_SS(8)[v] = pfac*F[8]; 731 #endif 732 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(9)) 733 erieval->LIBINT_T_SS_EREP_SS(9)[v] = pfac*F[9]; 734 #endif 735 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(10)) 736 erieval->LIBINT_T_SS_EREP_SS(10)[v] = pfac*F[10]; 737 #endif 738 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(11)) 739 erieval->LIBINT_T_SS_EREP_SS(11)[v] = pfac*F[11]; 740 #endif 741 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(12)) 742 erieval->LIBINT_T_SS_EREP_SS(12)[v] = pfac*F[12]; 743 #endif 744 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(13)) 745 erieval->LIBINT_T_SS_EREP_SS(13)[v] = pfac*F[13]; 746 #endif 747 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(14)) 748 erieval->LIBINT_T_SS_EREP_SS(14)[v] = pfac*F[14]; 749 #endif 750 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(15)) 751 erieval->LIBINT_T_SS_EREP_SS(15)[v] = pfac*F[15]; 752 #endif 753 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(16)) 754 erieval->LIBINT_T_SS_EREP_SS(16)[v] = pfac*F[16]; 755 #endif 756 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(17)) 757 erieval->LIBINT_T_SS_EREP_SS(17)[v] = pfac*F[17]; 758 #endif 759 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(18)) 760 erieval->LIBINT_T_SS_EREP_SS(18)[v] = pfac*F[18]; 761 #endif 762 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(19)) 763 erieval->LIBINT_T_SS_EREP_SS(19)[v] = pfac*F[19]; 764 #endif 765 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(20)) 766 erieval->LIBINT_T_SS_EREP_SS(20)[v] = pfac*F[20]; 767 #endif 768 } // end of v loop 769 770 } 771 } 772 } // end of primitive loops 773 } 774 775 template<typename LibintEval> 776 void prep_libint2(LibintEval* erievals, 777 const RandomShellSet<2u>& rsqset, 778 int norm_flag, 779 int deriv_order = 0) { 780 const uint veclen = rsqset.exp[0].size(); 781 const uint contrdepth = rsqset.exp[0][0].size(); 782 783 const unsigned int dummy_center1 = (LIBINT_SHELL_SET == LIBINT_SHELL_SET_STANDARD) ? 1 : 0; 784 const unsigned int dummy_center2 = (LIBINT_SHELL_SET == LIBINT_SHELL_SET_STANDARD) ? 3 : 2; 785 const double* A = &(rsqset.R[0][0]); 786 const double* B = &(rsqset.R[0][0]); 787 const double* C = &(rsqset.R[1][0]); 788 const double* D = &(rsqset.R[1][0]); 789 790 const uint* am = rsqset.l; 791 const unsigned int amtot = am[0] + am[1] + deriv_order; 792 static double F[LIBINT_MAX_AM*4 + 6]; 793 794 uint p01 = 0; 795 for (uint p0 = 0; p0 < contrdepth; p0++) { 796 for (uint p1 = 0; p1 < contrdepth; p1++, p01++) { 797 798 LibintEval* erieval = &erievals[p01]; 799 erieval->veclen = veclen; 800 #if LIBINT2_FLOP_COUNT 801 erieval->nflops = erievals[0].nflops; 802 #endif 803 804 for (uint v = 0; v < veclen; v++) { 805 806 const double alpha0 = (dummy_center1 == 0) ? 0.0 : rsqset.exp[0][v][p0]; 807 const double alpha1 = (dummy_center1 == 1) ? 0.0 : rsqset.exp[0][v][p0]; 808 const double alpha2 = (dummy_center2 == 2) ? 0.0 : rsqset.exp[1][v][p1]; 809 const double alpha3 = (dummy_center2 == 3) ? 0.0 : rsqset.exp[1][v][p1]; 810 811 const double c0 = (dummy_center1 == 0) ? 1.0 : rsqset.coef[0][v][p0]; 812 const double c1 = (dummy_center1 == 1) ? 1.0 : rsqset.coef[0][v][p0]; 813 const double c2 = (dummy_center2 == 2) ? 1.0 : rsqset.coef[1][v][p1]; 814 const double c3 = (dummy_center2 == 3) ? 1.0 : rsqset.coef[1][v][p1]; 815 816 const double gammap = alpha0 + alpha1; 817 const double oogammap = 1.0 / gammap; 818 const double rhop = alpha0 * alpha1 * oogammap; 819 const double Px = (alpha0 * A[0] + alpha1 * B[0]) * oogammap; 820 const double Py = (alpha0 * A[1] + alpha1 * B[1]) * oogammap; 821 const double Pz = (alpha0 * A[2] + alpha1 * B[2]) * oogammap; 822 const double PAx = Px - A[0]; 823 const double PAy = Py - A[1]; 824 const double PAz = Pz - A[2]; 825 const double PBx = Px - B[0]; 826 const double PBy = Py - B[1]; 827 const double PBz = Pz - B[2]; 828 const double AB_x = A[0] - B[0]; 829 const double AB_y = A[1] - B[1]; 830 const double AB_z = A[2] - B[2]; 831 const double AB2 = AB_x * AB_x + AB_y * AB_y + AB_z * AB_z; 832 833 #if LIBINT2_DEFINED(eri,PA_x) 834 erieval->PA_x[v] = PAx; 835 #endif 836 #if LIBINT2_DEFINED(eri,PA_y) 837 erieval->PA_y[v] = PAy; 838 #endif 839 #if LIBINT2_DEFINED(eri,PA_z) 840 erieval->PA_z[v] = PAz; 841 #endif 842 #if LIBINT2_DEFINED(eri,PB_x) 843 erieval->PB_x[v] = PBx; 844 #endif 845 #if LIBINT2_DEFINED(eri,PB_y) 846 erieval->PB_y[v] = PBy; 847 #endif 848 #if LIBINT2_DEFINED(eri,PB_z) 849 erieval->PB_z[v] = PBz; 850 #endif 851 852 #if LIBINT2_DEFINED(eri,AB_x) 853 erieval->AB_x[v] = AB_x; 854 #endif 855 #if LIBINT2_DEFINED(eri,AB_y) 856 erieval->AB_y[v] = AB_y; 857 #endif 858 #if LIBINT2_DEFINED(eri,AB_z) 859 erieval->AB_z[v] = AB_z; 860 #endif 861 #if LIBINT2_DEFINED(eri,oo2z) 862 erieval->oo2z[v] = 0.5*oogammap; 863 #endif 864 865 const double gammaq = alpha2 + alpha3; 866 const double oogammaq = 1.0 / gammaq; 867 const double rhoq = alpha2 * alpha3 * oogammaq; 868 const double gammapq = gammap * gammaq / (gammap + gammaq); 869 const double gammap_o_gammapgammaq = gammapq * oogammaq; 870 const double gammaq_o_gammapgammaq = gammapq * oogammap; 871 const double Qx = (alpha2 * C[0] + alpha3 * D[0]) * oogammaq; 872 const double Qy = (alpha2 * C[1] + alpha3 * D[1]) * oogammaq; 873 const double Qz = (alpha2 * C[2] + alpha3 * D[2]) * oogammaq; 874 const double QCx = Qx - C[0]; 875 const double QCy = Qy - C[1]; 876 const double QCz = Qz - C[2]; 877 const double QDx = Qx - D[0]; 878 const double QDy = Qy - D[1]; 879 const double QDz = Qz - D[2]; 880 const double CD_x = C[0] - D[0]; 881 const double CD_y = C[1] - D[1]; 882 const double CD_z = C[2] - D[2]; 883 const double CD2 = CD_x * CD_x + CD_y * CD_y + CD_z * CD_z; 884 885 #if LIBINT2_DEFINED(eri,QC_x) 886 erieval->QC_x[v] = QCx; 887 #endif 888 #if LIBINT2_DEFINED(eri,QC_y) 889 erieval->QC_y[v] = QCy; 890 #endif 891 #if LIBINT2_DEFINED(eri,QC_z) 892 erieval->QC_z[v] = QCz; 893 #endif 894 #if LIBINT2_DEFINED(eri,QD_x) 895 erieval->QD_x[v] = QDx; 896 #endif 897 #if LIBINT2_DEFINED(eri,QD_y) 898 erieval->QD_y[v] = QDy; 899 #endif 900 #if LIBINT2_DEFINED(eri,QD_z) 901 erieval->QD_z[v] = QDz; 902 #endif 903 904 #if LIBINT2_DEFINED(eri,CD_x) 905 erieval->CD_x[v] = CD_x; 906 #endif 907 #if LIBINT2_DEFINED(eri,CD_y) 908 erieval->CD_y[v] = CD_y; 909 #endif 910 #if LIBINT2_DEFINED(eri,CD_z) 911 erieval->CD_z[v] = CD_z; 912 #endif 913 #if LIBINT2_DEFINED(eri,oo2e) 914 erieval->oo2e[v] = 0.5*oogammaq; 915 #endif 916 917 // Prefactors for interelectron transfer relation 918 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_0_x) 919 erieval->TwoPRepITR_pfac0_0_0_x[v] = - (alpha1*AB_x + alpha3*CD_x)*oogammap; 920 #endif 921 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_0_y) 922 erieval->TwoPRepITR_pfac0_0_0_y[v] = - (alpha1*AB_y + alpha3*CD_y)*oogammap; 923 #endif 924 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_0_z) 925 erieval->TwoPRepITR_pfac0_0_0_z[v] = - (alpha1*AB_z + alpha3*CD_z)*oogammap; 926 #endif 927 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_0_x) 928 erieval->TwoPRepITR_pfac0_1_0_x[v] = - (alpha1*AB_x + alpha3*CD_x)*oogammaq; 929 #endif 930 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_0_y) 931 erieval->TwoPRepITR_pfac0_1_0_y[v] = - (alpha1*AB_y + alpha3*CD_y)*oogammaq; 932 #endif 933 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_0_z) 934 erieval->TwoPRepITR_pfac0_1_0_z[v] = - (alpha1*AB_z + alpha3*CD_z)*oogammaq; 935 #endif 936 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_1_x) 937 erieval->TwoPRepITR_pfac0_0_1_x[v] = (alpha0*AB_x + alpha2*CD_x)*oogammap; 938 #endif 939 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_1_y) 940 erieval->TwoPRepITR_pfac0_0_1_y[v] = (alpha0*AB_y + alpha2*CD_y)*oogammap; 941 #endif 942 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_1_z) 943 erieval->TwoPRepITR_pfac0_0_1_z[v] = (alpha0*AB_z + alpha2*CD_z)*oogammap; 944 #endif 945 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_1_x) 946 erieval->TwoPRepITR_pfac0_1_1_x[v] = (alpha0*AB_x + alpha2*CD_x)*oogammaq; 947 #endif 948 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_1_y) 949 erieval->TwoPRepITR_pfac0_1_1_y[v] = (alpha0*AB_y + alpha2*CD_y)*oogammaq; 950 #endif 951 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_1_z) 952 erieval->TwoPRepITR_pfac0_1_1_z[v] = (alpha0*AB_z + alpha2*CD_z)*oogammaq; 953 #endif 954 #if LIBINT2_DEFINED(eri,eoz) 955 erieval->eoz[v] = gammaq*oogammap; 956 #endif 957 #if LIBINT2_DEFINED(eri,zoe) 958 erieval->zoe[v] = gammap*oogammaq; 959 #endif 960 961 const double PQx = Px - Qx; 962 const double PQy = Py - Qy; 963 const double PQz = Pz - Qz; 964 const double PQ2 = PQx * PQx + PQy * PQy + PQz * PQz; 965 const double Wx = (gammap_o_gammapgammaq * Px + gammaq_o_gammapgammaq * Qx); 966 const double Wy = (gammap_o_gammapgammaq * Py + gammaq_o_gammapgammaq * Qy); 967 const double Wz = (gammap_o_gammapgammaq * Pz + gammaq_o_gammapgammaq * Qz); 968 969 #if LIBINT2_DEFINED(eri,WP_x) 970 erieval->WP_x[v] = Wx - Px; 971 #endif 972 #if LIBINT2_DEFINED(eri,WP_y) 973 erieval->WP_y[v] = Wy - Py; 974 #endif 975 #if LIBINT2_DEFINED(eri,WP_z) 976 erieval->WP_z[v] = Wz - Pz; 977 #endif 978 #if LIBINT2_DEFINED(eri,WQ_x) 979 erieval->WQ_x[v] = Wx - Qx; 980 #endif 981 #if LIBINT2_DEFINED(eri,WQ_y) 982 erieval->WQ_y[v] = Wy - Qy; 983 #endif 984 #if LIBINT2_DEFINED(eri,WQ_z) 985 erieval->WQ_z[v] = Wz - Qz; 986 #endif 987 #if LIBINT2_DEFINED(eri,oo2ze) 988 erieval->oo2ze[v] = 0.5/(gammap+gammaq); 989 #endif 990 #if LIBINT2_DEFINED(eri,roz) 991 erieval->roz[v] = gammapq*oogammap; 992 #endif 993 #if LIBINT2_DEFINED(eri,roe) 994 erieval->roe[v] = gammapq*oogammaq; 995 #endif 996 997 // prefactors for derivative ERI relations 998 if (deriv_order > 0) { 999 #if LIBINT2_DEFINED(eri,alpha1_rho_over_zeta2) 1000 erieval->alpha1_rho_over_zeta2[v] = alpha0 * gammapq / (gammap * gammap); 1001 #endif 1002 #if LIBINT2_DEFINED(eri,alpha2_rho_over_zeta2) 1003 erieval->alpha2_rho_over_zeta2[v] = alpha1 * gammapq / (gammap * gammap); 1004 #endif 1005 #if LIBINT2_DEFINED(eri,alpha3_rho_over_eta2) 1006 erieval->alpha3_rho_over_eta2[v] = alpha2 * gammapq / (gammaq * gammaq); 1007 #endif 1008 #if LIBINT2_DEFINED(eri,alpha4_rho_over_eta2) 1009 erieval->alpha4_rho_over_eta2[v] = alpha3 * gammapq / (gammaq * gammaq); 1010 #endif 1011 #if LIBINT2_DEFINED(eri,alpha1_over_zetapluseta) 1012 erieval->alpha1_over_zetapluseta[v] = alpha0 / (gammap + gammaq); 1013 #endif 1014 #if LIBINT2_DEFINED(eri,alpha2_over_zetapluseta) 1015 erieval->alpha2_over_zetapluseta[v] = alpha1 / (gammap + gammaq); 1016 #endif 1017 #if LIBINT2_DEFINED(eri,alpha3_over_zetapluseta) 1018 erieval->alpha3_over_zetapluseta[v] = alpha2 / (gammap + gammaq); 1019 #endif 1020 #if LIBINT2_DEFINED(eri,alpha4_over_zetapluseta) 1021 erieval->alpha4_over_zetapluseta[v] = alpha3 / (gammap + gammaq); 1022 #endif 1023 #if LIBINT2_DEFINED(eri,rho12_over_alpha1) 1024 erieval->rho12_over_alpha1[v] = rhop / alpha0; 1025 #endif 1026 #if LIBINT2_DEFINED(eri,rho12_over_alpha2) 1027 erieval->rho12_over_alpha2[v] = rhop / alpha1; 1028 #endif 1029 #if LIBINT2_DEFINED(eri,rho34_over_alpha3) 1030 erieval->rho34_over_alpha3[v] = rhoq / alpha2; 1031 #endif 1032 #if LIBINT2_DEFINED(eri,rho34_over_alpha4) 1033 erieval->rho34_over_alpha4[v] = rhoq / alpha3; 1034 #endif 1035 #if LIBINT2_DEFINED(eri,two_alpha0_bra) 1036 erieval->two_alpha0_bra[v] = 2.0 * alpha0; 1037 #endif 1038 #if LIBINT2_DEFINED(eri,two_alpha0_ket) 1039 erieval->two_alpha0_ket[v] = 2.0 * alpha1; 1040 #endif 1041 #if LIBINT2_DEFINED(eri,two_alpha1_bra) 1042 erieval->two_alpha1_bra[v] = 2.0 * alpha2; 1043 #endif 1044 #if LIBINT2_DEFINED(eri,two_alpha1_ket) 1045 erieval->two_alpha1_ket[v] = 2.0 * alpha3; 1046 #endif 1047 } 1048 1049 const double K1 = exp(- rhop * AB2 ); 1050 const double K2 = exp(- rhoq * CD2 ); 1051 const double two_times_M_PI_to_25 = 34.986836655249725693; 1052 double pfac = two_times_M_PI_to_25 * K1 * K2 / (gammap * gammaq * sqrt(gammap 1053 + gammaq)); 1054 pfac *= c0 * c1 * c2 * c3; 1055 1056 libint2::FmEval_Reference2<double>::eval(F,PQ2*gammapq,amtot); 1057 //fmeval_chebyshev.eval(F,PQ2*gammapq,amtot); 1058 //fmeval_taylor.eval(F,PQ2*gammapq,amtot); 1059 1060 // using dangerous macros from libint2.h 1061 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(0)) 1062 erieval->LIBINT_T_SS_EREP_SS(0)[v] = pfac*F[0]; 1063 #endif 1064 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(1)) 1065 erieval->LIBINT_T_SS_EREP_SS(1)[v] = pfac*F[1]; 1066 #endif 1067 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(2)) 1068 erieval->LIBINT_T_SS_EREP_SS(2)[v] = pfac*F[2]; 1069 #endif 1070 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(3)) 1071 erieval->LIBINT_T_SS_EREP_SS(3)[v] = pfac*F[3]; 1072 #endif 1073 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(4)) 1074 erieval->LIBINT_T_SS_EREP_SS(4)[v] = pfac*F[4]; 1075 #endif 1076 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(5)) 1077 erieval->LIBINT_T_SS_EREP_SS(5)[v] = pfac*F[5]; 1078 #endif 1079 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(6)) 1080 erieval->LIBINT_T_SS_EREP_SS(6)[v] = pfac*F[6]; 1081 #endif 1082 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(7)) 1083 erieval->LIBINT_T_SS_EREP_SS(7)[v] = pfac*F[7]; 1084 #endif 1085 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(8)) 1086 erieval->LIBINT_T_SS_EREP_SS(8)[v] = pfac*F[8]; 1087 #endif 1088 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(9)) 1089 erieval->LIBINT_T_SS_EREP_SS(9)[v] = pfac*F[9]; 1090 #endif 1091 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(10)) 1092 erieval->LIBINT_T_SS_EREP_SS(10)[v] = pfac*F[10]; 1093 #endif 1094 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(11)) 1095 erieval->LIBINT_T_SS_EREP_SS(11)[v] = pfac*F[11]; 1096 #endif 1097 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(12)) 1098 erieval->LIBINT_T_SS_EREP_SS(12)[v] = pfac*F[12]; 1099 #endif 1100 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(13)) 1101 erieval->LIBINT_T_SS_EREP_SS(13)[v] = pfac*F[13]; 1102 #endif 1103 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(14)) 1104 erieval->LIBINT_T_SS_EREP_SS(14)[v] = pfac*F[14]; 1105 #endif 1106 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(15)) 1107 erieval->LIBINT_T_SS_EREP_SS(15)[v] = pfac*F[15]; 1108 #endif 1109 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(16)) 1110 erieval->LIBINT_T_SS_EREP_SS(16)[v] = pfac*F[16]; 1111 #endif 1112 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(17)) 1113 erieval->LIBINT_T_SS_EREP_SS(17)[v] = pfac*F[17]; 1114 #endif 1115 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(18)) 1116 erieval->LIBINT_T_SS_EREP_SS(18)[v] = pfac*F[18]; 1117 #endif 1118 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(19)) 1119 erieval->LIBINT_T_SS_EREP_SS(19)[v] = pfac*F[19]; 1120 #endif 1121 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(20)) 1122 erieval->LIBINT_T_SS_EREP_SS(20)[v] = pfac*F[20]; 1123 #endif 1124 } // end of v loop 1125 1126 } 1127 } // end of primitive loops 1128 } 1129 1130