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