1 2.. _program_listing_file__Users_robertshaw_devfiles_libecpint_src_lib_ecpint.cpp: 3 4Program Listing for File ecpint.cpp 5=================================== 6 7|exhale_lsh| :ref:`Return to documentation for file <file__Users_robertshaw_devfiles_libecpint_src_lib_ecpint.cpp>` (``/Users/robertshaw/devfiles/libecpint/src/lib/ecpint.cpp``) 8 9.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS 10 11.. code-block:: cpp 12 13 /* 14 * Copyright (c) 2020 Robert Shaw 15 * This file is a part of Libecpint. 16 * 17 * Permission is hereby granted, free of charge, to any person obtaining 18 * a copy of this software and associated documentation files (the 19 * "Software"), to deal in the Software without restriction, including 20 * without limitation the rights to use, copy, modify, merge, publish, 21 * distribute, sublicense, and/or sell copies of the Software, and to 22 * permit persons to whom the Software is furnished to do so, subject to 23 * the following conditions: 24 * 25 * The above copyright notice and this permission notice shall be 26 * included in all copies or substantial portions of the Software. 27 * 28 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, 29 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF 30 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND 31 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE 32 * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION 33 * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION 34 * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. 35 */ 36 37 #include "ecpint.hpp" 38 #include <iostream> 39 #include <cmath> 40 #include <cassert> 41 #include "Faddeeva.hpp" 42 #include "mathutil.hpp" 43 #include "qgen.hpp" 44 #include <cassert> 45 46 namespace libecpint { 47 48 ECPIntegral::ECPIntegral(const int maxLB, const int maxLU, const int deriv) { 49 // Make sure library can perform requested integrals 50 assert(maxLB+deriv <= LIBECPINT_MAX_L); 51 assert(maxLU <= LIBECPINT_MAX_L); 52 53 // Initialise singletons 54 initFactorials(); 55 zero = nonzero = skipped = 0; 56 57 // Initialise angular and radial integrators 58 angInts.init(maxLB + deriv, maxLU); 59 angInts.compute(); 60 radInts.init(2*(maxLB+deriv) + maxLU, 1e-15, 256, 512); 61 }; 62 63 double ECPIntegral::calcC(const int a, const int m, const double A) const { 64 double value = 1.0 - 2*((a-m) % 2); 65 value *= std::pow(A, a-m); 66 value *= FAC[a]/(FAC[m] * FAC[a-m]); 67 return value; 68 } 69 70 void ECPIntegral::makeC(FiveIndex<double> &C, const int L, const double *A) const { 71 int z; double Ck, Cl; 72 int na = 0; 73 for (int x = L; x >= 0; x--) { 74 for (int y = L-x; y >= 0; y--) { 75 z = L - x - y; 76 77 for (int k = 0; k<= x; k++) { 78 Ck = calcC(x, k, A[0]); 79 for (int l = 0; l <= y; l++) { 80 Cl = calcC(y, l, A[1]); 81 for (int m = 0; m <= z; m++) C(0, na, k, l, m) = Ck * Cl * calcC(z, m, A[2]); 82 } 83 } 84 na++; 85 } 86 } 87 } 88 89 void ECPIntegral::type1( 90 const ECP &U, const GaussianShell &shellA, const GaussianShell &shellB, 91 const ShellPairData &data, const FiveIndex<double> &CA, const FiveIndex<double> &CB, 92 const RadialIntegral::Parameters & parameters, TwoIndex<double> &values) const { 93 94 int LA = data.LA; int LB = data.LB; 95 int maxLBasis = data.maxLBasis; 96 97 // Build radial integrals 98 int L = LA + LB; 99 TwoIndex<double> temp; 100 ThreeIndex<double> radials(L+1, L+1, 2*L+1); 101 for (int ix = 0; ix <= L; ix++) { 102 radInts.type1(ix, ix, ix % 2, U, shellA, shellB, data, parameters, temp); 103 for(int l = 0; l <= ix; l++) { 104 for (int m = -l; m <= l; m++) radials(ix, l, l+m) = temp(l, l+m); 105 } 106 } 107 108 // Unpack positions 109 double Ax = data.A[0]; double Ay = data.A[1]; double Az = data.A[2]; 110 double Bx = data.B[0]; double By = data.B[1]; double Bz = data.B[2]; 111 112 // Calculate chi_ab for all ab in shells 113 int z1, z2, lparity, mparity, msign, ix, k, l, m; 114 double C; 115 int na = 0, nb = 0; 116 for (int x1 = LA; x1 >= 0; x1--) { 117 for (int y1 = LA-x1; y1 >= 0; y1--) { 118 z1 = LA - x1 - y1; 119 nb = 0; 120 121 for (int x2 = LB; x2 >= 0; x2--) { 122 for (int y2 = LB-x2; y2 >= 0; y2--) { 123 z2 = LB - x2 - y2; 124 125 for (int k1 = 0; k1 <= x1; k1++) { 126 for (int k2 = 0; k2 <= x2; k2++) { 127 k = k1 + k2; 128 129 for (int l1 = 0; l1 <= y1; l1++) { 130 for (int l2 = 0; l2 <= y2; l2++) { 131 l = l1 + l2; 132 133 for (int m1 = 0; m1 <= z1; m1++) { 134 for (int m2 = 0; m2 <= z2; m2++){ 135 m = m1 + m2; 136 C = CA(0, na, k1, l1, m1) * CB(0, nb, k2, l2, m2); 137 if ( fabs(C) > 1e-14 ) { 138 // Build radial integrals 139 ix = k + l + m; 140 141 // Certain terms can be neglected as the angular integrals will always be zero 142 // See Flores06 appendix for details. 143 lparity = ix % 2; 144 msign = 1 - 2*(l%2); 145 mparity = (lparity + m) % 2; 146 147 for (int lam = lparity; lam <= ix; lam+=2) { 148 for (int mu = mparity; mu <= lam; mu+=2) 149 values(na, nb) += C * angInts.getIntegral(k, l, m, lam, msign*mu) * radials(ix, lam, lam+msign*mu); 150 } 151 152 } 153 } 154 } 155 } 156 } 157 } 158 } 159 160 values(na, nb) *= 4.0 * M_PI; 161 nb++; 162 } 163 } 164 165 na++; 166 } 167 } 168 169 } 170 171 void ECPIntegral::type2( 172 const int lam, const ECP& U, const GaussianShell &shellA, const GaussianShell &shellB, 173 const ShellPairData &data, const FiveIndex<double> &CA, const FiveIndex<double> &CB, 174 const RadialIntegral::Parameters & parameters, ThreeIndex<double> &values) const { 175 176 // Unpack some data for convenience 177 int LA = data.LA; 178 int LB = data.LB; 179 int L = LA + LB; 180 int maxLBasis = data.maxLBasis; 181 182 double Am = data.Am; double Bm = data.Bm; 183 184 if (data.A_on_ecp && data.B_on_ecp) { 185 186 // Both on ECP, simplest case - see Shaw2017 supplementary material 187 double prefactor = 4.0 * M_PI; 188 int npA = shellA.nprimitive(); 189 int npB = shellB.nprimitive(); 190 int npC = U.getN(); 191 192 double zA, zB, zC, dA, dB, dC, p; 193 int nC, z1, z2; 194 195 int na = 0; 196 for (int x1 = LA; x1 >= 0; x1--) { 197 for (int r1 = LA-x1; r1 >= 0; r1--) { 198 z1 = LA - x1 - r1; 199 200 int nb = 0; 201 for (int x2 = LB; x2 >= 0; x2--) { 202 for (int y2 = LB - x2; y2 >= 0; y2--) { 203 z2 = LB - x2 - y2; 204 205 double value = 0.0; 206 for (int c = 0; c < npC; c++) { 207 const GaussianECP& g = U.getGaussian(c); 208 if (g.l == lam) { 209 zC = g.a; 210 dC = g.d; 211 nC = g.n; 212 213 for (int a = 0; a < npA; a++) { 214 zA = shellA.exp(a); 215 dA = shellA.coef(a); 216 217 for (int b = 0; b < npB; b++) { 218 zB = shellB.exp(b); 219 dB = shellB.coef(b); 220 221 p = zA + zB + zC; 222 223 double o_root_p = 1.0 / sqrt(p); 224 int N = 2 + LA + LB + nC; 225 value += 0.5*dA*dB*dC*GAMMA[N]*FAST_POW[N+1](o_root_p); 226 } 227 } 228 } 229 } 230 231 for (int mu = -lam; mu <= lam; mu++) { 232 233 double angular = prefactor * angInts.getIntegral(x1, r1, z1, lam, mu, 0, 0) * angInts.getIntegral(x2, y2, z2, lam, mu, 0, 0); 234 values(na, nb, lam+mu) = angular * value; 235 } 236 nb++; 237 } 238 } 239 240 na++; 241 } 242 } 243 244 } else { 245 246 // At least one of the shells is not on the ECP, so spherical harmonics will be required 247 248 double xA = Am > 0 ? data.A[2] / Am : 0.0; 249 double xB = Bm > 0 ? data.B[2] / Bm : 0.0; 250 double phiA = atan2(data.A[1], data.A[0]); 251 double phiB = atan2(data.B[1], data.B[0]); 252 TwoIndex<double> SA = realSphericalHarmonics(lam+LA, xA, phiA); 253 TwoIndex<double> SB = realSphericalHarmonics(lam+LB, xB, phiB); 254 255 if (data.A_on_ecp) { 256 // Radial integrals need to be calculated by a different recursive scheme, or by quadrature 257 ThreeIndex<double> radials(L+1, lam + LA + 1, lam + LB + 1); 258 TwoIndex<double> temp; 259 std::fill(values.data.begin(), values.data.end(), 0.0); 260 261 for (int N = 0; N < L+1; N++) { 262 radInts.type2(lam, 0, lam + LA, 0, lam + LB, N, U, shellA, shellB, data, parameters, temp); 263 for (int l1 = 0; l1 < lam + LA + 1; l1++) 264 for (int l2 = 0; l2 < lam + LB + 1; l2++) 265 radials(N, l1, l2) = temp(l1, l2); 266 } 267 268 // a significant number of terms can be neglected a priori - see Shaw2017 supplementary material. 269 qgen::rolled_up_special(lam, LA, LB, radials, CB, SB, angInts, values); 270 271 } else if (data.B_on_ecp){ 272 // Same as above with A and B reversed 273 ThreeIndex<double> radials(L+1, lam + LB + 1, lam + LA + 1); 274 ThreeIndex<double> tmpValues(values.dims[1], values.dims[0], values.dims[2]); 275 std::fill(tmpValues.data.begin(), tmpValues.data.end(), 0.0); 276 TwoIndex<double> temp; 277 278 for (int N = 0; N < L+1; N++) { 279 radInts.type2(lam, 0, lam + LA, 0, lam + LB, N, U, shellA, shellB, data, parameters, temp); 280 for (int l1 = 0; l1 < lam + LB + 1; l1++) 281 for (int l2 = 0; l2 < lam + LA + 1; l2++) 282 radials(N, l1, l2) = temp(l2, l1); 283 } 284 285 // a significant number of terms can be neglected a priori - see Shaw2017 supplementary material. 286 qgen::rolled_up_special(lam, LB, LA, radials, CA, SA, angInts, tmpValues); 287 // transcribe back into values 288 for (int na = 0; na < values.dims[0]; na++) 289 for (int nb = 0; nb < values.dims[1]; nb++) 290 for (int nc = 0; nc < values.dims[2]; nc++) 291 values(na, nb, nc) = tmpValues(nb, na, nc); 292 } else { 293 294 // Neither is on the ECP, the full recursive scheme with generated integrals can be used 295 // Need LA <= LB, but symmetry means we can just swap the arguments if LB > LA. 296 if (LA <= LB) 297 QGEN[LA][LB][lam](U, shellA, shellB, CA, CB, SA, SB, Am, Bm, radInts, angInts, parameters, values); 298 else { 299 ThreeIndex<double> temp_values(data.ncartB, data.ncartA, 2*U.getL() + 1); 300 QGEN[LB][LA][lam](U, shellB, shellA, CB, CA, SB, SA, Bm, Am, radInts, angInts, parameters, temp_values); 301 for (int na = 0; na < data.ncartA; na++) 302 for (int nb = 0; nb < data.ncartB; nb++) 303 for (int nu = 0; nu < 2*U.getL() + 1; nu++) 304 values(na, nb, nu) = temp_values(nb, na, nu); 305 } 306 307 } 308 } 309 } 310 311 void ECPIntegral::estimate_type2( 312 const ECP& U, const GaussianShell &shellA, const GaussianShell &shellB, 313 const ShellPairData &data, double* results) const { 314 double sigma_a, sigma_b, min_eta, n2, an, bn, a_bound, b_bound, ab_bound; 315 double atilde, btilde, ztilde, Tk, Tk_0, xp; 316 317 double Na_0 = 0.5 * data.LA / M_EULER; 318 double Nb_0 = 0.5 * data.LB / M_EULER; 319 320 for (int l = 0; l <= U.getL(); l++) { 321 min_eta = U.min_exp_l[l]; 322 n2 = min_eta * min_eta; 323 an = shellA.min_exp + min_eta; 324 bn = shellB.min_exp + min_eta; 325 if (data.A2 < 1e-6) sigma_a = 0.5 * an / shellA.min_exp; 326 else sigma_a = 0.5 * data.LA * an * an / (shellA.min_exp * (n2*data.A2 + data.LA * an)); 327 if (data.B2 < 1e-6) sigma_b = 0.5 * bn / shellB.min_exp; 328 else sigma_b = 0.5 * data.LB * bn * bn / (shellB.min_exp * (n2*data.B2 + data.LB * bn)); 329 330 atilde = (1.0 - sigma_a) * shellA.min_exp; 331 btilde = (1.0 - sigma_b) * shellB.min_exp; 332 333 a_bound = 0.0; 334 for (int i = 0; i < shellA.exps.size(); i++) 335 a_bound += FAST_POW[data.LA](std::sqrt(Na_0 / (shellA.exps[i] * sigma_a))) * std::abs(shellA.coeffs[i]); 336 337 b_bound = 0.0; 338 for (int i = 0; i < shellB.exps.size(); i++) 339 b_bound += FAST_POW[data.LB](std::sqrt(Nb_0 / (shellB.exps[i] * sigma_b))) * std::abs(shellB.coeffs[i]); 340 341 double Tk_0 = 2.0 * atilde * btilde * data.Am * data.Bm; 342 ab_bound = 0.0; 343 xp = atilde*atilde*data.A2 + btilde*btilde*data.B2; 344 for (int k = U.l_starts[l]; k < U.l_starts[l+1]; k++) { 345 const GaussianECP& g = U.getGaussian(k); 346 ztilde = atilde + btilde + g.a; 347 Tk = Tk_0 / ztilde; 348 Tk = Tk > 1 ? 0.5 * std::exp(Tk) / Tk : SINH_1; 349 ab_bound += std::abs(g.d) * FAST_POW[3](std::sqrt(M_PI/g.a)) * std::exp(xp / ztilde) * Tk; 350 } 351 ab_bound *= std::exp(-atilde*data.A2 -btilde*data.B2); 352 results[l] = (2*l+1)*(2*l+1)* a_bound * b_bound * ab_bound; 353 } 354 } 355 356 void ECPIntegral::compute_shell_pair( 357 const ECP &U, const GaussianShell &shellA, const GaussianShell &shellB, 358 TwoIndex<double> &values, const int shiftA, const int shiftB) const { 359 360 ShellPairData data; 361 362 // Shift A and B to be relative to U 363 const double* C = U.center(); 364 data.A[0] = shellA.center()[0] - C[0]; 365 data.A[1] = shellA.center()[1] - C[1]; 366 data.A[2] = shellA.center()[2] - C[2]; 367 data.B[0] = shellB.center()[0] - C[0]; 368 data.B[1] = shellB.center()[1] - C[1]; 369 data.B[2] = shellB.center()[2] - C[2]; 370 371 // Construct data that will be reused everywhere, and takes account of derivative shifts 372 data.LA = shellA.am() + shiftA; 373 data.LB = shellB.am() + shiftB; 374 data.maxLBasis = data.LA > data.LB ? data.LA : data.LB; 375 data.ncartA = (data.LA+1)*(data.LA+2)/2; 376 data.ncartB = (data.LB+1)*(data.LB+2)/2; 377 378 data.A2 = data.A[0]*data.A[0] + data.A[1]*data.A[1] + data.A[2]*data.A[2]; 379 data.Am = sqrt(data.A2); 380 data.A_on_ecp = (data.Am < 1e-6); 381 data.B2 = data.B[0]*data.B[0] + data.B[1]*data.B[1] + data.B[2]*data.B[2]; 382 data.Bm = sqrt(data.B2); 383 data.B_on_ecp = (data.Bm < 1e-6); 384 double RAB[3] = {data.A[0] - data.B[0], data.A[1] - data.B[1], data.A[2] - data.B[2]}; 385 data.RAB2 = RAB[0]*RAB[0] + RAB[1]*RAB[1] + RAB[2]*RAB[2]; 386 data.RABm = sqrt(data.RAB2); 387 388 // Prepare the radial integrator 389 const auto radIntParameters = radInts.buildParameters(shellA, shellB, data); 390 391 // Construct coefficients 392 FiveIndex<double> CA(1, data.ncartA, data.LA+1, data.LA+1, data.LA+1); 393 FiveIndex<double> CB(1, data.ncartB, data.LB+1, data.LB+1, data.LB+1); 394 makeC(CA, data.LA, data.A); 395 makeC(CB, data.LB, data.B); 396 397 double screens[U.getL() + 1]; 398 estimate_type2(U, shellA, shellB, data, screens); 399 400 // Calculate type1 integrals, if necessary 401 values.assign(data.ncartA, data.ncartB, 0.0); 402 if (!U.noType1() && screens[U.getL()] > tolerance) 403 type1(U, shellA, shellB, data, CA, CB, radIntParameters, values); 404 405 std::vector<int> l_list; 406 for (int l = 0; l < U.getL(); l++) 407 if (screens[l] > tolerance) l_list.push_back(l); 408 409 // Now all the type2 integrals 410 ThreeIndex<double> t2vals(data.ncartA, data.ncartB, 2*U.getL() + 1); 411 for (int l : l_list) { 412 t2vals.fill(0.0); 413 type2(l, U, shellA, shellB, data, CA, CB, radIntParameters, t2vals); 414 415 for (int m = -l; m <= l; m++) { 416 for(int na = 0; na < data.ncartA; na++) { 417 for (int nb = 0; nb < data.ncartB; nb++) { 418 values(na, nb) += t2vals(na, nb, l+m); 419 } 420 } 421 } 422 } 423 } 424 425 void ECPIntegral::left_shell_derivative( 426 const ECP &U, const GaussianShell &shellA, const GaussianShell &shellB, 427 std::array<TwoIndex<double>, 3> &results) const { 428 int LA = shellA.am(); 429 int LB = shellB.am(); 430 431 int ncartB = (LB+1) * (LB+2) / 2; 432 int ncartA = (LA+1) * (LA+2) / 2; 433 int ncartA_minus = LA * (LA+1) / 2; 434 TwoIndex<double> Q_minus, Q_plus; 435 436 for (auto& r : results) r.assign(ncartA, ncartB, 0.0); 437 438 if (LA != 0) 439 compute_shell_pair(U, shellA, shellB, Q_minus, -1, 0); 440 441 // hack in the exponents to the coefficients 442 GaussianShell tempA = shellA.copy(); 443 for (int i = 0; i < tempA.nprimitive(); i++) 444 tempA.coeffs[i] *= tempA.exps[i]; 445 compute_shell_pair(U, tempA, shellB, Q_plus, 1, 0); 446 447 // Now compile the derivatives 448 if (LA != 0) { 449 int nA = 0; 450 int nA_minus, nA_plus; 451 for (int k=LA; k >= 0; k--) { 452 for (int l=LA-k; l>=0; l--) { 453 int m = LA - k - l; 454 455 for (int nB = 0; nB < ncartB; nB++) { 456 nA_plus = N_INDEX(l, m); 457 nA_minus = std::min(nA_plus, Q_minus.dims[0]-1); 458 results[0](nA, nB) = -k*Q_minus(nA_minus, nB) + 2.0*Q_plus(nA_plus, nB); 459 460 nA_minus = l > 0 ? N_INDEX(l-1, m) : 0; 461 nA_plus = N_INDEX(l+1, m); 462 results[1](nA, nB) = -l*Q_minus(nA_minus, nB) + 2.0*Q_plus(nA_plus, nB); 463 464 nA_minus = m > 0 ? N_INDEX(l, m-1) : 0; 465 nA_plus = N_INDEX(l, m+1); 466 results[2](nA, nB) = -m*Q_minus(nA_minus, nB) + 2.0*Q_plus(nA_plus, nB); 467 } 468 nA += 1; 469 } 470 } 471 } else { 472 for (int nB = 0; nB < ncartB; nB++) { 473 results[0](0, nB) = 2.0*Q_plus(0, nB); 474 results[1](0, nB) = 2.0*Q_plus(1, nB); 475 results[2](0, nB) = 2.0*Q_plus(2, nB); 476 } 477 } 478 } 479 480 void ECPIntegral::left_shell_second_derivative( 481 const ECP &U, const GaussianShell &shellA, const GaussianShell &shellB, 482 std::array<TwoIndex<double>, 6> &results) const { 483 int LA = shellA.am(); 484 int LB = shellB.am(); 485 486 int ncartB = (LB+1) * (LB+2) / 2; 487 int ncartA = (LA+1) * (LA+2) / 2; 488 int ncartA_minus = std::max(1, (LA-1) * (LA) / 2); 489 TwoIndex<double> Q_minus, Q_plus, Q_0; 490 491 for (auto& r : results) r.assign(ncartA, ncartB, 0.0); 492 493 if (LA > 1) 494 compute_shell_pair(U, shellA, shellB, Q_minus, -2, 0); 495 else 496 Q_minus.assign(ncartA_minus, ncartB, 0.0); 497 498 // hack in the exponents to the coefficients 499 GaussianShell tempA = shellA.copy(); 500 for (int i = 0; i < tempA.nprimitive(); i++) 501 tempA.coeffs[i] *= tempA.exps[i]; 502 compute_shell_pair(U, tempA, shellB, Q_0, 0, 0); 503 504 // and for the l+2 505 for (int i = 0; i < tempA.nprimitive(); i++) 506 tempA.coeffs[i] *= tempA.exps[i]; 507 compute_shell_pair(U, tempA, shellB, Q_plus, 2, 0); 508 509 // Now compile the derivatives 510 int nA = 0; 511 int nA_mm, nA_pp, nA_mp, nA_pm; 512 for (int k=LA; k >= 0; k--) { 513 for (int l=LA-k; l>=0; l--) { 514 int m = LA - k - l; 515 516 for (int nB = 0; nB < ncartB; nB++) { 517 nA_mp = nA_pp = N_INDEX(l, m); //dxx 518 nA_mm = std::min(nA_mp, Q_minus.dims[0]-1); 519 results[0](nA, nB) = k*(k-1)*Q_minus(nA_mm, nB) - 2.0*(2*k+1)*Q_0(nA_mp, nB) 520 +4.0*Q_plus(nA_pp, nB); 521 522 nA_pm = l > 0 ? N_INDEX(l-1, m) : 0; 523 nA_mm = k > 0 ? nA_pm : 0; //dxy 524 nA_pp = N_INDEX(l+1, m); 525 nA_mp = k > 0 ? nA_pp : 0; 526 results[1](nA, nB) = k*l*Q_minus(nA_mm, nB) - 2.0*k*Q_0(nA_mp, nB) 527 - 2.0*l*Q_0(nA_pm, nB) + 4.0*Q_plus(nA_pp, nB); 528 529 nA_pm = m > 0 ? N_INDEX(l, m-1) : 0; 530 nA_mm = k > 0 ? nA_pm : 0; //dxz 531 nA_pp = N_INDEX(l, m+1); 532 nA_mp = k > 0 ? nA_pp : 0; 533 results[2](nA, nB) = k*m*Q_minus(nA_mm, nB) - 2.0*k*Q_0(nA_mp, nB) 534 - 2.0*m*Q_0(nA_pm, nB) + 4.0*Q_plus(nA_pp, nB); 535 536 nA_mm = l > 1 ? N_INDEX(l-2, m) : 0; //dyy 537 nA_mp = N_INDEX(l, m); 538 nA_pp = N_INDEX(l+2,m); 539 results[3](nA, nB) = l*(l-1)*Q_minus(nA_mm, nB) - 2.0*(2*l+1)*Q_0(nA_mp, nB) 540 +4.0*Q_plus(nA_pp, nB); 541 542 nA_mm = l*m > 0 ? N_INDEX(l-1, m-1) : 0; //dyz 543 nA_mp = l > 0 ? N_INDEX(l-1, m+1) : 0; 544 nA_pm = m > 0 ? N_INDEX(l+1, m-1) : 0; 545 nA_pp = N_INDEX(l+1, m+1); 546 results[4](nA, nB) = l*m*Q_minus(nA_mm, nB) - 2.0*l*Q_0(nA_mp, nB) 547 - 2.0*m*Q_0(nA_pm, nB) + 4.0*Q_plus(nA_pp, nB); 548 549 nA_mm = m > 1 ? N_INDEX(l, m-2) : 0; //dzz 550 nA_mp = N_INDEX(l, m); 551 nA_pp = N_INDEX(l,m+2); 552 results[5](nA, nB) = m*(m-1)*Q_minus(nA_mm, nB) - 2.0*(2*m+1)*Q_0(nA_mp, nB) 553 +4.0*Q_plus(nA_pp, nB); 554 555 } 556 nA += 1; 557 } 558 } 559 } 560 561 void ECPIntegral::mixed_second_derivative( 562 const ECP &U, const GaussianShell &shellA, const GaussianShell &shellB, 563 std::array<TwoIndex<double>, 9> &results) const { 564 int LA = shellA.am(); 565 int LB = shellB.am(); 566 567 int ncartB = (LB+1) * (LB+2) / 2; 568 int ncartA = (LA+1) * (LA+2) / 2; 569 int ncartB_minus = std::max(1, (LB) * (LB+1) / 2); 570 int ncartA_minus = std::max(1, (LA) * (LA+1) / 2); 571 int ncartB_plus = (LB+2) * (LB+3) / 2; 572 int ncartA_plus = (LA+2) * (LA+3) / 2; 573 TwoIndex<double> Q_mm, Q_mp, Q_pm, Q_pp; 574 575 for (auto& r : results) r.assign(ncartA, ncartB, 0.0); 576 577 GaussianShell tempA = shellA.copy(); 578 for (int i = 0; i < tempA.nprimitive(); i++) 579 tempA.coeffs[i] *= tempA.exps[i]; 580 GaussianShell tempB = shellB.copy(); 581 for (int i = 0; i < tempB.nprimitive(); i++) 582 tempB.coeffs[i] *= tempB.exps[i]; 583 584 if (LA > 0) { 585 if (LB > 0) { 586 compute_shell_pair(U, shellA, shellB, Q_mm, -1, -1); 587 compute_shell_pair(U, tempA, shellB, Q_pm, 1, -1); 588 } else { 589 Q_mm.assign(ncartA_minus, ncartB_minus, 0.0); 590 Q_pm.assign(ncartA_plus, ncartB_minus, 0.0); 591 } 592 compute_shell_pair(U, shellA, tempB, Q_mp, -1, 1); 593 } else if (LB > 0) { 594 compute_shell_pair(U, tempA, shellB, Q_pm, 1, -1); 595 Q_mm.assign(ncartA_minus, ncartB_minus, 0.0); 596 Q_mp.assign(ncartA_minus, ncartB_plus, 0.0); 597 } else { 598 Q_mm.assign(ncartA_minus, ncartB_minus, 0.0); 599 Q_mp.assign(ncartA_minus, ncartB_plus, 0.0); 600 Q_pm.assign(ncartA_plus, ncartB_minus, 0.0); 601 } 602 compute_shell_pair(U, tempA, tempB, Q_pp, 1, 1); 603 604 // Now compile the derivatives 605 int nA = 0; 606 int nB = 0; 607 int nA_m[3], nA_p[3], nB_m[3], nB_p[3], AL[3], BL[3]; 608 for (int ka=LA; ka >= 0; ka--) { 609 for (int la=LA-ka; la>=0; la--) { 610 int ma = LA - ka - la; 611 AL[0]=ka; AL[1]=la; AL[2]=ma; 612 nA_p[0] = N_INDEX(la, ma); 613 nA_m[0] = std::min(nA_p[0], Q_mm.dims[0]-1); 614 nA_m[1] = la > 0 ? N_INDEX(la-1, ma) : 0; 615 nA_m[2] = ma > 0 ? N_INDEX(la, ma-1) : 0; 616 nA_p[1] = N_INDEX(la+1,ma); 617 nA_p[2] = N_INDEX(la, ma+1); 618 619 nB = 0; 620 for (int kb=LB; kb >= 0; kb--) { 621 for (int lb=LB-kb; lb>=0; lb--) { 622 int mb = LB - kb - lb; 623 nB_p[0] = N_INDEX(lb, mb); 624 nB_m[0] = std::min(nB_p[0], Q_mm.dims[1]-1); 625 nB_m[1] = lb > 0 ? N_INDEX(lb-1, mb) : 0; 626 nB_m[2] = mb > 0 ? N_INDEX(lb, mb-1) : 0; 627 nB_p[1] = N_INDEX(lb+1,mb); 628 nB_p[2] = N_INDEX(lb, mb+1); 629 BL[0]=kb; BL[1]=lb; BL[2]=mb; 630 631 for (int p = 0; p < 3; p++) { 632 for (int q = 0; q < 3; q++) { 633 results[3*p+q](nA, nB) = AL[p]*BL[q]*Q_mm(nA_m[p], nB_m[q]) - 2.0*BL[q]*Q_pm(nA_p[p], nB_m[q]) 634 - 2.0*AL[p]*Q_mp(nA_m[p], nB_p[q]) + 4.0*Q_pp(nA_p[p], nB_p[q]); 635 } 636 } 637 638 nB += 1; 639 } 640 } 641 nA += 1; 642 } 643 } 644 } 645 646 void ECPIntegral::compute_shell_pair_derivative( 647 const ECP &U, const GaussianShell &shellA, const GaussianShell &shellB, 648 std::array<TwoIndex<double>, 9> &results) const { 649 // First we check centres 650 double A[3], B[3], C[3]; 651 for (int i = 0; i < 3; i++) { 652 A[i] = shellA.center()[i]; 653 B[i] = shellB.center()[i]; 654 C[i] = U.center()[i]; 655 } 656 657 double dAC = std::abs(A[0] - C[0]) + std::abs(A[1] - C[1]) + std::abs(A[2] - C[2]); 658 double dBC = std::abs(B[0] - C[0]) + std::abs(B[1] - C[1]) + std::abs(B[2] - C[2]); 659 660 // Calculate shell derivatives 661 std::array<TwoIndex<double>, 3> QA, QB; 662 if (dAC > 1e-6) 663 left_shell_derivative(U, shellA, shellB, QA); 664 if (dBC > 1e-6) 665 left_shell_derivative(U, shellB, shellA, QB); 666 667 // initialise results matrices 668 int ncartA = (shellA.am()+1) * (shellA.am()+2) / 2; 669 int ncartB = (shellB.am()+1) * (shellB.am()+2) / 2; 670 671 // Now construct the nuclear derivs 672 if (dAC > 1e-6) { 673 results[0] = QA[0]; 674 results[1] = QA[1]; 675 results[2] = QA[2]; 676 if (dBC > 1e-6) { 677 results[3] = QB[0].transpose(); 678 results[4] = QB[1].transpose(); 679 results[5] = QB[2].transpose(); 680 for (int i = 6; i < 9; i++) results[i].assign(ncartA, ncartB, 0.0); 681 for (int nA = 0; nA < ncartA; nA++) { 682 for (int nB = 0; nB < ncartB; nB++){ 683 results[6](nA, nB) = -1.0 * (results[0](nA, nB) + results[3](nA, nB)); 684 results[7](nA, nB) = -1.0 * (results[1](nA, nB) + results[4](nA, nB)); 685 results[8](nA, nB) = -1.0 * (results[2](nA, nB) + results[5](nA, nB)); 686 } 687 } 688 } else { 689 results[3] = results[0]; results[3].multiply(-1.0); 690 results[4] = results[1]; results[4].multiply(-1.0); 691 results[5] = results[2]; results[5].multiply(-1.0); 692 for (int i = 6; i < 9; i++) results[i].assign(ncartA, ncartB, 0.0); 693 } 694 } else if (dBC > 1e-6) { 695 results[3] = QB[0].transpose(); 696 results[4] = QB[1].transpose(); 697 results[5] = QB[2].transpose(); 698 results[0] = results[3]; results[0].multiply(-1.0); 699 results[1] = results[4]; results[1].multiply(-1.0); 700 results[2] = results[5]; results[2].multiply(-1.0); 701 for (int i = 6; i < 9; i++) results[i].assign(ncartA, ncartB, 0.0); 702 } else { 703 // else everything is zero 704 for (auto& r : results) r.assign(ncartA, ncartB, 0.0); 705 } 706 } 707 708 void ECPIntegral::compute_shell_pair_second_derivative( 709 const ECP &U, const GaussianShell &shellA, const GaussianShell &shellB, 710 std::array<TwoIndex<double>, 45> &results) const { 711 // First we check centres 712 double A[3], B[3], C[3]; 713 for (int i = 0; i < 3; i++) { 714 A[i] = shellA.center()[i]; 715 B[i] = shellB.center()[i]; 716 C[i] = U.center()[i]; 717 } 718 719 double dAC = std::abs(A[0] - C[0]) + std::abs(A[1] - C[1]) + std::abs(A[2] - C[2]); 720 double dBC = std::abs(B[0] - C[0]) + std::abs(B[1] - C[1]) + std::abs(B[2] - C[2]); 721 722 // Calculate shell derivatives 723 std::array<TwoIndex<double>, 6> QAA, QBB; 724 std::array<TwoIndex<double>, 9> QAB; 725 726 if (dAC > 1e-6) { 727 left_shell_second_derivative(U, shellA, shellB, QAA); 728 if (dBC > 1e-6) { 729 left_shell_second_derivative(U, shellB, shellA, QBB); 730 mixed_second_derivative(U, shellA, shellB, QAB); 731 } 732 } else if (dBC > 1e-6) { 733 left_shell_second_derivative(U, shellB, shellA, QBB); 734 } 735 736 // initialise results matrices 737 int ncartA = (shellA.am()+1) * (shellA.am()+2) / 2; 738 int ncartB = (shellB.am()+1) * (shellB.am()+2) / 2; 739 for (auto& r : results) r.assign(ncartA, ncartB, 0.0); 740 741 // Now construct the nuclear derivs 742 int jaas[9] = {0, 1, 2, 1, 3, 4, 2, 4, 5}; 743 int jbbs[9] = {0, 3, 6, 1, 4, 7, 2, 5, 8}; 744 int jaa, jbb; 745 if (dAC > 1e-6) { 746 //AA (xx, xy, xz, yy, yz, zz) 747 for (int i = 0; i < 6; i++) results[i] = QAA[i]; 748 749 if (dBC > 1e-6) { 750 // AB (xx, xy, xz, yx, yy, yz, zx, zy, zz) 751 for (int i = 6; i < 15; i++) results[i] = QAB[i-6]; 752 //BB (xx, xy, xz, yy, yz, zz) 753 for (int i = 24; i < 30; i++) results[i] = QBB[i-24].transpose(); 754 755 for (int nA = 0; nA < ncartA; nA++) { 756 for (int nB = 0; nB < ncartB; nB++){ 757 for (int j = 0; j < 9; j++) { 758 jaa = jaas[j]; 759 jbb = jbbs[j]; 760 761 // AC (xx, xy, xz, yx, yy, yz, zx, zy, zz) 762 results[15+j](nA, nB) = -1.0*(QAA[jaa](nA, nB) + QAB[j](nA, nB)); 763 764 // BC (xx, xy, xz, yx, yy, yz, zx, zy, zz) 765 results[30+j](nA, nB) = -1.0*(QBB[jaa](nB, nA) + QAB[jbb](nA, nB)); 766 767 // CC (xx, xy, xz, yy, yz, zz) 768 results[39+jaa](nA, nB) = -results[30+j](nA, nB) -results[15+j](nA, nB); 769 } 770 } 771 } 772 } else { 773 // AB (xx, xy, xz, yx, yy, yz, zx, zy, zz) 774 for (int i = 6; i < 15; i++) { 775 results[i] = QAA[jaas[i-6]]; 776 results[i].multiply(-1.0); 777 } 778 //BB (xx, xy, xz, yy, yz, zz) 779 for (int i = 24; i < 30; i++) results[i] = QAA[i-24]; 780 } 781 } else if (dBC > 1e-6) { 782 //BB (xx, xy, xz, yy, yz, zz) 783 for (int i = 24; i < 30; i++) results[i] = QBB[i-24].transpose(); 784 // AB (xx, xy, xz, yx, yy, yz, zx, zy, zz) 785 for (int i = 6; i < 15; i++) { 786 results[i] = QBB[jaas[i-6]].transpose(); 787 results[i].multiply(-1.0); 788 } 789 //AA (xx, xy, xz, yy, yz, zz) 790 for (int i = 0; i < 6; i++) results[i] = QBB[i].transpose(); 791 } 792 } 793 794 } 795