29 const AOShell& shell_col)
const {
37 if (lmax_col > 6 || lmax_row > 6) {
38 throw std::runtime_error(
39 "Orbitals higher than i are not yet implemented. This should not have "
43 const Eigen::Vector3d& pos_row = shell_row.
getPos();
44 const Eigen::Vector3d& pos_col = shell_col.
getPos();
45 const Eigen::Vector3d diff = pos_row - pos_col;
46 const double distsq = diff.squaredNorm();
48 const double kmodulus =
k_.squaredNorm();
58 Eigen::MatrixXcd cartesian = Eigen::MatrixXcd::Zero(
62 for (
const auto& gaussian_row : shell_row) {
65 const double decay_row = gaussian_row.getDecay();
67 for (
const auto& gaussian_col : shell_col) {
69 const double decay_col = gaussian_col.getDecay();
73 const double fak = 0.5 / (decay_row + decay_col);
74 const double fak2 = 2.0 * fak;
75 double exparg = fak2 * decay_row * decay_col * distsq;
78 Eigen::MatrixXcd olk = Eigen::MatrixXcd::Zero(nrows, ncols);
80 using COMPLEX = std::complex<double>;
83 PmA.real() = fak2 * (decay_row * pos_row + decay_col * pos_col) - pos_row;
84 PmA.imag() = fak *
k_;
86 PmB.real() = fak2 * (decay_row * pos_row + decay_col * pos_col) - pos_col;
87 PmB.imag() = fak *
k_;
89 const COMPLEX cfak(fak, 0.0);
90 const COMPLEX cfak2(fak2, 0.0);
93 COMPLEX ssol(std::pow(4.0 * decay_row * decay_col, 0.75) *
94 std::pow(fak2, 1.5) * std::exp(-exparg),
98 double kdotr_row =
k_.dot(pos_row);
99 double kdotr_col =
k_.dot(pos_col);
101 fak2 * (-0.25) * (kmodulus),
102 fak2 * (decay_row) * (kdotr_row) + fak2 * (decay_col) * (kdotr_col));
104 olk(0, 0) = ssol * (std::exp(kexparg));
108 olk(
Cart::x, 0) = PmA(0) * olk(0, 0);
109 olk(
Cart::y, 0) = PmA(1) * olk(0, 0);
110 olk(
Cart::z, 0) = PmA(2) * olk(0, 0);
116 COMPLEX term = (cfak) * (olk(0, 0));
141 COMPLEX term_xx = (cfak) * (olk(
Cart::xx, 0));
142 COMPLEX term_yy = (cfak) * (olk(
Cart::yy, 0));
143 COMPLEX term_zz = (cfak) * (olk(
Cart::zz, 0));
163 COMPLEX term_xxx = (cfak) * (olk(
Cart::xxx, 0));
164 COMPLEX term_yyy = (cfak) * (olk(
Cart::yyy, 0));
165 COMPLEX term_zzz = (cfak) * (olk(
Cart::zzz, 0));
191 COMPLEX term_xxxx = (cfak) * (olk(
Cart::xxxx, 0));
192 COMPLEX term_xyyy = (cfak) * (olk(
Cart::xyyy, 0));
193 COMPLEX term_xzzz = (cfak) * (olk(
Cart::xzzz, 0));
194 COMPLEX term_yyyy = (cfak) * (olk(
Cart::yyyy, 0));
195 COMPLEX term_yyzz = (cfak) * (olk(
Cart::yyzz, 0));
196 COMPLEX term_yzzz = (cfak) * (olk(
Cart::yzzz, 0));
197 COMPLEX term_zzzz = (cfak) * (olk(
Cart::zzzz, 0));
231 olk(0,
Cart::x) = PmB(0) * olk(0, 0);
232 olk(0,
Cart::y) = PmB(1) * olk(0, 0);
233 olk(0,
Cart::z) = PmB(2) * olk(0, 0);
238 for (
Index i = 1; i < n_orbitals[lmax_row]; i++) {
241 PmB(0) * olk(i, 0) + double(nx[i]) * cfak * olk(i_less_x[i], 0);
244 PmB(1) * olk(i, 0) + double(ny[i]) * cfak * olk(i_less_y[i], 0);
247 PmB(2) * olk(i, 0) + double(nz[i]) * cfak * olk(i_less_z[i], 0);
255 COMPLEX term = cfak * olk(0, 0);
266 for (
Index i = 1; i < n_orbitals[lmax_row]; i++) {
267 COMPLEX term_loc = cfak * olk(i, 0);
269 double(nx[i]) * cfak * olk(i_less_x[i],
Cart::x) +
272 double(nx[i]) * cfak * olk(i_less_x[i],
Cart::y);
274 double(nx[i]) * cfak * olk(i_less_x[i],
Cart::z);
276 double(ny[i]) * cfak * olk(i_less_y[i],
Cart::y) +
279 double(ny[i]) * cfak * olk(i_less_y[i],
Cart::z);
281 double(nz[i]) * cfak * olk(i_less_z[i],
Cart::z) +
307 for (
Index i = 1; i < n_orbitals[lmax_row]; i++) {
308 COMPLEX term_x = 2.0 * cfak * olk(i,
Cart::x);
309 COMPLEX term_y = 2.0 * cfak * olk(i,
Cart::y);
310 COMPLEX term_z = 2.0 * cfak * olk(i,
Cart::z);
313 double(nx[i]) * cfak * olk(i_less_x[i],
Cart::xx) + term_x;
315 double(ny[i]) * cfak * olk(i_less_y[i],
Cart::xx);
317 double(nz[i]) * cfak * olk(i_less_z[i],
Cart::xx);
319 double(nx[i]) * cfak * olk(i_less_x[i],
Cart::yy);
321 double(nx[i]) * cfak * olk(i_less_x[i],
Cart::yz);
323 double(nx[i]) * cfak * olk(i_less_x[i],
Cart::zz);
326 double(ny[i]) * cfak * olk(i_less_y[i],
Cart::yy) + term_y;
328 double(nz[i]) * cfak * olk(i_less_z[i],
Cart::yy);
330 double(ny[i]) * cfak * olk(i_less_y[i],
Cart::zz);
333 double(nz[i]) * cfak * olk(i_less_z[i],
Cart::zz) + term_z;
342 COMPLEX term_xx = cfak * olk(0,
Cart::xx);
343 COMPLEX term_yy = cfak * olk(0,
Cart::yy);
344 COMPLEX term_zz = cfak * olk(0,
Cart::zz);
363 for (
Index i = 1; i < n_orbitals[lmax_row]; i++) {
364 COMPLEX term_xx_loc = cfak * olk(i,
Cart::xx);
365 COMPLEX term_yy_loc = cfak * olk(i,
Cart::yy);
366 COMPLEX term_zz_loc = cfak * olk(i,
Cart::zz);
369 double(nx[i]) * cfak * olk(i_less_x[i],
Cart::xxx) +
373 double(ny[i]) * cfak * olk(i_less_y[i],
Cart::xxx);
376 double(nz[i]) * cfak * olk(i_less_z[i],
Cart::xxx);
379 double(nx[i]) * cfak * olk(i_less_x[i],
Cart::xyy) + term_yy_loc;
382 double(ny[i]) * cfak * olk(i_less_y[i],
Cart::xxz);
385 double(nx[i]) * cfak * olk(i_less_x[i],
Cart::xzz) + term_zz_loc;
388 double(nx[i]) * cfak * olk(i_less_x[i],
Cart::yyy);
391 double(nx[i]) * cfak * olk(i_less_x[i],
Cart::yyz);
394 double(nx[i]) * cfak * olk(i_less_x[i],
Cart::yzz);
397 double(nx[i]) * cfak * olk(i_less_x[i],
Cart::zzz);
400 double(ny[i]) * cfak * olk(i_less_y[i],
Cart::yyy) +
404 double(nz[i]) * cfak * olk(i_less_z[i],
Cart::yyy);
407 double(ny[i]) * cfak * olk(i_less_y[i],
Cart::yzz) + term_zz_loc;
410 double(ny[i]) * cfak * olk(i_less_y[i],
Cart::zzz);
413 double(nz[i]) * cfak * olk(i_less_z[i],
Cart::zzz) +
423 COMPLEX term_xxx = cfak * olk(0,
Cart::xxx);
424 COMPLEX term_yyy = cfak * olk(0,
Cart::yyy);
425 COMPLEX term_zzz = cfak * olk(0,
Cart::zzz);
450 for (
Index i = 1; i < n_orbitals[lmax_row]; i++) {
451 COMPLEX term_xxx_loc = cfak * olk(i,
Cart::xxx);
452 COMPLEX term_yyy_loc = cfak * olk(i,
Cart::yyy);
453 COMPLEX term_zzz_loc = cfak * olk(i,
Cart::zzz);
456 double(nx[i]) * cfak * olk(i_less_x[i],
Cart::xxxx) +
460 double(ny[i]) * cfak * olk(i_less_y[i],
Cart::xxxx);
463 double(nz[i]) * cfak * olk(i_less_z[i],
Cart::xxxx);
466 double(ny[i]) * cfak * olk(i_less_y[i],
Cart::xxxy) +
470 double(ny[i]) * cfak * olk(i_less_y[i],
Cart::xxxz);
473 double(nz[i]) * cfak * olk(i_less_z[i],
Cart::xxxz) +
477 double(nx[i]) * cfak * olk(i_less_x[i],
Cart::xyyy) +
481 double(nz[i]) * cfak * olk(i_less_z[i],
Cart::xxyy);
484 double(ny[i]) * cfak * olk(i_less_y[i],
Cart::xxzz);
487 double(nx[i]) * cfak * olk(i_less_x[i],
Cart::xzzz) +
491 double(nx[i]) * cfak * olk(i_less_x[i],
Cart::yyyy);
494 double(nx[i]) * cfak * olk(i_less_x[i],
Cart::yyyz);
497 double(nx[i]) * cfak * olk(i_less_x[i],
Cart::yyzz);
500 double(nx[i]) * cfak * olk(i_less_x[i],
Cart::yzzz);
503 double(nx[i]) * cfak * olk(i_less_x[i],
Cart::zzzz);
506 double(ny[i]) * cfak * olk(i_less_y[i],
Cart::yyyy) +
510 double(nz[i]) * cfak * olk(i_less_z[i],
Cart::yyyy);
513 double(nz[i]) * cfak * olk(i_less_z[i],
Cart::yyyz) +
517 double(ny[i]) * cfak * olk(i_less_y[i],
Cart::yzzz) +
521 double(ny[i]) * cfak * olk(i_less_y[i],
Cart::zzzz);
524 double(nz[i]) * cfak * olk(i_less_z[i],
Cart::zzzz) +
534 COMPLEX term_xxxx = cfak * olk(0,
Cart::xxxx);
535 COMPLEX term_xyyy = cfak * olk(0,
Cart::xyyy);
536 COMPLEX term_xzzz = cfak * olk(0,
Cart::xzzz);
537 COMPLEX term_yyyy = cfak * olk(0,
Cart::yyyy);
538 COMPLEX term_yyzz = cfak * olk(0,
Cart::yyzz);
539 COMPLEX term_yzzz = cfak * olk(0,
Cart::yzzz);
540 COMPLEX term_zzzz = cfak * olk(0,
Cart::zzzz);
573 for (
Index i = 1; i < n_orbitals[lmax_row]; i++) {
574 COMPLEX term_xxxx_loc = cfak * olk(i,
Cart::xxxx);
575 COMPLEX term_xyyy_loc = cfak * olk(i,
Cart::xyyy);
576 COMPLEX term_xzzz_loc = cfak * olk(i,
Cart::xzzz);
577 COMPLEX term_yyyy_loc = cfak * olk(i,
Cart::yyyy);
578 COMPLEX term_yyzz_loc = cfak * olk(i,
Cart::yyzz);
579 COMPLEX term_yzzz_loc = cfak * olk(i,
Cart::yzzz);
580 COMPLEX term_zzzz_loc = cfak * olk(i,
Cart::zzzz);
583 double(nx[i]) * cfak * olk(i_less_x[i],
Cart::xxxxx) +
587 double(ny[i]) * cfak * olk(i_less_y[i],
Cart::xxxxx);
590 double(nz[i]) * cfak * olk(i_less_z[i],
Cart::xxxxx);
593 double(ny[i]) * cfak * olk(i_less_y[i],
Cart::xxxxy) +
597 double(ny[i]) * cfak * olk(i_less_y[i],
Cart::xxxxz);
600 double(nz[i]) * cfak * olk(i_less_z[i],
Cart::xxxxz) +
604 double(nx[i]) * cfak * olk(i_less_x[i],
Cart::xxyyy) +
608 double(nz[i]) * cfak * olk(i_less_z[i],
Cart::xxxyy);
611 double(ny[i]) * cfak * olk(i_less_y[i],
Cart::xxxzz);
614 double(nx[i]) * cfak * olk(i_less_x[i],
Cart::xxzzz) +
618 double(nx[i]) * cfak * olk(i_less_x[i],
Cart::xyyyy) +
622 double(nz[i]) * cfak * olk(i_less_z[i],
Cart::xxyyy);
625 double(nx[i]) * cfak * olk(i_less_x[i],
Cart::xyyzz) +
629 double(ny[i]) * cfak * olk(i_less_y[i],
Cart::xxzzz);
632 double(nx[i]) * cfak * olk(i_less_x[i],
Cart::xzzzz) +
636 double(nx[i]) * cfak * olk(i_less_x[i],
Cart::yyyyy);
639 double(nx[i]) * cfak * olk(i_less_x[i],
Cart::yyyyz);
642 double(nx[i]) * cfak * olk(i_less_x[i],
Cart::yyyzz);
645 double(nx[i]) * cfak * olk(i_less_x[i],
Cart::yyzzz);
648 double(nx[i]) * cfak * olk(i_less_x[i],
Cart::yzzzz);
651 double(nx[i]) * cfak * olk(i_less_x[i],
Cart::zzzzz);
654 double(ny[i]) * cfak * olk(i_less_y[i],
Cart::yyyyy) +
658 double(nz[i]) * cfak * olk(i_less_z[i],
Cart::yyyyy);
661 double(nz[i]) * cfak * olk(i_less_z[i],
Cart::yyyyz) +
665 double(ny[i]) * cfak * olk(i_less_y[i],
Cart::yyzzz) +
669 double(ny[i]) * cfak * olk(i_less_y[i],
Cart::yzzzz) +
673 double(ny[i]) * cfak * olk(i_less_y[i],
Cart::zzzzz);
676 double(nz[i]) * cfak * olk(i_less_z[i],
Cart::zzzzz) +
687 shell_col.getCartesianNumFunc());