votca 2026-dev
Loading...
Searching...
No Matches
aomultipole.cc
Go to the documentation of this file.
1/*
2 * Copyright 2009-2020 The VOTCA Development Team
3 * (http://www.votca.org)
4 *
5 * Licensed under the Apache License, Version 2.0 (the "License")
6 *
7 * You may not use this file except in compliance with the License.
8 * You may obtain a copy of the License at
9 *
10 * http://www.apache.org/licenses/LICENSE-2.0
11 *
12 *Unless required by applicable law or agreed to in writing, software
13 * distributed under the License is distributed on an "AS IS" BASIS,
14 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
15 * See the License for the specific language governing permissions and
16 * limitations under the License.
17 *
18 */
19
20// Local VOTCA includes
21#include "votca/xtp/aopotential.h"
22#include "votca/xtp/aotransform.h"
23#include "votca/xtp/qmmolecule.h"
24
25namespace votca {
26namespace xtp {
27
28namespace {
29constexpr double kMaxExpArg = 50.0;
30
31// How a Cartesian function t is built in the recursions: along direction d
32// (its first nonzero exponent, as in the explicit recursions: xy from y
33// along x) from parent = t - 1_d, with n = exponent of d in the parent and
34// grand = parent - 1_d (if n > 0); l is the angular momentum of t.
35struct BuildStep {
36 int d;
37 int parent;
38 int n;
39 int grand;
40 int l;
41};
42
43const std::array<BuildStep, 165>& BuildSteps() {
44 static const std::array<BuildStep, 165> steps = [] {
45 const std::array<int, 165> nx = AOTransform::nx();
46 const std::array<int, 165> ny = AOTransform::ny();
47 const std::array<int, 165> nz = AOTransform::nz();
48 const std::array<std::array<int, 165>, 3> less = {AOTransform::i_less_x(),
49 AOTransform::i_less_y(),
50 AOTransform::i_less_z()};
51 const std::array<const std::array<int, 165>*, 3> n = {&nx, &ny, &nz};
52 std::array<BuildStep, 165> result{};
53 for (std::size_t t = 1; t < 165; ++t) {
54 BuildStep step;
55 step.d = (nx[t] > 0) ? 0 : ((ny[t] > 0) ? 1 : 2);
56 step.parent = less[std::size_t(step.d)][t];
57 step.n = (*n[std::size_t(step.d)])[std::size_t(step.parent)];
58 step.grand = (step.n > 0)
59 ? less[std::size_t(step.d)][std::size_t(step.parent)]
60 : 0;
61 step.l = nx[t] + ny[t] + nz[t];
62 result[t] = step;
63 }
64 return result;
65 }();
66 return steps;
67}
68
69// Potential integrals of a point dipole mu, contracted with mu during the
70// recursion: one Obara-Saika recursion for mu . (dipole integrals) instead
71// of one per Cartesian component. The recursion is linear in the dipole
72// component k, and k enters only through the start values (PmC_k) and the
73// term delta_kd * nuc(m+1), so contracting replaces those by PmC . mu and
74// mu_d * nuc(m+1). Same terms and same build order as the per-component
75// recursion; on return dipc(:, :, 0) holds the contracted integrals.
76void ContractedDipoleIntegrals(const Eigen::Tensor<double, 3>& nuc,
77 Eigen::Tensor<double, 3>& dipc, Index lmax_row,
78 Index lmax_col, const Eigen::Vector3d& PmA,
79 const Eigen::Vector3d& PmB,
80 const Eigen::Vector3d& PmC, double fak,
81 double prefactor_dip, const Eigen::VectorXd& FmU,
82 const Eigen::Vector3d& mu) {
83 // the index tables are fixed: built once, not copied per call
84 static const std::array<int, 9> n_orbitals = AOTransform::n_orbitals();
85 static const std::array<std::array<int, 165>, 3> n_of = {
86 AOTransform::nx(), AOTransform::ny(), AOTransform::nz()};
87 static const std::array<std::array<int, 165>, 3> less = {
88 AOTransform::i_less_x(), AOTransform::i_less_y(),
89 AOTransform::i_less_z()};
90 const std::array<BuildStep, 165>& steps = BuildSteps();
91
92 const Index lsum = lmax_row + lmax_col;
93 const Index nrows = dipc.dimension(0);
94 const Index s_col = nrows; // stride of the column index
95 const Index s_m = nrows * dipc.dimension(1); // stride of m
96 double* D = dipc.data();
97 const double* N = nuc.data();
98 const Index rows = n_orbitals[lmax_row];
99
100 const double start = PmC.dot(mu) * prefactor_dip;
101 for (Index m = 0; m < lsum + 1; m++) {
102 D[m * s_m] = start * FmU[m + 1];
103 }
104 // rows: (t, s) from (t - 1_d, s)
105 for (Index t = 1; t < rows; t++) {
106 const BuildStep& st = steps[std::size_t(t)];
107 const double pa = PmA(st.d);
108 const double pc = PmC(st.d);
109 const double md = mu(st.d);
110 const double nf = st.n * fak;
111 for (Index m = 0; m < lsum - st.l + 1; m++) {
112 const Index o = m * s_m;
113 D[t + o] = pa * D[st.parent + o] - pc * D[st.parent + o + s_m] +
114 md * N[st.parent + o + s_m] +
115 nf * (D[st.grand + o] - D[st.grand + o + s_m]);
116 }
117 }
118 // columns: (i, c) from (i, c - 1_d)
119 for (Index c = 1; c < n_orbitals[lmax_col]; c++) {
120 const BuildStep& st = steps[std::size_t(c)];
121 const double pb = PmB(st.d);
122 const double pc = PmC(st.d);
123 const double md = mu(st.d);
124 const double nqf = st.n * fak;
125 const std::array<int, 165>& n_i = n_of[std::size_t(st.d)];
126 const std::array<int, 165>& less_i = less[std::size_t(st.d)];
127 for (Index m = 0; m < lmax_col - st.l + 1; m++) {
128 const Index o = m * s_m;
129 const double* Dq = D + st.parent * s_col + o;
130 const double* Dqq = D + st.grand * s_col + o;
131 const double* Nq = N + st.parent * s_col + o;
132 double* Dc = D + c * s_col + o;
133 for (Index i = 0; i < rows; i++) {
134 double value = pb * Dq[i] - pc * Dq[i + s_m] + md * Nq[i + s_m] +
135 nqf * (Dqq[i] - Dqq[i + s_m]);
136 const int ni = n_i[std::size_t(i)];
137 if (ni > 0) {
138 const int il = less_i[std::size_t(i)];
139 value += ni * fak * (Dq[il] - Dq[il + s_m]);
140 }
141 Dc[i] = value;
142 }
143 }
144 }
145}
146} // namespace
147
148void AOMultipole::FillBlock(Eigen::Block<Eigen::MatrixXd>& matrix,
149 const AOShell& shell_row,
150 const AOShell& shell_col) const {
151
152 const double pi = boost::math::constants::pi<double>();
153
154 // shell info, only lmax tells how far to go
155 Index lmax_row = Index(shell_row.getL());
156 Index lmax_col = Index(shell_col.getL());
157 Index lsum = lmax_row + lmax_col;
158 // set size of internal block for recursion
159 Index nrows = AOTransform::getBlockSize(lmax_row);
160 Index ncols = AOTransform::getBlockSize(lmax_col);
161
162 std::array<int, 9> n_orbitals = AOTransform::n_orbitals();
163 std::array<int, 165> nx = AOTransform::nx();
164 std::array<int, 165> ny = AOTransform::ny();
165 std::array<int, 165> nz = AOTransform::nz();
166 std::array<int, 165> i_less_x = AOTransform::i_less_x();
167 std::array<int, 165> i_less_y = AOTransform::i_less_y();
168 std::array<int, 165> i_less_z = AOTransform::i_less_z();
169
170 // get shell positions
171 const Eigen::Vector3d& pos_row = shell_row.getPos();
172 const Eigen::Vector3d& pos_col = shell_col.getPos();
173 const Eigen::Vector3d diff = pos_row - pos_col;
174 // initialize some helper
175
176 double distsq = diff.squaredNorm();
177
178 Eigen::MatrixXd cartesian = Eigen::MatrixXd::Zero(
179 shell_row.getCartesianNumFunc(), shell_col.getCartesianNumFunc());
180
181 // recursion buffers; every entry a recursion reads is written first, so
182 // they are not zeroed per site (checked by filling them with NaN)
183 Eigen::Tensor<double, 3> nuc3(nrows, ncols, lsum + 1);
184 Eigen::Tensor<double, 4> dip4(nrows, ncols, 3, lsum + 1);
185 Eigen::Tensor<double, 3> dipc(nrows, ncols, lsum + 1);
186 Eigen::Tensor<double, 4> quad4(nrows, ncols, 5, lsum + 1);
187 Eigen::MatrixXd multipole(nrows, ncols);
188 Eigen::VectorXd FmU(lsum + 3); // Boys function values, rank <= 2
189
190 // iterate over Gaussians in this shell_row
191 for (const auto& gaussian_row : shell_row) {
192 // iterate over Gaussians in this shell_col
193 // get decay constant
194 const double decay_row = gaussian_row.getDecay();
195
196 for (const auto& gaussian_col : shell_col) {
197 // get decay constant
198 const double decay_col = gaussian_col.getDecay();
199
200 const double zeta = decay_row + decay_col;
201 const double fak = 0.5 / zeta;
202 const double fak2 = 2.0 * fak;
203 const double xi = decay_row * decay_col * fak2;
204
205 double exparg = xi * distsq;
206 // The whole primitive pair carries exp(-exparg); beyond this it is
207 // below 2e-22 for every site and cannot change a matrix element.
208 if (exparg > kMaxExpArg) {
209 continue;
210 }
211
212 // some helpers
213 const Eigen::Vector3d PmA =
214 fak2 * (decay_row * pos_row + decay_col * pos_col) - pos_row;
215 const Eigen::Vector3d PmB =
216 fak2 * (decay_row * pos_row + decay_col * pos_col) - pos_col;
217
218 // summed over all sites, normalised once per primitive pair
219 multipole.setZero();
220 for (const SiteData& site : sites_) {
221 const Index rank = site.rank;
222 const double charge = site.charge;
223 const Eigen::Vector3d& dipole = site.dipole;
224 const Eigen::Matrix3d& quadrupole = site.quadrupole;
225 const Eigen::Vector3d PmC =
226 fak2 * (decay_row * pos_row + decay_col * pos_col) - site.pos;
227
228 const double U = zeta * PmC.squaredNorm();
229
230 // +3 quadrupole, +2 dipole, +1 nuclear attraction integrals
231 AOTransform::XIntegrate(lsum + rank + 1, U, FmU.data());
232
233 // (s-s element normiert )
234 double prefactor = 4. * sqrt(2. / pi) *
235 pow(decay_row * decay_col, .75) * fak2 *
236 exp(-exparg);
237 for (Index m = 0; m < lsum + 1; m++) {
238 nuc3(0, 0, m) = prefactor * FmU[m];
239 }
240 //------------------------------------------------------
241
242 // Integrals p - s
243 if (lmax_row > 0) {
244 for (Index m = 0; m < lsum; m++) {
245 nuc3(Cart::x, 0, m) =
246 PmA(0) * nuc3(0, 0, m) - PmC(0) * nuc3(0, 0, m + 1);
247 nuc3(Cart::y, 0, m) =
248 PmA(1) * nuc3(0, 0, m) - PmC(1) * nuc3(0, 0, m + 1);
249 nuc3(Cart::z, 0, m) =
250 PmA(2) * nuc3(0, 0, m) - PmC(2) * nuc3(0, 0, m + 1);
251 }
252 }
253 //------------------------------------------------------
254
255 // Integrals d - s
256 if (lmax_row > 1) {
257 for (Index m = 0; m < lsum - 1; m++) {
258 double term = fak * (nuc3(0, 0, m) - nuc3(0, 0, m + 1));
259 nuc3(Cart::xx, 0, m) = PmA(0) * nuc3(Cart::x, 0, m) -
260 PmC(0) * nuc3(Cart::x, 0, m + 1) + term;
261 nuc3(Cart::xy, 0, m) =
262 PmA(0) * nuc3(Cart::y, 0, m) - PmC(0) * nuc3(Cart::y, 0, m + 1);
263 nuc3(Cart::xz, 0, m) =
264 PmA(0) * nuc3(Cart::z, 0, m) - PmC(0) * nuc3(Cart::z, 0, m + 1);
265 nuc3(Cart::yy, 0, m) = PmA(1) * nuc3(Cart::y, 0, m) -
266 PmC(1) * nuc3(Cart::y, 0, m + 1) + term;
267 nuc3(Cart::yz, 0, m) =
268 PmA(1) * nuc3(Cart::z, 0, m) - PmC(1) * nuc3(Cart::z, 0, m + 1);
269 nuc3(Cart::zz, 0, m) = PmA(2) * nuc3(Cart::z, 0, m) -
270 PmC(2) * nuc3(Cart::z, 0, m + 1) + term;
271 }
272 }
273 //------------------------------------------------------
274
275 // Integrals f - s
276 if (lmax_row > 2) {
277 for (Index m = 0; m < lsum - 2; m++) {
278 nuc3(Cart::xxx, 0, m) =
279 PmA(0) * nuc3(Cart::xx, 0, m) -
280 PmC(0) * nuc3(Cart::xx, 0, m + 1) +
281 2 * fak * (nuc3(Cart::x, 0, m) - nuc3(Cart::x, 0, m + 1));
282 nuc3(Cart::xxy, 0, m) = PmA(1) * nuc3(Cart::xx, 0, m) -
283 PmC(1) * nuc3(Cart::xx, 0, m + 1);
284 nuc3(Cart::xxz, 0, m) = PmA(2) * nuc3(Cart::xx, 0, m) -
285 PmC(2) * nuc3(Cart::xx, 0, m + 1);
286 nuc3(Cart::xyy, 0, m) = PmA(0) * nuc3(Cart::yy, 0, m) -
287 PmC(0) * nuc3(Cart::yy, 0, m + 1);
288 nuc3(Cart::xyz, 0, m) = PmA(0) * nuc3(Cart::yz, 0, m) -
289 PmC(0) * nuc3(Cart::yz, 0, m + 1);
290 nuc3(Cart::xzz, 0, m) = PmA(0) * nuc3(Cart::zz, 0, m) -
291 PmC(0) * nuc3(Cart::zz, 0, m + 1);
292 nuc3(Cart::yyy, 0, m) =
293 PmA(1) * nuc3(Cart::yy, 0, m) -
294 PmC(1) * nuc3(Cart::yy, 0, m + 1) +
295 2 * fak * (nuc3(Cart::y, 0, m) - nuc3(Cart::y, 0, m + 1));
296 nuc3(Cart::yyz, 0, m) = PmA(2) * nuc3(Cart::yy, 0, m) -
297 PmC(2) * nuc3(Cart::yy, 0, m + 1);
298 nuc3(Cart::yzz, 0, m) = PmA(1) * nuc3(Cart::zz, 0, m) -
299 PmC(1) * nuc3(Cart::zz, 0, m + 1);
300 nuc3(Cart::zzz, 0, m) =
301 PmA(2) * nuc3(Cart::zz, 0, m) -
302 PmC(2) * nuc3(Cart::zz, 0, m + 1) +
303 2 * fak * (nuc3(Cart::z, 0, m) - nuc3(Cart::z, 0, m + 1));
304 }
305 }
306 //------------------------------------------------------
307
308 // Integrals g - s
309 if (lmax_row > 3) {
310 for (Index m = 0; m < lsum - 3; m++) {
311 double term_xx =
312 fak * (nuc3(Cart::xx, 0, m) - nuc3(Cart::xx, 0, m + 1));
313 double term_yy =
314 fak * (nuc3(Cart::yy, 0, m) - nuc3(Cart::yy, 0, m + 1));
315 double term_zz =
316 fak * (nuc3(Cart::zz, 0, m) - nuc3(Cart::zz, 0, m + 1));
317 nuc3(Cart::xxxx, 0, m) = PmA(0) * nuc3(Cart::xxx, 0, m) -
318 PmC(0) * nuc3(Cart::xxx, 0, m + 1) +
319 3 * term_xx;
320 nuc3(Cart::xxxy, 0, m) = PmA(1) * nuc3(Cart::xxx, 0, m) -
321 PmC(1) * nuc3(Cart::xxx, 0, m + 1);
322 nuc3(Cart::xxxz, 0, m) = PmA(2) * nuc3(Cart::xxx, 0, m) -
323 PmC(2) * nuc3(Cart::xxx, 0, m + 1);
324 nuc3(Cart::xxyy, 0, m) = PmA(0) * nuc3(Cart::xyy, 0, m) -
325 PmC(0) * nuc3(Cart::xyy, 0, m + 1) +
326 term_yy;
327 nuc3(Cart::xxyz, 0, m) = PmA(1) * nuc3(Cart::xxz, 0, m) -
328 PmC(1) * nuc3(Cart::xxz, 0, m + 1);
329 nuc3(Cart::xxzz, 0, m) = PmA(0) * nuc3(Cart::xzz, 0, m) -
330 PmC(0) * nuc3(Cart::xzz, 0, m + 1) +
331 term_zz;
332 nuc3(Cart::xyyy, 0, m) = PmA(0) * nuc3(Cart::yyy, 0, m) -
333 PmC(0) * nuc3(Cart::yyy, 0, m + 1);
334 nuc3(Cart::xyyz, 0, m) = PmA(0) * nuc3(Cart::yyz, 0, m) -
335 PmC(0) * nuc3(Cart::yyz, 0, m + 1);
336 nuc3(Cart::xyzz, 0, m) = PmA(0) * nuc3(Cart::yzz, 0, m) -
337 PmC(0) * nuc3(Cart::yzz, 0, m + 1);
338 nuc3(Cart::xzzz, 0, m) = PmA(0) * nuc3(Cart::zzz, 0, m) -
339 PmC(0) * nuc3(Cart::zzz, 0, m + 1);
340 nuc3(Cart::yyyy, 0, m) = PmA(1) * nuc3(Cart::yyy, 0, m) -
341 PmC(1) * nuc3(Cart::yyy, 0, m + 1) +
342 3 * term_yy;
343 nuc3(Cart::yyyz, 0, m) = PmA(2) * nuc3(Cart::yyy, 0, m) -
344 PmC(2) * nuc3(Cart::yyy, 0, m + 1);
345 nuc3(Cart::yyzz, 0, m) = PmA(1) * nuc3(Cart::yzz, 0, m) -
346 PmC(1) * nuc3(Cart::yzz, 0, m + 1) +
347 term_zz;
348 nuc3(Cart::yzzz, 0, m) = PmA(1) * nuc3(Cart::zzz, 0, m) -
349 PmC(1) * nuc3(Cart::zzz, 0, m + 1);
350 nuc3(Cart::zzzz, 0, m) = PmA(2) * nuc3(Cart::zzz, 0, m) -
351 PmC(2) * nuc3(Cart::zzz, 0, m + 1) +
352 3 * term_zz;
353 }
354 }
355 //------------------------------------------------------
356
357 if (lmax_col > 0) {
358
359 // Integrals s - p
360 for (Index m = 0; m < lmax_col; m++) {
361 nuc3(0, Cart::x, m) =
362 PmB(0) * nuc3(0, 0, m) - PmC(0) * nuc3(0, 0, m + 1);
363 nuc3(0, Cart::y, m) =
364 PmB(1) * nuc3(0, 0, m) - PmC(1) * nuc3(0, 0, m + 1);
365 nuc3(0, Cart::z, m) =
366 PmB(2) * nuc3(0, 0, m) - PmC(2) * nuc3(0, 0, m + 1);
367 }
368 //------------------------------------------------------
369
370 // Integrals p - p
371 if (lmax_row > 0) {
372 for (Index m = 0; m < lmax_col; m++) {
373 double term = fak * (nuc3(0, 0, m) - nuc3(0, 0, m + 1));
374 for (Index i = 1; i < 4; i++) {
375 nuc3(i, Cart::x, m) = PmB(0) * nuc3(i, 0, m) -
376 PmC(0) * nuc3(i, 0, m + 1) + nx[i] * term;
377 nuc3(i, Cart::y, m) = PmB(1) * nuc3(i, 0, m) -
378 PmC(1) * nuc3(i, 0, m + 1) + ny[i] * term;
379 nuc3(i, Cart::z, m) = PmB(2) * nuc3(i, 0, m) -
380 PmC(2) * nuc3(i, 0, m + 1) + nz[i] * term;
381 }
382 }
383 }
384 //------------------------------------------------------
385
386 // Integrals d - p f - p g - p
387 for (Index m = 0; m < lmax_col; m++) {
388 for (Index i = 4; i < n_orbitals[lmax_row]; i++) {
389 int nx_i = nx[i];
390 int ny_i = ny[i];
391 int nz_i = nz[i];
392 int ilx_i = i_less_x[i];
393 int ily_i = i_less_y[i];
394 int ilz_i = i_less_z[i];
395 nuc3(i, Cart::x, m) =
396 PmB(0) * nuc3(i, 0, m) - PmC(0) * nuc3(i, 0, m + 1) +
397 nx_i * fak * (nuc3(ilx_i, 0, m) - nuc3(ilx_i, 0, m + 1));
398 nuc3(i, Cart::y, m) =
399 PmB(1) * nuc3(i, 0, m) - PmC(1) * nuc3(i, 0, m + 1) +
400 ny_i * fak * (nuc3(ily_i, 0, m) - nuc3(ily_i, 0, m + 1));
401 nuc3(i, Cart::z, m) =
402 PmB(2) * nuc3(i, 0, m) - PmC(2) * nuc3(i, 0, m + 1) +
403 nz_i * fak * (nuc3(ilz_i, 0, m) - nuc3(ilz_i, 0, m + 1));
404 }
405 }
406 //------------------------------------------------------
407
408 } // end if (lmax_col > 0)
409
410 if (lmax_col > 1) {
411
412 // Integrals s - d
413 for (Index m = 0; m < lmax_col - 1; m++) {
414 double term = fak * (nuc3(0, 0, m) - nuc3(0, 0, m + 1));
415 nuc3(0, Cart::xx, m) = PmB(0) * nuc3(0, Cart::x, m) -
416 PmC(0) * nuc3(0, Cart::x, m + 1) + term;
417 nuc3(0, Cart::xy, m) =
418 PmB(0) * nuc3(0, Cart::y, m) - PmC(0) * nuc3(0, Cart::y, m + 1);
419 nuc3(0, Cart::xz, m) =
420 PmB(0) * nuc3(0, Cart::z, m) - PmC(0) * nuc3(0, Cart::z, m + 1);
421 nuc3(0, Cart::yy, m) = PmB(1) * nuc3(0, Cart::y, m) -
422 PmC(1) * nuc3(0, Cart::y, m + 1) + term;
423 nuc3(0, Cart::yz, m) =
424 PmB(1) * nuc3(0, Cart::z, m) - PmC(1) * nuc3(0, Cart::z, m + 1);
425 nuc3(0, Cart::zz, m) = PmB(2) * nuc3(0, Cart::z, m) -
426 PmC(2) * nuc3(0, Cart::z, m + 1) + term;
427 }
428 //------------------------------------------------------
429
430 // Integrals p - d d - d f - d g - d
431 for (Index m = 0; m < lmax_col - 1; m++) {
432 for (Index i = 1; i < n_orbitals[lmax_row]; i++) {
433 int nx_i = nx[i];
434 int ny_i = ny[i];
435 int nz_i = nz[i];
436 int ilx_i = i_less_x[i];
437 int ily_i = i_less_y[i];
438 int ilz_i = i_less_z[i];
439 double term = fak * (nuc3(i, 0, m) - nuc3(i, 0, m + 1));
440 nuc3(i, Cart::xx, m) =
441 PmB(0) * nuc3(i, Cart::x, m) -
442 PmC(0) * nuc3(i, Cart::x, m + 1) +
443 nx_i * fak *
444 (nuc3(ilx_i, Cart::x, m) - nuc3(ilx_i, Cart::x, m + 1)) +
445 term;
446 nuc3(i, Cart::xy, m) =
447 PmB(0) * nuc3(i, Cart::y, m) -
448 PmC(0) * nuc3(i, Cart::y, m + 1) +
449 nx_i * fak *
450 (nuc3(ilx_i, Cart::y, m) - nuc3(ilx_i, Cart::y, m + 1));
451 nuc3(i, Cart::xz, m) =
452 PmB(0) * nuc3(i, Cart::z, m) -
453 PmC(0) * nuc3(i, Cart::z, m + 1) +
454 nx_i * fak *
455 (nuc3(ilx_i, Cart::z, m) - nuc3(ilx_i, Cart::z, m + 1));
456 nuc3(i, Cart::yy, m) =
457 PmB(1) * nuc3(i, Cart::y, m) -
458 PmC(1) * nuc3(i, Cart::y, m + 1) +
459 ny_i * fak *
460 (nuc3(ily_i, Cart::y, m) - nuc3(ily_i, Cart::y, m + 1)) +
461 term;
462 nuc3(i, Cart::yz, m) =
463 PmB(1) * nuc3(i, Cart::z, m) -
464 PmC(1) * nuc3(i, Cart::z, m + 1) +
465 ny_i * fak *
466 (nuc3(ily_i, Cart::z, m) - nuc3(ily_i, Cart::z, m + 1));
467 nuc3(i, Cart::zz, m) =
468 PmB(2) * nuc3(i, Cart::z, m) -
469 PmC(2) * nuc3(i, Cart::z, m + 1) +
470 nz_i * fak *
471 (nuc3(ilz_i, Cart::z, m) - nuc3(ilz_i, Cart::z, m + 1)) +
472 term;
473 }
474 }
475 //------------------------------------------------------
476
477 } // end if (lmax_col > 1)
478
479 if (lmax_col > 2) {
480
481 // Integrals s - f
482 for (Index m = 0; m < lmax_col - 2; m++) {
483 nuc3(0, Cart::xxx, m) =
484 PmB(0) * nuc3(0, Cart::xx, m) -
485 PmC(0) * nuc3(0, Cart::xx, m + 1) +
486 2 * fak * (nuc3(0, Cart::x, m) - nuc3(0, Cart::x, m + 1));
487 nuc3(0, Cart::xxy, m) = PmB(1) * nuc3(0, Cart::xx, m) -
488 PmC(1) * nuc3(0, Cart::xx, m + 1);
489 nuc3(0, Cart::xxz, m) = PmB(2) * nuc3(0, Cart::xx, m) -
490 PmC(2) * nuc3(0, Cart::xx, m + 1);
491 nuc3(0, Cart::xyy, m) = PmB(0) * nuc3(0, Cart::yy, m) -
492 PmC(0) * nuc3(0, Cart::yy, m + 1);
493 nuc3(0, Cart::xyz, m) = PmB(0) * nuc3(0, Cart::yz, m) -
494 PmC(0) * nuc3(0, Cart::yz, m + 1);
495 nuc3(0, Cart::xzz, m) = PmB(0) * nuc3(0, Cart::zz, m) -
496 PmC(0) * nuc3(0, Cart::zz, m + 1);
497 nuc3(0, Cart::yyy, m) =
498 PmB(1) * nuc3(0, Cart::yy, m) -
499 PmC(1) * nuc3(0, Cart::yy, m + 1) +
500 2 * fak * (nuc3(0, Cart::y, m) - nuc3(0, Cart::y, m + 1));
501 nuc3(0, Cart::yyz, m) = PmB(2) * nuc3(0, Cart::yy, m) -
502 PmC(2) * nuc3(0, Cart::yy, m + 1);
503 nuc3(0, Cart::yzz, m) = PmB(1) * nuc3(0, Cart::zz, m) -
504 PmC(1) * nuc3(0, Cart::zz, m + 1);
505 nuc3(0, Cart::zzz, m) =
506 PmB(2) * nuc3(0, Cart::zz, m) -
507 PmC(2) * nuc3(0, Cart::zz, m + 1) +
508 2 * fak * (nuc3(0, Cart::z, m) - nuc3(0, Cart::z, m + 1));
509 }
510 //------------------------------------------------------
511
512 // Integrals p - f d - f f - f g - f
513 for (Index m = 0; m < lmax_col - 2; m++) {
514 for (Index i = 1; i < n_orbitals[lmax_row]; i++) {
515 int nx_i = nx[i];
516 int ny_i = ny[i];
517 int nz_i = nz[i];
518 int ilx_i = i_less_x[i];
519 int ily_i = i_less_y[i];
520 int ilz_i = i_less_z[i];
521 double term_x =
522 2 * fak * (nuc3(i, Cart::x, m) - nuc3(i, Cart::x, m + 1));
523 double term_y =
524 2 * fak * (nuc3(i, Cart::y, m) - nuc3(i, Cart::y, m + 1));
525 double term_z =
526 2 * fak * (nuc3(i, Cart::z, m) - nuc3(i, Cart::z, m + 1));
527 nuc3(i, Cart::xxx, m) = PmB(0) * nuc3(i, Cart::xx, m) -
528 PmC(0) * nuc3(i, Cart::xx, m + 1) +
529 nx_i * fak *
530 (nuc3(ilx_i, Cart::xx, m) -
531 nuc3(ilx_i, Cart::xx, m + 1)) +
532 term_x;
533 nuc3(i, Cart::xxy, m) =
534 PmB(1) * nuc3(i, Cart::xx, m) -
535 PmC(1) * nuc3(i, Cart::xx, m + 1) +
536 ny_i * fak *
537 (nuc3(ily_i, Cart::xx, m) - nuc3(ily_i, Cart::xx, m + 1));
538 nuc3(i, Cart::xxz, m) =
539 PmB(2) * nuc3(i, Cart::xx, m) -
540 PmC(2) * nuc3(i, Cart::xx, m + 1) +
541 nz_i * fak *
542 (nuc3(ilz_i, Cart::xx, m) - nuc3(ilz_i, Cart::xx, m + 1));
543 nuc3(i, Cart::xyy, m) =
544 PmB(0) * nuc3(i, Cart::yy, m) -
545 PmC(0) * nuc3(i, Cart::yy, m + 1) +
546 nx_i * fak *
547 (nuc3(ilx_i, Cart::yy, m) - nuc3(ilx_i, Cart::yy, m + 1));
548 nuc3(i, Cart::xyz, m) =
549 PmB(0) * nuc3(i, Cart::yz, m) -
550 PmC(0) * nuc3(i, Cart::yz, m + 1) +
551 nx_i * fak *
552 (nuc3(ilx_i, Cart::yz, m) - nuc3(ilx_i, Cart::yz, m + 1));
553 nuc3(i, Cart::xzz, m) =
554 PmB(0) * nuc3(i, Cart::zz, m) -
555 PmC(0) * nuc3(i, Cart::zz, m + 1) +
556 nx_i * fak *
557 (nuc3(ilx_i, Cart::zz, m) - nuc3(ilx_i, Cart::zz, m + 1));
558 nuc3(i, Cart::yyy, m) = PmB(1) * nuc3(i, Cart::yy, m) -
559 PmC(1) * nuc3(i, Cart::yy, m + 1) +
560 ny_i * fak *
561 (nuc3(ily_i, Cart::yy, m) -
562 nuc3(ily_i, Cart::yy, m + 1)) +
563 term_y;
564 nuc3(i, Cart::yyz, m) =
565 PmB(2) * nuc3(i, Cart::yy, m) -
566 PmC(2) * nuc3(i, Cart::yy, m + 1) +
567 nz_i * fak *
568 (nuc3(ilz_i, Cart::yy, m) - nuc3(ilz_i, Cart::yy, m + 1));
569 nuc3(i, Cart::yzz, m) =
570 PmB(1) * nuc3(i, Cart::zz, m) -
571 PmC(1) * nuc3(i, Cart::zz, m + 1) +
572 ny_i * fak *
573 (nuc3(ily_i, Cart::zz, m) - nuc3(ily_i, Cart::zz, m + 1));
574 nuc3(i, Cart::zzz, m) = PmB(2) * nuc3(i, Cart::zz, m) -
575 PmC(2) * nuc3(i, Cart::zz, m + 1) +
576 nz_i * fak *
577 (nuc3(ilz_i, Cart::zz, m) -
578 nuc3(ilz_i, Cart::zz, m + 1)) +
579 term_z;
580 }
581 }
582 //------------------------------------------------------
583
584 } // end if (lmax_col > 2)
585
586 if (lmax_col > 3) {
587
588 // Integrals s - g
589 for (Index m = 0; m < lmax_col - 3; m++) {
590 double term_xx =
591 fak * (nuc3(0, Cart::xx, m) - nuc3(0, Cart::xx, m + 1));
592 double term_yy =
593 fak * (nuc3(0, Cart::yy, m) - nuc3(0, Cart::yy, m + 1));
594 double term_zz =
595 fak * (nuc3(0, Cart::zz, m) - nuc3(0, Cart::zz, m + 1));
596 nuc3(0, Cart::xxxx, m) = PmB(0) * nuc3(0, Cart::xxx, m) -
597 PmC(0) * nuc3(0, Cart::xxx, m + 1) +
598 3 * term_xx;
599 nuc3(0, Cart::xxxy, m) = PmB(1) * nuc3(0, Cart::xxx, m) -
600 PmC(1) * nuc3(0, Cart::xxx, m + 1);
601 nuc3(0, Cart::xxxz, m) = PmB(2) * nuc3(0, Cart::xxx, m) -
602 PmC(2) * nuc3(0, Cart::xxx, m + 1);
603 nuc3(0, Cart::xxyy, m) = PmB(0) * nuc3(0, Cart::xyy, m) -
604 PmC(0) * nuc3(0, Cart::xyy, m + 1) +
605 term_yy;
606 nuc3(0, Cart::xxyz, m) = PmB(1) * nuc3(0, Cart::xxz, m) -
607 PmC(1) * nuc3(0, Cart::xxz, m + 1);
608 nuc3(0, Cart::xxzz, m) = PmB(0) * nuc3(0, Cart::xzz, m) -
609 PmC(0) * nuc3(0, Cart::xzz, m + 1) +
610 term_zz;
611 nuc3(0, Cart::xyyy, m) = PmB(0) * nuc3(0, Cart::yyy, m) -
612 PmC(0) * nuc3(0, Cart::yyy, m + 1);
613 nuc3(0, Cart::xyyz, m) = PmB(0) * nuc3(0, Cart::yyz, m) -
614 PmC(0) * nuc3(0, Cart::yyz, m + 1);
615 nuc3(0, Cart::xyzz, m) = PmB(0) * nuc3(0, Cart::yzz, m) -
616 PmC(0) * nuc3(0, Cart::yzz, m + 1);
617 nuc3(0, Cart::xzzz, m) = PmB(0) * nuc3(0, Cart::zzz, m) -
618 PmC(0) * nuc3(0, Cart::zzz, m + 1);
619 nuc3(0, Cart::yyyy, m) = PmB(1) * nuc3(0, Cart::yyy, m) -
620 PmC(1) * nuc3(0, Cart::yyy, m + 1) +
621 3 * term_yy;
622 nuc3(0, Cart::yyyz, m) = PmB(2) * nuc3(0, Cart::yyy, m) -
623 PmC(2) * nuc3(0, Cart::yyy, m + 1);
624 nuc3(0, Cart::yyzz, m) = PmB(1) * nuc3(0, Cart::yzz, m) -
625 PmC(1) * nuc3(0, Cart::yzz, m + 1) +
626 term_zz;
627 nuc3(0, Cart::yzzz, m) = PmB(1) * nuc3(0, Cart::zzz, m) -
628 PmC(1) * nuc3(0, Cart::zzz, m + 1);
629 nuc3(0, Cart::zzzz, m) = PmB(2) * nuc3(0, Cart::zzz, m) -
630 PmC(2) * nuc3(0, Cart::zzz, m + 1) +
631 3 * term_zz;
632 }
633 //------------------------------------------------------
634
635 // Integrals p - g d - g f - g g - g
636 for (Index m = 0; m < lmax_col - 3; m++) {
637 for (Index i = 1; i < n_orbitals[lmax_row]; i++) {
638 int nx_i = nx[i];
639 int ny_i = ny[i];
640 int nz_i = nz[i];
641 int ilx_i = i_less_x[i];
642 int ily_i = i_less_y[i];
643 int ilz_i = i_less_z[i];
644 double term_xx =
645 fak * (nuc3(i, Cart::xx, m) - nuc3(i, Cart::xx, m + 1));
646 double term_yy =
647 fak * (nuc3(i, Cart::yy, m) - nuc3(i, Cart::yy, m + 1));
648 double term_zz =
649 fak * (nuc3(i, Cart::zz, m) - nuc3(i, Cart::zz, m + 1));
650 nuc3(i, Cart::xxxx, m) = PmB(0) * nuc3(i, Cart::xxx, m) -
651 PmC(0) * nuc3(i, Cart::xxx, m + 1) +
652 nx_i * fak *
653 (nuc3(ilx_i, Cart::xxx, m) -
654 nuc3(ilx_i, Cart::xxx, m + 1)) +
655 3 * term_xx;
656 nuc3(i, Cart::xxxy, m) = PmB(1) * nuc3(i, Cart::xxx, m) -
657 PmC(1) * nuc3(i, Cart::xxx, m + 1) +
658 ny_i * fak *
659 (nuc3(ily_i, Cart::xxx, m) -
660 nuc3(ily_i, Cart::xxx, m + 1));
661 nuc3(i, Cart::xxxz, m) = PmB(2) * nuc3(i, Cart::xxx, m) -
662 PmC(2) * nuc3(i, Cart::xxx, m + 1) +
663 nz_i * fak *
664 (nuc3(ilz_i, Cart::xxx, m) -
665 nuc3(ilz_i, Cart::xxx, m + 1));
666 nuc3(i, Cart::xxyy, m) = PmB(0) * nuc3(i, Cart::xyy, m) -
667 PmC(0) * nuc3(i, Cart::xyy, m + 1) +
668 nx_i * fak *
669 (nuc3(ilx_i, Cart::xyy, m) -
670 nuc3(ilx_i, Cart::xyy, m + 1)) +
671 term_yy;
672 nuc3(i, Cart::xxyz, m) = PmB(1) * nuc3(i, Cart::xxz, m) -
673 PmC(1) * nuc3(i, Cart::xxz, m + 1) +
674 ny_i * fak *
675 (nuc3(ily_i, Cart::xxz, m) -
676 nuc3(ily_i, Cart::xxz, m + 1));
677 nuc3(i, Cart::xxzz, m) = PmB(0) * nuc3(i, Cart::xzz, m) -
678 PmC(0) * nuc3(i, Cart::xzz, m + 1) +
679 nx_i * fak *
680 (nuc3(ilx_i, Cart::xzz, m) -
681 nuc3(ilx_i, Cart::xzz, m + 1)) +
682 term_zz;
683 nuc3(i, Cart::xyyy, m) = PmB(0) * nuc3(i, Cart::yyy, m) -
684 PmC(0) * nuc3(i, Cart::yyy, m + 1) +
685 nx_i * fak *
686 (nuc3(ilx_i, Cart::yyy, m) -
687 nuc3(ilx_i, Cart::yyy, m + 1));
688 nuc3(i, Cart::xyyz, m) = PmB(0) * nuc3(i, Cart::yyz, m) -
689 PmC(0) * nuc3(i, Cart::yyz, m + 1) +
690 nx_i * fak *
691 (nuc3(ilx_i, Cart::yyz, m) -
692 nuc3(ilx_i, Cart::yyz, m + 1));
693 nuc3(i, Cart::xyzz, m) = PmB(0) * nuc3(i, Cart::yzz, m) -
694 PmC(0) * nuc3(i, Cart::yzz, m + 1) +
695 nx_i * fak *
696 (nuc3(ilx_i, Cart::yzz, m) -
697 nuc3(ilx_i, Cart::yzz, m + 1));
698 nuc3(i, Cart::xzzz, m) = PmB(0) * nuc3(i, Cart::zzz, m) -
699 PmC(0) * nuc3(i, Cart::zzz, m + 1) +
700 nx_i * fak *
701 (nuc3(ilx_i, Cart::zzz, m) -
702 nuc3(ilx_i, Cart::zzz, m + 1));
703 nuc3(i, Cart::yyyy, m) = PmB(1) * nuc3(i, Cart::yyy, m) -
704 PmC(1) * nuc3(i, Cart::yyy, m + 1) +
705 ny_i * fak *
706 (nuc3(ily_i, Cart::yyy, m) -
707 nuc3(ily_i, Cart::yyy, m + 1)) +
708 3 * term_yy;
709 nuc3(i, Cart::yyyz, m) = PmB(2) * nuc3(i, Cart::yyy, m) -
710 PmC(2) * nuc3(i, Cart::yyy, m + 1) +
711 nz_i * fak *
712 (nuc3(ilz_i, Cart::yyy, m) -
713 nuc3(ilz_i, Cart::yyy, m + 1));
714 nuc3(i, Cart::yyzz, m) = PmB(1) * nuc3(i, Cart::yzz, m) -
715 PmC(1) * nuc3(i, Cart::yzz, m + 1) +
716 ny_i * fak *
717 (nuc3(ily_i, Cart::yzz, m) -
718 nuc3(ily_i, Cart::yzz, m + 1)) +
719 term_zz;
720 nuc3(i, Cart::yzzz, m) = PmB(1) * nuc3(i, Cart::zzz, m) -
721 PmC(1) * nuc3(i, Cart::zzz, m + 1) +
722 ny_i * fak *
723 (nuc3(ily_i, Cart::zzz, m) -
724 nuc3(ily_i, Cart::zzz, m + 1));
725 nuc3(i, Cart::zzzz, m) = PmB(2) * nuc3(i, Cart::zzz, m) -
726 PmC(2) * nuc3(i, Cart::zzz, m + 1) +
727 nz_i * fak *
728 (nuc3(ilz_i, Cart::zzz, m) -
729 nuc3(ilz_i, Cart::zzz, m + 1)) +
730 3 * term_zz;
731 }
732 }
733 //------------------------------------------------------
734 } // end if (lmax_col > 3)
735
736 multipole.noalias() +=
737 charge * Eigen::Map<Eigen::MatrixXd>(nuc3.data(), nrows, ncols);
738
739 if (rank == 1) {
740 // no quadrupole: the dipole integrals are needed only contracted
741 ContractedDipoleIntegrals(nuc3, dipc, lmax_row, lmax_col, PmA, PmB,
742 PmC, fak, 2. * zeta * prefactor, FmU,
743 dipole);
744 multipole += Eigen::Map<Eigen::MatrixXd>(dipc.data(), nrows, ncols);
745 } else if (rank > 1) {
746
747 // (s-s element normiert )
748 double prefactor_dip_ = 2. * zeta * prefactor;
749 for (Index m = 0; m < lsum + 1; m++) {
750 dip4(0, 0, 0, m) = PmC(0) * prefactor_dip_ * FmU[m + 1];
751 dip4(0, 0, 1, m) = PmC(1) * prefactor_dip_ * FmU[m + 1];
752 dip4(0, 0, 2, m) = PmC(2) * prefactor_dip_ * FmU[m + 1];
753 }
754 //------------------------------------------------------
755
756 // Integrals p - s
757 if (lmax_row > 0) {
758 for (Index m = 0; m < lsum; m++) {
759 for (Index k = 0; k < 3; k++) {
760 dip4(Cart::x, 0, k, m) = PmA(0) * dip4(0, 0, k, m) -
761 PmC(0) * dip4(0, 0, k, m + 1) +
762 (k == 0) * nuc3(0, 0, m + 1);
763 dip4(Cart::y, 0, k, m) = PmA(1) * dip4(0, 0, k, m) -
764 PmC(1) * dip4(0, 0, k, m + 1) +
765 (k == 1) * nuc3(0, 0, m + 1);
766 dip4(Cart::z, 0, k, m) = PmA(2) * dip4(0, 0, k, m) -
767 PmC(2) * dip4(0, 0, k, m + 1) +
768 (k == 2) * nuc3(0, 0, m + 1);
769 }
770 }
771 }
772 //------------------------------------------------------
773
774 // Integrals d - s
775 if (lmax_row > 1) {
776 for (Index m = 0; m < lsum - 1; m++) {
777 for (Index k = 0; k < 3; k++) {
778 double term = fak * (dip4(0, 0, k, m) - dip4(0, 0, k, m + 1));
779 dip4(Cart::xx, 0, k, m) = PmA(0) * dip4(Cart::x, 0, k, m) -
780 PmC(0) * dip4(Cart::x, 0, k, m + 1) +
781 (k == 0) * nuc3(Cart::x, 0, m + 1) +
782 term;
783 dip4(Cart::xy, 0, k, m) = PmA(0) * dip4(Cart::y, 0, k, m) -
784 PmC(0) * dip4(Cart::y, 0, k, m + 1) +
785 (k == 0) * nuc3(Cart::y, 0, m + 1);
786 dip4(Cart::xz, 0, k, m) = PmA(0) * dip4(Cart::z, 0, k, m) -
787 PmC(0) * dip4(Cart::z, 0, k, m + 1) +
788 (k == 0) * nuc3(Cart::z, 0, m + 1);
789 dip4(Cart::yy, 0, k, m) = PmA(1) * dip4(Cart::y, 0, k, m) -
790 PmC(1) * dip4(Cart::y, 0, k, m + 1) +
791 (k == 1) * nuc3(Cart::y, 0, m + 1) +
792 term;
793 dip4(Cart::yz, 0, k, m) = PmA(1) * dip4(Cart::z, 0, k, m) -
794 PmC(1) * dip4(Cart::z, 0, k, m + 1) +
795 (k == 1) * nuc3(Cart::z, 0, m + 1);
796 dip4(Cart::zz, 0, k, m) = PmA(2) * dip4(Cart::z, 0, k, m) -
797 PmC(2) * dip4(Cart::z, 0, k, m + 1) +
798 (k == 2) * nuc3(Cart::z, 0, m + 1) +
799 term;
800 }
801 }
802 }
803 //------------------------------------------------------
804
805 // Integrals f - s
806 if (lmax_row > 2) {
807 for (Index m = 0; m < lsum - 2; m++) {
808 for (Index k = 0; k < 3; k++) {
809 dip4(Cart::xxx, 0, k, m) =
810 PmA(0) * dip4(Cart::xx, 0, k, m) -
811 PmC(0) * dip4(Cart::xx, 0, k, m + 1) +
812 (k == 0) * nuc3(Cart::xx, 0, m + 1) +
813 2 * fak *
814 (dip4(Cart::x, 0, k, m) - dip4(Cart::x, 0, k, m + 1));
815 dip4(Cart::xxy, 0, k, m) =
816 PmA(1) * dip4(Cart::xx, 0, k, m) -
817 PmC(1) * dip4(Cart::xx, 0, k, m + 1) +
818 (k == 1) * nuc3(Cart::xx, 0, m + 1);
819 dip4(Cart::xxz, 0, k, m) =
820 PmA(2) * dip4(Cart::xx, 0, k, m) -
821 PmC(2) * dip4(Cart::xx, 0, k, m + 1) +
822 (k == 2) * nuc3(Cart::xx, 0, m + 1);
823 dip4(Cart::xyy, 0, k, m) =
824 PmA(0) * dip4(Cart::yy, 0, k, m) -
825 PmC(0) * dip4(Cart::yy, 0, k, m + 1) +
826 (k == 0) * nuc3(Cart::yy, 0, m + 1);
827 dip4(Cart::xyz, 0, k, m) =
828 PmA(0) * dip4(Cart::yz, 0, k, m) -
829 PmC(0) * dip4(Cart::yz, 0, k, m + 1) +
830 (k == 0) * nuc3(Cart::yz, 0, m + 1);
831 dip4(Cart::xzz, 0, k, m) =
832 PmA(0) * dip4(Cart::zz, 0, k, m) -
833 PmC(0) * dip4(Cart::zz, 0, k, m + 1) +
834 (k == 0) * nuc3(Cart::zz, 0, m + 1);
835 dip4(Cart::yyy, 0, k, m) =
836 PmA(1) * dip4(Cart::yy, 0, k, m) -
837 PmC(1) * dip4(Cart::yy, 0, k, m + 1) +
838 (k == 1) * nuc3(Cart::yy, 0, m + 1) +
839 2 * fak *
840 (dip4(Cart::y, 0, k, m) - dip4(Cart::y, 0, k, m + 1));
841 dip4(Cart::yyz, 0, k, m) =
842 PmA(2) * dip4(Cart::yy, 0, k, m) -
843 PmC(2) * dip4(Cart::yy, 0, k, m + 1) +
844 (k == 2) * nuc3(Cart::yy, 0, m + 1);
845 dip4(Cart::yzz, 0, k, m) =
846 PmA(1) * dip4(Cart::zz, 0, k, m) -
847 PmC(1) * dip4(Cart::zz, 0, k, m + 1) +
848 (k == 1) * nuc3(Cart::zz, 0, m + 1);
849 dip4(Cart::zzz, 0, k, m) =
850 PmA(2) * dip4(Cart::zz, 0, k, m) -
851 PmC(2) * dip4(Cart::zz, 0, k, m + 1) +
852 (k == 2) * nuc3(Cart::zz, 0, m + 1) +
853 2 * fak *
854 (dip4(Cart::z, 0, k, m) - dip4(Cart::z, 0, k, m + 1));
855 }
856 }
857 }
858 //------------------------------------------------------
859
860 // Integrals g - s
861 if (lmax_row > 3) {
862 for (Index m = 0; m < lsum - 3; m++) {
863 for (Index k = 0; k < 3; k++) {
864 double term_xx = fak * (dip4(Cart::xx, 0, k, m) -
865 dip4(Cart::xx, 0, k, m + 1));
866 double term_yy = fak * (dip4(Cart::yy, 0, k, m) -
867 dip4(Cart::yy, 0, k, m + 1));
868 double term_zz = fak * (dip4(Cart::zz, 0, k, m) -
869 dip4(Cart::zz, 0, k, m + 1));
870 dip4(Cart::xxxx, 0, k, m) =
871 PmA(0) * dip4(Cart::xxx, 0, k, m) -
872 PmC(0) * dip4(Cart::xxx, 0, k, m + 1) +
873 (k == 0) * nuc3(Cart::xxx, 0, m + 1) + 3 * term_xx;
874 dip4(Cart::xxxy, 0, k, m) =
875 PmA(1) * dip4(Cart::xxx, 0, k, m) -
876 PmC(1) * dip4(Cart::xxx, 0, k, m + 1) +
877 (k == 1) * nuc3(Cart::xxx, 0, m + 1);
878 dip4(Cart::xxxz, 0, k, m) =
879 PmA(2) * dip4(Cart::xxx, 0, k, m) -
880 PmC(2) * dip4(Cart::xxx, 0, k, m + 1) +
881 (k == 2) * nuc3(Cart::xxx, 0, m + 1);
882 dip4(Cart::xxyy, 0, k, m) =
883 PmA(0) * dip4(Cart::xyy, 0, k, m) -
884 PmC(0) * dip4(Cart::xyy, 0, k, m + 1) +
885 (k == 0) * nuc3(Cart::xyy, 0, m + 1) + term_yy;
886 dip4(Cart::xxyz, 0, k, m) =
887 PmA(1) * dip4(Cart::xxz, 0, k, m) -
888 PmC(1) * dip4(Cart::xxz, 0, k, m + 1) +
889 (k == 1) * nuc3(Cart::xxz, 0, m + 1);
890 dip4(Cart::xxzz, 0, k, m) =
891 PmA(0) * dip4(Cart::xzz, 0, k, m) -
892 PmC(0) * dip4(Cart::xzz, 0, k, m + 1) +
893 (k == 0) * nuc3(Cart::xzz, 0, m + 1) + term_zz;
894 dip4(Cart::xyyy, 0, k, m) =
895 PmA(0) * dip4(Cart::yyy, 0, k, m) -
896 PmC(0) * dip4(Cart::yyy, 0, k, m + 1) +
897 (k == 0) * nuc3(Cart::yyy, 0, m + 1);
898 dip4(Cart::xyyz, 0, k, m) =
899 PmA(0) * dip4(Cart::yyz, 0, k, m) -
900 PmC(0) * dip4(Cart::yyz, 0, k, m + 1) +
901 (k == 0) * nuc3(Cart::yyz, 0, m + 1);
902 dip4(Cart::xyzz, 0, k, m) =
903 PmA(0) * dip4(Cart::yzz, 0, k, m) -
904 PmC(0) * dip4(Cart::yzz, 0, k, m + 1) +
905 (k == 0) * nuc3(Cart::yzz, 0, m + 1);
906 dip4(Cart::xzzz, 0, k, m) =
907 PmA(0) * dip4(Cart::zzz, 0, k, m) -
908 PmC(0) * dip4(Cart::zzz, 0, k, m + 1) +
909 (k == 0) * nuc3(Cart::zzz, 0, m + 1);
910 dip4(Cart::yyyy, 0, k, m) =
911 PmA(1) * dip4(Cart::yyy, 0, k, m) -
912 PmC(1) * dip4(Cart::yyy, 0, k, m + 1) +
913 (k == 1) * nuc3(Cart::yyy, 0, m + 1) + 3 * term_yy;
914 dip4(Cart::yyyz, 0, k, m) =
915 PmA(2) * dip4(Cart::yyy, 0, k, m) -
916 PmC(2) * dip4(Cart::yyy, 0, k, m + 1) +
917 (k == 2) * nuc3(Cart::yyy, 0, m + 1);
918 dip4(Cart::yyzz, 0, k, m) =
919 PmA(1) * dip4(Cart::yzz, 0, k, m) -
920 PmC(1) * dip4(Cart::yzz, 0, k, m + 1) +
921 (k == 1) * nuc3(Cart::yzz, 0, m + 1) + term_zz;
922 dip4(Cart::yzzz, 0, k, m) =
923 PmA(1) * dip4(Cart::zzz, 0, k, m) -
924 PmC(1) * dip4(Cart::zzz, 0, k, m + 1) +
925 (k == 1) * nuc3(Cart::zzz, 0, m + 1);
926 dip4(Cart::zzzz, 0, k, m) =
927 PmA(2) * dip4(Cart::zzz, 0, k, m) -
928 PmC(2) * dip4(Cart::zzz, 0, k, m + 1) +
929 (k == 2) * nuc3(Cart::zzz, 0, m + 1) + 3 * term_zz;
930 }
931 }
932 }
933 //------------------------------------------------------
934
935 if (lmax_col > 0) {
936
937 // Integrals s - p
938 for (Index m = 0; m < lmax_col; m++) {
939 for (Index k = 0; k < 3; k++) {
940 dip4(0, Cart::x, k, m) = PmB(0) * dip4(0, 0, k, m) -
941 PmC(0) * dip4(0, 0, k, m + 1) +
942 (k == 0) * nuc3(0, 0, m + 1);
943 dip4(0, Cart::y, k, m) = PmB(1) * dip4(0, 0, k, m) -
944 PmC(1) * dip4(0, 0, k, m + 1) +
945 (k == 1) * nuc3(0, 0, m + 1);
946 dip4(0, Cart::z, k, m) = PmB(2) * dip4(0, 0, k, m) -
947 PmC(2) * dip4(0, 0, k, m + 1) +
948 (k == 2) * nuc3(0, 0, m + 1);
949 }
950 }
951 //------------------------------------------------------
952
953 // Integrals p - p
954 if (lmax_row > 0) {
955 for (Index m = 0; m < lmax_col; m++) {
956 for (Index i = 1; i < 4; i++) {
957 for (Index k = 0; k < 3; k++) {
958 double term =
959 fak * (dip4(0, 0, k, m) - dip4(0, 0, k, m + 1));
960 dip4(i, Cart::x, k, m) = PmB(0) * dip4(i, 0, k, m) -
961 PmC(0) * dip4(i, 0, k, m + 1) +
962 (k == 0) * nuc3(i, 0, m + 1) +
963 nx[i] * term;
964 dip4(i, Cart::y, k, m) = PmB(1) * dip4(i, 0, k, m) -
965 PmC(1) * dip4(i, 0, k, m + 1) +
966 (k == 1) * nuc3(i, 0, m + 1) +
967 ny[i] * term;
968 dip4(i, Cart::z, k, m) = PmB(2) * dip4(i, 0, k, m) -
969 PmC(2) * dip4(i, 0, k, m + 1) +
970 (k == 2) * nuc3(i, 0, m + 1) +
971 nz[i] * term;
972 }
973 }
974 }
975 }
976 //------------------------------------------------------
977
978 // Integrals d - p f - p g - p
979 for (Index m = 0; m < lmax_col; m++) {
980 for (Index i = 4; i < n_orbitals[lmax_row]; i++) {
981 int nx_i = nx[i];
982 int ny_i = ny[i];
983 int nz_i = nz[i];
984 int ilx_i = i_less_x[i];
985 int ily_i = i_less_y[i];
986 int ilz_i = i_less_z[i];
987 for (Index k = 0; k < 3; k++) {
988 dip4(i, Cart::x, k, m) =
989 PmB(0) * dip4(i, 0, k, m) -
990 PmC(0) * dip4(i, 0, k, m + 1) +
991 (k == 0) * nuc3(i, 0, m + 1) +
992 nx_i * fak *
993 (dip4(ilx_i, 0, k, m) - dip4(ilx_i, 0, k, m + 1));
994 dip4(i, Cart::y, k, m) =
995 PmB(1) * dip4(i, 0, k, m) -
996 PmC(1) * dip4(i, 0, k, m + 1) +
997 (k == 1) * nuc3(i, 0, m + 1) +
998 ny_i * fak *
999 (dip4(ily_i, 0, k, m) - dip4(ily_i, 0, k, m + 1));
1000 dip4(i, Cart::z, k, m) =
1001 PmB(2) * dip4(i, 0, k, m) -
1002 PmC(2) * dip4(i, 0, k, m + 1) +
1003 (k == 2) * nuc3(i, 0, m + 1) +
1004 nz_i * fak *
1005 (dip4(ilz_i, 0, k, m) - dip4(ilz_i, 0, k, m + 1));
1006 }
1007 }
1008 }
1009 //------------------------------------------------------
1010
1011 } // end if (lmax_col > 0)
1012
1013 if (lmax_col > 1) {
1014
1015 // Integrals s - d
1016 for (Index m = 0; m < lmax_col - 1; m++) {
1017 for (Index k = 0; k < 3; k++) {
1018 double term = fak * (dip4(0, 0, k, m) - dip4(0, 0, k, m + 1));
1019 dip4(0, Cart::xx, k, m) = PmB(0) * dip4(0, Cart::x, k, m) -
1020 PmC(0) * dip4(0, Cart::x, k, m + 1) +
1021 (k == 0) * nuc3(0, Cart::x, m + 1) +
1022 term;
1023 dip4(0, Cart::xy, k, m) = PmB(0) * dip4(0, Cart::y, k, m) -
1024 PmC(0) * dip4(0, Cart::y, k, m + 1) +
1025 (k == 0) * nuc3(0, Cart::y, m + 1);
1026 dip4(0, Cart::xz, k, m) = PmB(0) * dip4(0, Cart::z, k, m) -
1027 PmC(0) * dip4(0, Cart::z, k, m + 1) +
1028 (k == 0) * nuc3(0, Cart::z, m + 1);
1029 dip4(0, Cart::yy, k, m) = PmB(1) * dip4(0, Cart::y, k, m) -
1030 PmC(1) * dip4(0, Cart::y, k, m + 1) +
1031 (k == 1) * nuc3(0, Cart::y, m + 1) +
1032 term;
1033 dip4(0, Cart::yz, k, m) = PmB(1) * dip4(0, Cart::z, k, m) -
1034 PmC(1) * dip4(0, Cart::z, k, m + 1) +
1035 (k == 1) * nuc3(0, Cart::z, m + 1);
1036 dip4(0, Cart::zz, k, m) = PmB(2) * dip4(0, Cart::z, k, m) -
1037 PmC(2) * dip4(0, Cart::z, k, m + 1) +
1038 (k == 2) * nuc3(0, Cart::z, m + 1) +
1039 term;
1040 }
1041 }
1042 //------------------------------------------------------
1043
1044 // Integrals p - d d - d f - d g - d
1045 for (Index m = 0; m < lmax_col - 1; m++) {
1046 for (Index i = 1; i < n_orbitals[lmax_row]; i++) {
1047 int nx_i = nx[i];
1048 int ny_i = ny[i];
1049 int nz_i = nz[i];
1050 int ilx_i = i_less_x[i];
1051 int ily_i = i_less_y[i];
1052 int ilz_i = i_less_z[i];
1053 for (Index k = 0; k < 3; k++) {
1054 double term = fak * (dip4(i, 0, k, m) - dip4(i, 0, k, m + 1));
1055 dip4(i, Cart::xx, k, m) =
1056 PmB(0) * dip4(i, Cart::x, k, m) -
1057 PmC(0) * dip4(i, Cart::x, k, m + 1) +
1058 (k == 0) * nuc3(i, Cart::x, m + 1) +
1059 nx_i * fak *
1060 (dip4(ilx_i, Cart::x, k, m) -
1061 dip4(ilx_i, Cart::x, k, m + 1)) +
1062 term;
1063 dip4(i, Cart::xy, k, m) =
1064 PmB(0) * dip4(i, Cart::y, k, m) -
1065 PmC(0) * dip4(i, Cart::y, k, m + 1) +
1066 (k == 0) * nuc3(i, Cart::y, m + 1) +
1067 nx_i * fak *
1068 (dip4(ilx_i, Cart::y, k, m) -
1069 dip4(ilx_i, Cart::y, k, m + 1));
1070 dip4(i, Cart::xz, k, m) =
1071 PmB(0) * dip4(i, Cart::z, k, m) -
1072 PmC(0) * dip4(i, Cart::z, k, m + 1) +
1073 (k == 0) * nuc3(i, Cart::z, m + 1) +
1074 nx_i * fak *
1075 (dip4(ilx_i, Cart::z, k, m) -
1076 dip4(ilx_i, Cart::z, k, m + 1));
1077 dip4(i, Cart::yy, k, m) =
1078 PmB(1) * dip4(i, Cart::y, k, m) -
1079 PmC(1) * dip4(i, Cart::y, k, m + 1) +
1080 (k == 1) * nuc3(i, Cart::y, m + 1) +
1081 ny_i * fak *
1082 (dip4(ily_i, Cart::y, k, m) -
1083 dip4(ily_i, Cart::y, k, m + 1)) +
1084 term;
1085 dip4(i, Cart::yz, k, m) =
1086 PmB(1) * dip4(i, Cart::z, k, m) -
1087 PmC(1) * dip4(i, Cart::z, k, m + 1) +
1088 (k == 1) * nuc3(i, Cart::z, m + 1) +
1089 ny_i * fak *
1090 (dip4(ily_i, Cart::z, k, m) -
1091 dip4(ily_i, Cart::z, k, m + 1));
1092 dip4(i, Cart::zz, k, m) =
1093 PmB(2) * dip4(i, Cart::z, k, m) -
1094 PmC(2) * dip4(i, Cart::z, k, m + 1) +
1095 (k == 2) * nuc3(i, Cart::z, m + 1) +
1096 nz_i * fak *
1097 (dip4(ilz_i, Cart::z, k, m) -
1098 dip4(ilz_i, Cart::z, k, m + 1)) +
1099 term;
1100 }
1101 }
1102 }
1103 //------------------------------------------------------
1104
1105 } // end if (lmax_col > 1)
1106
1107 if (lmax_col > 2) {
1108
1109 // Integrals s - f
1110 for (Index m = 0; m < lmax_col - 2; m++) {
1111 for (Index k = 0; k < 3; k++) {
1112 dip4(0, Cart::xxx, k, m) =
1113 PmB(0) * dip4(0, Cart::xx, k, m) -
1114 PmC(0) * dip4(0, Cart::xx, k, m + 1) +
1115 (k == 0) * nuc3(0, Cart::xx, m + 1) +
1116 2 * fak *
1117 (dip4(0, Cart::x, k, m) - dip4(0, Cart::x, k, m + 1));
1118 dip4(0, Cart::xxy, k, m) =
1119 PmB(1) * dip4(0, Cart::xx, k, m) -
1120 PmC(1) * dip4(0, Cart::xx, k, m + 1) +
1121 (k == 1) * nuc3(0, Cart::xx, m + 1);
1122 dip4(0, Cart::xxz, k, m) =
1123 PmB(2) * dip4(0, Cart::xx, k, m) -
1124 PmC(2) * dip4(0, Cart::xx, k, m + 1) +
1125 (k == 2) * nuc3(0, Cart::xx, m + 1);
1126 dip4(0, Cart::xyy, k, m) =
1127 PmB(0) * dip4(0, Cart::yy, k, m) -
1128 PmC(0) * dip4(0, Cart::yy, k, m + 1) +
1129 (k == 0) * nuc3(0, Cart::yy, m + 1);
1130 dip4(0, Cart::xyz, k, m) =
1131 PmB(0) * dip4(0, Cart::yz, k, m) -
1132 PmC(0) * dip4(0, Cart::yz, k, m + 1) +
1133 (k == 0) * nuc3(0, Cart::yz, m + 1);
1134 dip4(0, Cart::xzz, k, m) =
1135 PmB(0) * dip4(0, Cart::zz, k, m) -
1136 PmC(0) * dip4(0, Cart::zz, k, m + 1) +
1137 (k == 0) * nuc3(0, Cart::zz, m + 1);
1138 dip4(0, Cart::yyy, k, m) =
1139 PmB(1) * dip4(0, Cart::yy, k, m) -
1140 PmC(1) * dip4(0, Cart::yy, k, m + 1) +
1141 (k == 1) * nuc3(0, Cart::yy, m + 1) +
1142 2 * fak *
1143 (dip4(0, Cart::y, k, m) - dip4(0, Cart::y, k, m + 1));
1144 dip4(0, Cart::yyz, k, m) =
1145 PmB(2) * dip4(0, Cart::yy, k, m) -
1146 PmC(2) * dip4(0, Cart::yy, k, m + 1) +
1147 (k == 2) * nuc3(0, Cart::yy, m + 1);
1148 dip4(0, Cart::yzz, k, m) =
1149 PmB(1) * dip4(0, Cart::zz, k, m) -
1150 PmC(1) * dip4(0, Cart::zz, k, m + 1) +
1151 (k == 1) * nuc3(0, Cart::zz, m + 1);
1152 dip4(0, Cart::zzz, k, m) =
1153 PmB(2) * dip4(0, Cart::zz, k, m) -
1154 PmC(2) * dip4(0, Cart::zz, k, m + 1) +
1155 (k == 2) * nuc3(0, Cart::zz, m + 1) +
1156 2 * fak *
1157 (dip4(0, Cart::z, k, m) - dip4(0, Cart::z, k, m + 1));
1158 }
1159 }
1160 //------------------------------------------------------
1161
1162 // Integrals p - f d - f f - f g - f
1163 for (Index m = 0; m < lmax_col - 2; m++) {
1164 for (Index i = 1; i < n_orbitals[lmax_row]; i++) {
1165 int nx_i = nx[i];
1166 int ny_i = ny[i];
1167 int nz_i = nz[i];
1168 int ilx_i = i_less_x[i];
1169 int ily_i = i_less_y[i];
1170 int ilz_i = i_less_z[i];
1171 for (Index k = 0; k < 3; k++) {
1172 double term_x =
1173 2 * fak *
1174 (dip4(i, Cart::x, k, m) - dip4(i, Cart::x, k, m + 1));
1175 double term_y =
1176 2 * fak *
1177 (dip4(i, Cart::y, k, m) - dip4(i, Cart::y, k, m + 1));
1178 double term_z =
1179 2 * fak *
1180 (dip4(i, Cart::z, k, m) - dip4(i, Cart::z, k, m + 1));
1181 dip4(i, Cart::xxx, k, m) =
1182 PmB(0) * dip4(i, Cart::xx, k, m) -
1183 PmC(0) * dip4(i, Cart::xx, k, m + 1) +
1184 (k == 0) * nuc3(i, Cart::xx, m + 1) +
1185 nx_i * fak *
1186 (dip4(ilx_i, Cart::xx, k, m) -
1187 dip4(ilx_i, Cart::xx, k, m + 1)) +
1188 term_x;
1189 dip4(i, Cart::xxy, k, m) =
1190 PmB(1) * dip4(i, Cart::xx, k, m) -
1191 PmC(1) * dip4(i, Cart::xx, k, m + 1) +
1192 (k == 1) * nuc3(i, Cart::xx, m + 1) +
1193 ny_i * fak *
1194 (dip4(ily_i, Cart::xx, k, m) -
1195 dip4(ily_i, Cart::xx, k, m + 1));
1196 dip4(i, Cart::xxz, k, m) =
1197 PmB(2) * dip4(i, Cart::xx, k, m) -
1198 PmC(2) * dip4(i, Cart::xx, k, m + 1) +
1199 (k == 2) * nuc3(i, Cart::xx, m + 1) +
1200 nz_i * fak *
1201 (dip4(ilz_i, Cart::xx, k, m) -
1202 dip4(ilz_i, Cart::xx, k, m + 1));
1203 dip4(i, Cart::xyy, k, m) =
1204 PmB(0) * dip4(i, Cart::yy, k, m) -
1205 PmC(0) * dip4(i, Cart::yy, k, m + 1) +
1206 (k == 0) * nuc3(i, Cart::yy, m + 1) +
1207 nx_i * fak *
1208 (dip4(ilx_i, Cart::yy, k, m) -
1209 dip4(ilx_i, Cart::yy, k, m + 1));
1210 dip4(i, Cart::xyz, k, m) =
1211 PmB(0) * dip4(i, Cart::yz, k, m) -
1212 PmC(0) * dip4(i, Cart::yz, k, m + 1) +
1213 (k == 0) * nuc3(i, Cart::yz, m + 1) +
1214 nx_i * fak *
1215 (dip4(ilx_i, Cart::yz, k, m) -
1216 dip4(ilx_i, Cart::yz, k, m + 1));
1217 dip4(i, Cart::xzz, k, m) =
1218 PmB(0) * dip4(i, Cart::zz, k, m) -
1219 PmC(0) * dip4(i, Cart::zz, k, m + 1) +
1220 (k == 0) * nuc3(i, Cart::zz, m + 1) +
1221 nx_i * fak *
1222 (dip4(ilx_i, Cart::zz, k, m) -
1223 dip4(ilx_i, Cart::zz, k, m + 1));
1224 dip4(i, Cart::yyy, k, m) =
1225 PmB(1) * dip4(i, Cart::yy, k, m) -
1226 PmC(1) * dip4(i, Cart::yy, k, m + 1) +
1227 (k == 1) * nuc3(i, Cart::yy, m + 1) +
1228 ny_i * fak *
1229 (dip4(ily_i, Cart::yy, k, m) -
1230 dip4(ily_i, Cart::yy, k, m + 1)) +
1231 term_y;
1232 dip4(i, Cart::yyz, k, m) =
1233 PmB(2) * dip4(i, Cart::yy, k, m) -
1234 PmC(2) * dip4(i, Cart::yy, k, m + 1) +
1235 (k == 2) * nuc3(i, Cart::yy, m + 1) +
1236 nz_i * fak *
1237 (dip4(ilz_i, Cart::yy, k, m) -
1238 dip4(ilz_i, Cart::yy, k, m + 1));
1239 dip4(i, Cart::yzz, k, m) =
1240 PmB(1) * dip4(i, Cart::zz, k, m) -
1241 PmC(1) * dip4(i, Cart::zz, k, m + 1) +
1242 (k == 1) * nuc3(i, Cart::zz, m + 1) +
1243 ny_i * fak *
1244 (dip4(ily_i, Cart::zz, k, m) -
1245 dip4(ily_i, Cart::zz, k, m + 1));
1246 dip4(i, Cart::zzz, k, m) =
1247 PmB(2) * dip4(i, Cart::zz, k, m) -
1248 PmC(2) * dip4(i, Cart::zz, k, m + 1) +
1249 (k == 2) * nuc3(i, Cart::zz, m + 1) +
1250 nz_i * fak *
1251 (dip4(ilz_i, Cart::zz, k, m) -
1252 dip4(ilz_i, Cart::zz, k, m + 1)) +
1253 term_z;
1254 }
1255 }
1256 }
1257 //------------------------------------------------------
1258
1259 } // end if (lmax_col > 2)
1260
1261 if (lmax_col > 3) {
1262
1263 // Integrals s - g
1264 for (Index m = 0; m < lmax_col - 3; m++) {
1265 for (Index k = 0; k < 3; k++) {
1266 double term_xx = fak * (dip4(0, Cart::xx, k, m) -
1267 dip4(0, Cart::xx, k, m + 1));
1268 double term_yy = fak * (dip4(0, Cart::yy, k, m) -
1269 dip4(0, Cart::yy, k, m + 1));
1270 double term_zz = fak * (dip4(0, Cart::zz, k, m) -
1271 dip4(0, Cart::zz, k, m + 1));
1272 dip4(0, Cart::xxxx, k, m) =
1273 PmB(0) * dip4(0, Cart::xxx, k, m) -
1274 PmC(0) * dip4(0, Cart::xxx, k, m + 1) +
1275 (k == 0) * nuc3(0, Cart::xxx, m + 1) + 3 * term_xx;
1276 dip4(0, Cart::xxxy, k, m) =
1277 PmB(1) * dip4(0, Cart::xxx, k, m) -
1278 PmC(1) * dip4(0, Cart::xxx, k, m + 1) +
1279 (k == 1) * nuc3(0, Cart::xxx, m + 1);
1280 dip4(0, Cart::xxxz, k, m) =
1281 PmB(2) * dip4(0, Cart::xxx, k, m) -
1282 PmC(2) * dip4(0, Cart::xxx, k, m + 1) +
1283 (k == 2) * nuc3(0, Cart::xxx, m + 1);
1284 dip4(0, Cart::xxyy, k, m) =
1285 PmB(0) * dip4(0, Cart::xyy, k, m) -
1286 PmC(0) * dip4(0, Cart::xyy, k, m + 1) +
1287 (k == 0) * nuc3(0, Cart::xyy, m + 1) + term_yy;
1288 dip4(0, Cart::xxyz, k, m) =
1289 PmB(1) * dip4(0, Cart::xxz, k, m) -
1290 PmC(1) * dip4(0, Cart::xxz, k, m + 1) +
1291 (k == 1) * nuc3(0, Cart::xxz, m + 1);
1292 dip4(0, Cart::xxzz, k, m) =
1293 PmB(0) * dip4(0, Cart::xzz, k, m) -
1294 PmC(0) * dip4(0, Cart::xzz, k, m + 1) +
1295 (k == 0) * nuc3(0, Cart::xzz, m + 1) + term_zz;
1296 dip4(0, Cart::xyyy, k, m) =
1297 PmB(0) * dip4(0, Cart::yyy, k, m) -
1298 PmC(0) * dip4(0, Cart::yyy, k, m + 1) +
1299 (k == 0) * nuc3(0, Cart::yyy, m + 1);
1300 dip4(0, Cart::xyyz, k, m) =
1301 PmB(0) * dip4(0, Cart::yyz, k, m) -
1302 PmC(0) * dip4(0, Cart::yyz, k, m + 1) +
1303 (k == 0) * nuc3(0, Cart::yyz, m + 1);
1304 dip4(0, Cart::xyzz, k, m) =
1305 PmB(0) * dip4(0, Cart::yzz, k, m) -
1306 PmC(0) * dip4(0, Cart::yzz, k, m + 1) +
1307 (k == 0) * nuc3(0, Cart::yzz, m + 1);
1308 dip4(0, Cart::xzzz, k, m) =
1309 PmB(0) * dip4(0, Cart::zzz, k, m) -
1310 PmC(0) * dip4(0, Cart::zzz, k, m + 1) +
1311 (k == 0) * nuc3(0, Cart::zzz, m + 1);
1312 dip4(0, Cart::yyyy, k, m) =
1313 PmB(1) * dip4(0, Cart::yyy, k, m) -
1314 PmC(1) * dip4(0, Cart::yyy, k, m + 1) +
1315 (k == 1) * nuc3(0, Cart::yyy, m + 1) + 3 * term_yy;
1316 dip4(0, Cart::yyyz, k, m) =
1317 PmB(2) * dip4(0, Cart::yyy, k, m) -
1318 PmC(2) * dip4(0, Cart::yyy, k, m + 1) +
1319 (k == 2) * nuc3(0, Cart::yyy, m + 1);
1320 dip4(0, Cart::yyzz, k, m) =
1321 PmB(1) * dip4(0, Cart::yzz, k, m) -
1322 PmC(1) * dip4(0, Cart::yzz, k, m + 1) +
1323 (k == 1) * nuc3(0, Cart::yzz, m + 1) + term_zz;
1324 dip4(0, Cart::yzzz, k, m) =
1325 PmB(1) * dip4(0, Cart::zzz, k, m) -
1326 PmC(1) * dip4(0, Cart::zzz, k, m + 1) +
1327 (k == 1) * nuc3(0, Cart::zzz, m + 1);
1328 dip4(0, Cart::zzzz, k, m) =
1329 PmB(2) * dip4(0, Cart::zzz, k, m) -
1330 PmC(2) * dip4(0, Cart::zzz, k, m + 1) +
1331 (k == 2) * nuc3(0, Cart::zzz, m + 1) + 3 * term_zz;
1332 }
1333 }
1334 //------------------------------------------------------
1335
1336 // Integrals p - g d - g f - g g - g
1337 for (Index m = 0; m < lmax_col - 3; m++) {
1338 for (Index i = 1; i < n_orbitals[lmax_row]; i++) {
1339 int nx_i = nx[i];
1340 int ny_i = ny[i];
1341 int nz_i = nz[i];
1342 int ilx_i = i_less_x[i];
1343 int ily_i = i_less_y[i];
1344 int ilz_i = i_less_z[i];
1345 for (Index k = 0; k < 3; k++) {
1346 double term_xx = fak * (dip4(i, Cart::xx, k, m) -
1347 dip4(i, Cart::xx, k, m + 1));
1348 double term_yy = fak * (dip4(i, Cart::yy, k, m) -
1349 dip4(i, Cart::yy, k, m + 1));
1350 double term_zz = fak * (dip4(i, Cart::zz, k, m) -
1351 dip4(i, Cart::zz, k, m + 1));
1352 dip4(i, Cart::xxxx, k, m) =
1353 PmB(0) * dip4(i, Cart::xxx, k, m) -
1354 PmC(0) * dip4(i, Cart::xxx, k, m + 1) +
1355 (k == 0) * nuc3(i, Cart::xxx, m + 1) +
1356 nx_i * fak *
1357 (dip4(ilx_i, Cart::xxx, k, m) -
1358 dip4(ilx_i, Cart::xxx, k, m + 1)) +
1359 3 * term_xx;
1360 dip4(i, Cart::xxxy, k, m) =
1361 PmB(1) * dip4(i, Cart::xxx, k, m) -
1362 PmC(1) * dip4(i, Cart::xxx, k, m + 1) +
1363 (k == 1) * nuc3(i, Cart::xxx, m + 1) +
1364 ny_i * fak *
1365 (dip4(ily_i, Cart::xxx, k, m) -
1366 dip4(ily_i, Cart::xxx, k, m + 1));
1367 dip4(i, Cart::xxxz, k, m) =
1368 PmB(2) * dip4(i, Cart::xxx, k, m) -
1369 PmC(2) * dip4(i, Cart::xxx, k, m + 1) +
1370 (k == 2) * nuc3(i, Cart::xxx, m + 1) +
1371 nz_i * fak *
1372 (dip4(ilz_i, Cart::xxx, k, m) -
1373 dip4(ilz_i, Cart::xxx, k, m + 1));
1374 dip4(i, Cart::xxyy, k, m) =
1375 PmB(0) * dip4(i, Cart::xyy, k, m) -
1376 PmC(0) * dip4(i, Cart::xyy, k, m + 1) +
1377 (k == 0) * nuc3(i, Cart::xyy, m + 1) +
1378 nx_i * fak *
1379 (dip4(ilx_i, Cart::xyy, k, m) -
1380 dip4(ilx_i, Cart::xyy, k, m + 1)) +
1381 term_yy;
1382 dip4(i, Cart::xxyz, k, m) =
1383 PmB(1) * dip4(i, Cart::xxz, k, m) -
1384 PmC(1) * dip4(i, Cart::xxz, k, m + 1) +
1385 (k == 1) * nuc3(i, Cart::xxz, m + 1) +
1386 ny_i * fak *
1387 (dip4(ily_i, Cart::xxz, k, m) -
1388 dip4(ily_i, Cart::xxz, k, m + 1));
1389 dip4(i, Cart::xxzz, k, m) =
1390 PmB(0) * dip4(i, Cart::xzz, k, m) -
1391 PmC(0) * dip4(i, Cart::xzz, k, m + 1) +
1392 (k == 0) * nuc3(i, Cart::xzz, m + 1) +
1393 nx_i * fak *
1394 (dip4(ilx_i, Cart::xzz, k, m) -
1395 dip4(ilx_i, Cart::xzz, k, m + 1)) +
1396 term_zz;
1397 dip4(i, Cart::xyyy, k, m) =
1398 PmB(0) * dip4(i, Cart::yyy, k, m) -
1399 PmC(0) * dip4(i, Cart::yyy, k, m + 1) +
1400 (k == 0) * nuc3(i, Cart::yyy, m + 1) +
1401 nx_i * fak *
1402 (dip4(ilx_i, Cart::yyy, k, m) -
1403 dip4(ilx_i, Cart::yyy, k, m + 1));
1404 dip4(i, Cart::xyyz, k, m) =
1405 PmB(0) * dip4(i, Cart::yyz, k, m) -
1406 PmC(0) * dip4(i, Cart::yyz, k, m + 1) +
1407 (k == 0) * nuc3(i, Cart::yyz, m + 1) +
1408 nx_i * fak *
1409 (dip4(ilx_i, Cart::yyz, k, m) -
1410 dip4(ilx_i, Cart::yyz, k, m + 1));
1411 dip4(i, Cart::xyzz, k, m) =
1412 PmB(0) * dip4(i, Cart::yzz, k, m) -
1413 PmC(0) * dip4(i, Cart::yzz, k, m + 1) +
1414 (k == 0) * nuc3(i, Cart::yzz, m + 1) +
1415 nx_i * fak *
1416 (dip4(ilx_i, Cart::yzz, k, m) -
1417 dip4(ilx_i, Cart::yzz, k, m + 1));
1418 dip4(i, Cart::xzzz, k, m) =
1419 PmB(0) * dip4(i, Cart::zzz, k, m) -
1420 PmC(0) * dip4(i, Cart::zzz, k, m + 1) +
1421 (k == 0) * nuc3(i, Cart::zzz, m + 1) +
1422 nx_i * fak *
1423 (dip4(ilx_i, Cart::zzz, k, m) -
1424 dip4(ilx_i, Cart::zzz, k, m + 1));
1425 dip4(i, Cart::yyyy, k, m) =
1426 PmB(1) * dip4(i, Cart::yyy, k, m) -
1427 PmC(1) * dip4(i, Cart::yyy, k, m + 1) +
1428 (k == 1) * nuc3(i, Cart::yyy, m + 1) +
1429 ny_i * fak *
1430 (dip4(ily_i, Cart::yyy, k, m) -
1431 dip4(ily_i, Cart::yyy, k, m + 1)) +
1432 3 * term_yy;
1433 dip4(i, Cart::yyyz, k, m) =
1434 PmB(2) * dip4(i, Cart::yyy, k, m) -
1435 PmC(2) * dip4(i, Cart::yyy, k, m + 1) +
1436 (k == 2) * nuc3(i, Cart::yyy, m + 1) +
1437 nz_i * fak *
1438 (dip4(ilz_i, Cart::yyy, k, m) -
1439 dip4(ilz_i, Cart::yyy, k, m + 1));
1440 dip4(i, Cart::yyzz, k, m) =
1441 PmB(1) * dip4(i, Cart::yzz, k, m) -
1442 PmC(1) * dip4(i, Cart::yzz, k, m + 1) +
1443 (k == 1) * nuc3(i, Cart::yzz, m + 1) +
1444 ny_i * fak *
1445 (dip4(ily_i, Cart::yzz, k, m) -
1446 dip4(ily_i, Cart::yzz, k, m + 1)) +
1447 term_zz;
1448 dip4(i, Cart::yzzz, k, m) =
1449 PmB(1) * dip4(i, Cart::zzz, k, m) -
1450 PmC(1) * dip4(i, Cart::zzz, k, m + 1) +
1451 (k == 1) * nuc3(i, Cart::zzz, m + 1) +
1452 ny_i * fak *
1453 (dip4(ily_i, Cart::zzz, k, m) -
1454 dip4(ily_i, Cart::zzz, k, m + 1));
1455 dip4(i, Cart::zzzz, k, m) =
1456 PmB(2) * dip4(i, Cart::zzz, k, m) -
1457 PmC(2) * dip4(i, Cart::zzz, k, m + 1) +
1458 (k == 2) * nuc3(i, Cart::zzz, m + 1) +
1459 nz_i * fak *
1460 (dip4(ilz_i, Cart::zzz, k, m) -
1461 dip4(ilz_i, Cart::zzz, k, m + 1)) +
1462 3 * term_zz;
1463 }
1464 }
1465 }
1466 //------------------------------------------------------
1467
1468 } // end if (lmax_col > 3)
1469
1470 multipole += dipole.x() *
1471 Eigen::Map<Eigen::MatrixXd>(dip4.data(), nrows, ncols);
1472 size_t offset = nrows * ncols;
1473 multipole += dipole.y() * Eigen::Map<Eigen::MatrixXd>(
1474 dip4.data() + offset, nrows, ncols);
1475 multipole += dipole.z() * Eigen::Map<Eigen::MatrixXd>(
1476 dip4.data() + 2 * offset, nrows, ncols);
1477
1478 if (rank > 1) {
1479
1480 double fact = 1. / 3.;
1481 std::array<double, 5> fac0 = {fact, fact, 0., 2. * fact, 0.};
1482 std::array<double, 5> fac1 = {fact, 0., fact, 0., 2. * fact};
1483 std::array<double, 5> fac2 = {0., fact, fact, -2. * fact,
1484 -2. * fact};
1485
1486 std::array<int, 5> ind0 = {1, 2, 0, 0, 0};
1487 std::array<int, 5> ind1 = {0, 0, 2, 0, 1};
1488 std::array<int, 5> ind2 = {0, 0, 1, 2, 2};
1489
1490 // (s-s element normiert )
1491 double prefactor_quad = (4. * zeta * zeta * prefactor) / 3.;
1492 for (Index m = 0; m < lsum + 1; m++) {
1493 quad4(0, 0, 0, m) = PmC(0) * PmC(1) * prefactor_quad * FmU[m + 2];
1494 quad4(0, 0, 1, m) = PmC(0) * PmC(2) * prefactor_quad * FmU[m + 2];
1495 quad4(0, 0, 2, m) = PmC(1) * PmC(2) * prefactor_quad * FmU[m + 2];
1496 quad4(0, 0, 3, m) = (PmC(0) * PmC(0) - PmC(2) * PmC(2)) *
1497 prefactor_quad * FmU[m + 2];
1498 quad4(0, 0, 4, m) = (PmC(1) * PmC(1) - PmC(2) * PmC(2)) *
1499 prefactor_quad * FmU[m + 2];
1500 }
1501 //------------------------------------------------------
1502
1503 // Integrals p - s
1504 if (lmax_row > 0) {
1505 for (Index m = 0; m < lsum; m++) {
1506 for (Index k = 0; k < 5; k++) {
1507 quad4(Cart::x, 0, k, m) =
1508 PmA(0) * quad4(0, 0, k, m) -
1509 PmC(0) * quad4(0, 0, k, m + 1) +
1510 fac0[k] * dip4(0, 0, ind0[k], m + 1);
1511 quad4(Cart::y, 0, k, m) =
1512 PmA(1) * quad4(0, 0, k, m) -
1513 PmC(1) * quad4(0, 0, k, m + 1) +
1514 fac1[k] * dip4(0, 0, ind1[k], m + 1);
1515 quad4(Cart::z, 0, k, m) =
1516 PmA(2) * quad4(0, 0, k, m) -
1517 PmC(2) * quad4(0, 0, k, m + 1) +
1518 fac2[k] * dip4(0, 0, ind2[k], m + 1);
1519 }
1520 }
1521 }
1522 //------------------------------------------------------
1523
1524 // Integrals d - s
1525 if (lmax_row > 1) {
1526 for (Index m = 0; m < lsum - 1; m++) {
1527 for (Index k = 0; k < 5; k++) {
1528 double term =
1529 fak * (quad4(0, 0, k, m) - quad4(0, 0, k, m + 1));
1530 quad4(Cart::xx, 0, k, m) =
1531 PmA(0) * quad4(Cart::x, 0, k, m) -
1532 PmC(0) * quad4(Cart::x, 0, k, m + 1) +
1533 fac0[k] * dip4(Cart::x, 0, ind0[k], m + 1) + term;
1534 quad4(Cart::xy, 0, k, m) =
1535 PmA(0) * quad4(Cart::y, 0, k, m) -
1536 PmC(0) * quad4(Cart::y, 0, k, m + 1) +
1537 fac0[k] * dip4(Cart::y, 0, ind0[k], m + 1);
1538 quad4(Cart::xz, 0, k, m) =
1539 PmA(0) * quad4(Cart::z, 0, k, m) -
1540 PmC(0) * quad4(Cart::z, 0, k, m + 1) +
1541 fac0[k] * dip4(Cart::z, 0, ind0[k], m + 1);
1542 quad4(Cart::yy, 0, k, m) =
1543 PmA(1) * quad4(Cart::y, 0, k, m) -
1544 PmC(1) * quad4(Cart::y, 0, k, m + 1) +
1545 fac1[k] * dip4(Cart::y, 0, ind1[k], m + 1) + term;
1546 quad4(Cart::yz, 0, k, m) =
1547 PmA(1) * quad4(Cart::z, 0, k, m) -
1548 PmC(1) * quad4(Cart::z, 0, k, m + 1) +
1549 fac1[k] * dip4(Cart::z, 0, ind1[k], m + 1);
1550 quad4(Cart::zz, 0, k, m) =
1551 PmA(2) * quad4(Cart::z, 0, k, m) -
1552 PmC(2) * quad4(Cart::z, 0, k, m + 1) +
1553 fac2[k] * dip4(Cart::z, 0, ind2[k], m + 1) + term;
1554 }
1555 }
1556 }
1557 //------------------------------------------------------
1558
1559 // Integrals f - s
1560 if (lmax_row > 2) {
1561 for (Index m = 0; m < lsum - 2; m++) {
1562 for (Index k = 0; k < 5; k++) {
1563 quad4(Cart::xxx, 0, k, m) =
1564 PmA(0) * quad4(Cart::xx, 0, k, m) -
1565 PmC(0) * quad4(Cart::xx, 0, k, m + 1) +
1566 fac0[k] * dip4(Cart::xx, 0, ind0[k], m + 1) +
1567 2 * fak *
1568 (quad4(Cart::x, 0, k, m) -
1569 quad4(Cart::x, 0, k, m + 1));
1570 quad4(Cart::xxy, 0, k, m) =
1571 PmA(1) * quad4(Cart::xx, 0, k, m) -
1572 PmC(1) * quad4(Cart::xx, 0, k, m + 1) +
1573 fac1[k] * dip4(Cart::xx, 0, ind1[k], m + 1);
1574 quad4(Cart::xxz, 0, k, m) =
1575 PmA(2) * quad4(Cart::xx, 0, k, m) -
1576 PmC(2) * quad4(Cart::xx, 0, k, m + 1) +
1577 fac2[k] * dip4(Cart::xx, 0, ind2[k], m + 1);
1578 quad4(Cart::xyy, 0, k, m) =
1579 PmA(0) * quad4(Cart::yy, 0, k, m) -
1580 PmC(0) * quad4(Cart::yy, 0, k, m + 1) +
1581 fac0[k] * dip4(Cart::yy, 0, ind0[k], m + 1);
1582 quad4(Cart::xyz, 0, k, m) =
1583 PmA(0) * quad4(Cart::yz, 0, k, m) -
1584 PmC(0) * quad4(Cart::yz, 0, k, m + 1) +
1585 fac0[k] * dip4(Cart::yz, 0, ind0[k], m + 1);
1586 quad4(Cart::xzz, 0, k, m) =
1587 PmA(0) * quad4(Cart::zz, 0, k, m) -
1588 PmC(0) * quad4(Cart::zz, 0, k, m + 1) +
1589 fac0[k] * dip4(Cart::zz, 0, ind0[k], m + 1);
1590 quad4(Cart::yyy, 0, k, m) =
1591 PmA(1) * quad4(Cart::yy, 0, k, m) -
1592 PmC(1) * quad4(Cart::yy, 0, k, m + 1) +
1593 fac1[k] * dip4(Cart::yy, 0, ind1[k], m + 1) +
1594 2 * fak *
1595 (quad4(Cart::y, 0, k, m) -
1596 quad4(Cart::y, 0, k, m + 1));
1597 quad4(Cart::yyz, 0, k, m) =
1598 PmA(2) * quad4(Cart::yy, 0, k, m) -
1599 PmC(2) * quad4(Cart::yy, 0, k, m + 1) +
1600 fac2[k] * dip4(Cart::yy, 0, ind2[k], m + 1);
1601 quad4(Cart::yzz, 0, k, m) =
1602 PmA(1) * quad4(Cart::zz, 0, k, m) -
1603 PmC(1) * quad4(Cart::zz, 0, k, m + 1) +
1604 fac1[k] * dip4(Cart::zz, 0, ind1[k], m + 1);
1605 quad4(Cart::zzz, 0, k, m) =
1606 PmA(2) * quad4(Cart::zz, 0, k, m) -
1607 PmC(2) * quad4(Cart::zz, 0, k, m + 1) +
1608 fac2[k] * dip4(Cart::zz, 0, ind2[k], m + 1) +
1609 2 * fak *
1610 (quad4(Cart::z, 0, k, m) -
1611 quad4(Cart::z, 0, k, m + 1));
1612 }
1613 }
1614 }
1615 //------------------------------------------------------
1616
1617 // Integrals g - s
1618 if (lmax_row > 3) {
1619 for (Index m = 0; m < lsum - 3; m++) {
1620 for (Index k = 0; k < 5; k++) {
1621 double term_xx = fak * (quad4(Cart::xx, 0, k, m) -
1622 quad4(Cart::xx, 0, k, m + 1));
1623 double term_yy = fak * (quad4(Cart::yy, 0, k, m) -
1624 quad4(Cart::yy, 0, k, m + 1));
1625 double term_zz = fak * (quad4(Cart::zz, 0, k, m) -
1626 quad4(Cart::zz, 0, k, m + 1));
1627 quad4(Cart::xxxx, 0, k, m) =
1628 PmA(0) * quad4(Cart::xxx, 0, k, m) -
1629 PmC(0) * quad4(Cart::xxx, 0, k, m + 1) +
1630 fac0[k] * dip4(Cart::xxx, 0, ind0[k], m + 1) +
1631 3 * term_xx;
1632 quad4(Cart::xxxy, 0, k, m) =
1633 PmA(1) * quad4(Cart::xxx, 0, k, m) -
1634 PmC(1) * quad4(Cart::xxx, 0, k, m + 1) +
1635 fac1[k] * dip4(Cart::xxx, 0, ind1[k], m + 1);
1636 quad4(Cart::xxxz, 0, k, m) =
1637 PmA(2) * quad4(Cart::xxx, 0, k, m) -
1638 PmC(2) * quad4(Cart::xxx, 0, k, m + 1) +
1639 fac2[k] * dip4(Cart::xxx, 0, ind2[k], m + 1);
1640 quad4(Cart::xxyy, 0, k, m) =
1641 PmA(0) * quad4(Cart::xyy, 0, k, m) -
1642 PmC(0) * quad4(Cart::xyy, 0, k, m + 1) +
1643 fac0[k] * dip4(Cart::xyy, 0, ind0[k], m + 1) + term_yy;
1644 quad4(Cart::xxyz, 0, k, m) =
1645 PmA(1) * quad4(Cart::xxz, 0, k, m) -
1646 PmC(1) * quad4(Cart::xxz, 0, k, m + 1) +
1647 fac1[k] * dip4(Cart::xxz, 0, ind1[k], m + 1);
1648 quad4(Cart::xxzz, 0, k, m) =
1649 PmA(0) * quad4(Cart::xzz, 0, k, m) -
1650 PmC(0) * quad4(Cart::xzz, 0, k, m + 1) +
1651 fac0[k] * dip4(Cart::xzz, 0, ind0[k], m + 1) + term_zz;
1652 quad4(Cart::xyyy, 0, k, m) =
1653 PmA(0) * quad4(Cart::yyy, 0, k, m) -
1654 PmC(0) * quad4(Cart::yyy, 0, k, m + 1) +
1655 fac0[k] * dip4(Cart::yyy, 0, ind0[k], m + 1);
1656 quad4(Cart::xyyz, 0, k, m) =
1657 PmA(0) * quad4(Cart::yyz, 0, k, m) -
1658 PmC(0) * quad4(Cart::yyz, 0, k, m + 1) +
1659 fac0[k] * dip4(Cart::yyz, 0, ind0[k], m + 1);
1660 quad4(Cart::xyzz, 0, k, m) =
1661 PmA(0) * quad4(Cart::yzz, 0, k, m) -
1662 PmC(0) * quad4(Cart::yzz, 0, k, m + 1) +
1663 fac0[k] * dip4(Cart::yzz, 0, ind0[k], m + 1);
1664 quad4(Cart::xzzz, 0, k, m) =
1665 PmA(0) * quad4(Cart::zzz, 0, k, m) -
1666 PmC(0) * quad4(Cart::zzz, 0, k, m + 1) +
1667 fac0[k] * dip4(Cart::zzz, 0, ind0[k], m + 1);
1668 quad4(Cart::yyyy, 0, k, m) =
1669 PmA(1) * quad4(Cart::yyy, 0, k, m) -
1670 PmC(1) * quad4(Cart::yyy, 0, k, m + 1) +
1671 fac1[k] * dip4(Cart::yyy, 0, ind1[k], m + 1) +
1672 3 * term_yy;
1673 quad4(Cart::yyyz, 0, k, m) =
1674 PmA(2) * quad4(Cart::yyy, 0, k, m) -
1675 PmC(2) * quad4(Cart::yyy, 0, k, m + 1) +
1676 fac2[k] * dip4(Cart::yyy, 0, ind2[k], m + 1);
1677 quad4(Cart::yyzz, 0, k, m) =
1678 PmA(1) * quad4(Cart::yzz, 0, k, m) -
1679 PmC(1) * quad4(Cart::yzz, 0, k, m + 1) +
1680 fac1[k] * dip4(Cart::yzz, 0, ind1[k], m + 1) + term_zz;
1681 quad4(Cart::yzzz, 0, k, m) =
1682 PmA(1) * quad4(Cart::zzz, 0, k, m) -
1683 PmC(1) * quad4(Cart::zzz, 0, k, m + 1) +
1684 fac1[k] * dip4(Cart::zzz, 0, ind1[k], m + 1);
1685 quad4(Cart::zzzz, 0, k, m) =
1686 PmA(2) * quad4(Cart::zzz, 0, k, m) -
1687 PmC(2) * quad4(Cart::zzz, 0, k, m + 1) +
1688 fac2[k] * dip4(Cart::zzz, 0, ind2[k], m + 1) +
1689 3 * term_zz;
1690 }
1691 }
1692 }
1693 //------------------------------------------------------
1694
1695 if (lmax_col > 0) {
1696
1697 // Integrals s - p
1698 for (Index m = 0; m < lmax_col; m++) {
1699 for (Index k = 0; k < 5; k++) {
1700 quad4(0, Cart::x, k, m) =
1701 PmB(0) * quad4(0, 0, k, m) -
1702 PmC(0) * quad4(0, 0, k, m + 1) +
1703 fac0[k] * dip4(0, 0, ind0[k], m + 1);
1704 quad4(0, Cart::y, k, m) =
1705 PmB(1) * quad4(0, 0, k, m) -
1706 PmC(1) * quad4(0, 0, k, m + 1) +
1707 fac1[k] * dip4(0, 0, ind1[k], m + 1);
1708 quad4(0, Cart::z, k, m) =
1709 PmB(2) * quad4(0, 0, k, m) -
1710 PmC(2) * quad4(0, 0, k, m + 1) +
1711 fac2[k] * dip4(0, 0, ind2[k], m + 1);
1712 }
1713 }
1714 //------------------------------------------------------
1715
1716 // Integrals p - p
1717 if (lmax_row > 0) {
1718 for (Index m = 0; m < lmax_col; m++) {
1719 for (Index i = 1; i < 4; i++) {
1720 for (Index k = 0; k < 5; k++) {
1721 double term =
1722 fak * (quad4(0, 0, k, m) - quad4(0, 0, k, m + 1));
1723 quad4(i, Cart::x, k, m) =
1724 PmB(0) * quad4(i, 0, k, m) -
1725 PmC(0) * quad4(i, 0, k, m + 1) +
1726 fac0[k] * dip4(i, 0, ind0[k], m + 1) + nx[i] * term;
1727 quad4(i, Cart::y, k, m) =
1728 PmB(1) * quad4(i, 0, k, m) -
1729 PmC(1) * quad4(i, 0, k, m + 1) +
1730 fac1[k] * dip4(i, 0, ind1[k], m + 1) + ny[i] * term;
1731 quad4(i, Cart::z, k, m) =
1732 PmB(2) * quad4(i, 0, k, m) -
1733 PmC(2) * quad4(i, 0, k, m + 1) +
1734 fac2[k] * dip4(i, 0, ind2[k], m + 1) + nz[i] * term;
1735 }
1736 }
1737 }
1738 }
1739 //------------------------------------------------------
1740
1741 // Integrals d - p f - p g - p
1742 for (Index m = 0; m < lmax_col; m++) {
1743 for (Index i = 4; i < n_orbitals[lmax_row]; i++) {
1744 int nx_i = nx[i];
1745 int ny_i = ny[i];
1746 int nz_i = nz[i];
1747 int ilx_i = i_less_x[i];
1748 int ily_i = i_less_y[i];
1749 int ilz_i = i_less_z[i];
1750 for (Index k = 0; k < 5; k++) {
1751 quad4(i, Cart::x, k, m) =
1752 PmB(0) * quad4(i, 0, k, m) -
1753 PmC(0) * quad4(i, 0, k, m + 1) +
1754 fac0[k] * dip4(i, 0, ind0[k], m + 1) +
1755 nx_i * fak *
1756 (quad4(ilx_i, 0, k, m) - quad4(ilx_i, 0, k, m + 1));
1757 quad4(i, Cart::y, k, m) =
1758 PmB(1) * quad4(i, 0, k, m) -
1759 PmC(1) * quad4(i, 0, k, m + 1) +
1760 fac1[k] * dip4(i, 0, ind1[k], m + 1) +
1761 ny_i * fak *
1762 (quad4(ily_i, 0, k, m) - quad4(ily_i, 0, k, m + 1));
1763 quad4(i, Cart::z, k, m) =
1764 PmB(2) * quad4(i, 0, k, m) -
1765 PmC(2) * quad4(i, 0, k, m + 1) +
1766 fac2[k] * dip4(i, 0, ind2[k], m + 1) +
1767 nz_i * fak *
1768 (quad4(ilz_i, 0, k, m) - quad4(ilz_i, 0, k, m + 1));
1769 }
1770 }
1771 }
1772 //------------------------------------------------------
1773
1774 } // end if (lmax_col > 0)
1775
1776 if (lmax_col > 1) {
1777
1778 // Integrals s - d
1779 for (Index m = 0; m < lmax_col - 1; m++) {
1780 for (Index k = 0; k < 5; k++) {
1781 double term =
1782 fak * (quad4(0, 0, k, m) - quad4(0, 0, k, m + 1));
1783 quad4(0, Cart::xx, k, m) =
1784 PmB(0) * quad4(0, Cart::x, k, m) -
1785 PmC(0) * quad4(0, Cart::x, k, m + 1) +
1786 fac0[k] * dip4(0, Cart::x, ind0[k], m + 1) + term;
1787 quad4(0, Cart::xy, k, m) =
1788 PmB(0) * quad4(0, Cart::y, k, m) -
1789 PmC(0) * quad4(0, Cart::y, k, m + 1) +
1790 fac0[k] * dip4(0, Cart::y, ind0[k], m + 1);
1791 quad4(0, Cart::xz, k, m) =
1792 PmB(0) * quad4(0, Cart::z, k, m) -
1793 PmC(0) * quad4(0, Cart::z, k, m + 1) +
1794 fac0[k] * dip4(0, Cart::z, ind0[k], m + 1);
1795 quad4(0, Cart::yy, k, m) =
1796 PmB(1) * quad4(0, Cart::y, k, m) -
1797 PmC(1) * quad4(0, Cart::y, k, m + 1) +
1798 fac1[k] * dip4(0, Cart::y, ind1[k], m + 1) + term;
1799 quad4(0, Cart::yz, k, m) =
1800 PmB(1) * quad4(0, Cart::z, k, m) -
1801 PmC(1) * quad4(0, Cart::z, k, m + 1) +
1802 fac1[k] * dip4(0, Cart::z, ind1[k], m + 1);
1803 quad4(0, Cart::zz, k, m) =
1804 PmB(2) * quad4(0, Cart::z, k, m) -
1805 PmC(2) * quad4(0, Cart::z, k, m + 1) +
1806 fac2[k] * dip4(0, Cart::z, ind2[k], m + 1) + term;
1807 }
1808 }
1809 //------------------------------------------------------
1810
1811 // Integrals p - d d - d f - d g - d
1812 for (Index m = 0; m < lmax_col - 1; m++) {
1813 for (Index i = 1; i < n_orbitals[lmax_row]; i++) {
1814 int nx_i = nx[i];
1815 int ny_i = ny[i];
1816 int nz_i = nz[i];
1817 int ilx_i = i_less_x[i];
1818 int ily_i = i_less_y[i];
1819 int ilz_i = i_less_z[i];
1820 for (Index k = 0; k < 5; k++) {
1821 double term =
1822 fak * (quad4(i, 0, k, m) - quad4(i, 0, k, m + 1));
1823 quad4(i, Cart::xx, k, m) =
1824 PmB(0) * quad4(i, Cart::x, k, m) -
1825 PmC(0) * quad4(i, Cart::x, k, m + 1) +
1826 fac0[k] * dip4(i, Cart::x, ind0[k], m + 1) +
1827 nx_i * fak *
1828 (quad4(ilx_i, Cart::x, k, m) -
1829 quad4(ilx_i, Cart::x, k, m + 1)) +
1830 term;
1831 quad4(i, Cart::xy, k, m) =
1832 PmB(0) * quad4(i, Cart::y, k, m) -
1833 PmC(0) * quad4(i, Cart::y, k, m + 1) +
1834 fac0[k] * dip4(i, Cart::y, ind0[k], m + 1) +
1835 nx_i * fak *
1836 (quad4(ilx_i, Cart::y, k, m) -
1837 quad4(ilx_i, Cart::y, k, m + 1));
1838 quad4(i, Cart::xz, k, m) =
1839 PmB(0) * quad4(i, Cart::z, k, m) -
1840 PmC(0) * quad4(i, Cart::z, k, m + 1) +
1841 fac0[k] * dip4(i, Cart::z, ind0[k], m + 1) +
1842 nx_i * fak *
1843 (quad4(ilx_i, Cart::z, k, m) -
1844 quad4(ilx_i, Cart::z, k, m + 1));
1845 quad4(i, Cart::yy, k, m) =
1846 PmB(1) * quad4(i, Cart::y, k, m) -
1847 PmC(1) * quad4(i, Cart::y, k, m + 1) +
1848 fac1[k] * dip4(i, Cart::y, ind1[k], m + 1) +
1849 ny_i * fak *
1850 (quad4(ily_i, Cart::y, k, m) -
1851 quad4(ily_i, Cart::y, k, m + 1)) +
1852 term;
1853 quad4(i, Cart::yz, k, m) =
1854 PmB(1) * quad4(i, Cart::z, k, m) -
1855 PmC(1) * quad4(i, Cart::z, k, m + 1) +
1856 fac1[k] * dip4(i, Cart::z, ind1[k], m + 1) +
1857 ny_i * fak *
1858 (quad4(ily_i, Cart::z, k, m) -
1859 quad4(ily_i, Cart::z, k, m + 1));
1860 quad4(i, Cart::zz, k, m) =
1861 PmB(2) * quad4(i, Cart::z, k, m) -
1862 PmC(2) * quad4(i, Cart::z, k, m + 1) +
1863 fac2[k] * dip4(i, Cart::z, ind2[k], m + 1) +
1864 nz_i * fak *
1865 (quad4(ilz_i, Cart::z, k, m) -
1866 quad4(ilz_i, Cart::z, k, m + 1)) +
1867 term;
1868 }
1869 }
1870 }
1871 //------------------------------------------------------
1872
1873 } // end if (lmax_col > 1)
1874
1875 if (lmax_col > 2) {
1876
1877 // Integrals s - f
1878 for (Index m = 0; m < lmax_col - 2; m++) {
1879 for (Index k = 0; k < 5; k++) {
1880 quad4(0, Cart::xxx, k, m) =
1881 PmB(0) * quad4(0, Cart::xx, k, m) -
1882 PmC(0) * quad4(0, Cart::xx, k, m + 1) +
1883 fac0[k] * dip4(0, Cart::xx, ind0[k], m + 1) +
1884 2 * fak *
1885 (quad4(0, Cart::x, k, m) -
1886 quad4(0, Cart::x, k, m + 1));
1887 quad4(0, Cart::xxy, k, m) =
1888 PmB(1) * quad4(0, Cart::xx, k, m) -
1889 PmC(1) * quad4(0, Cart::xx, k, m + 1) +
1890 fac1[k] * dip4(0, Cart::xx, ind1[k], m + 1);
1891 quad4(0, Cart::xxz, k, m) =
1892 PmB(2) * quad4(0, Cart::xx, k, m) -
1893 PmC(2) * quad4(0, Cart::xx, k, m + 1) +
1894 fac2[k] * dip4(0, Cart::xx, ind2[k], m + 1);
1895 quad4(0, Cart::xyy, k, m) =
1896 PmB(0) * quad4(0, Cart::yy, k, m) -
1897 PmC(0) * quad4(0, Cart::yy, k, m + 1) +
1898 fac0[k] * dip4(0, Cart::yy, ind0[k], m + 1);
1899 quad4(0, Cart::xyz, k, m) =
1900 PmB(0) * quad4(0, Cart::yz, k, m) -
1901 PmC(0) * quad4(0, Cart::yz, k, m + 1) +
1902 fac0[k] * dip4(0, Cart::yz, ind0[k], m + 1);
1903 quad4(0, Cart::xzz, k, m) =
1904 PmB(0) * quad4(0, Cart::zz, k, m) -
1905 PmC(0) * quad4(0, Cart::zz, k, m + 1) +
1906 fac0[k] * dip4(0, Cart::zz, ind0[k], m + 1);
1907 quad4(0, Cart::yyy, k, m) =
1908 PmB(1) * quad4(0, Cart::yy, k, m) -
1909 PmC(1) * quad4(0, Cart::yy, k, m + 1) +
1910 fac1[k] * dip4(0, Cart::yy, ind1[k], m + 1) +
1911 2 * fak *
1912 (quad4(0, Cart::y, k, m) -
1913 quad4(0, Cart::y, k, m + 1));
1914 quad4(0, Cart::yyz, k, m) =
1915 PmB(2) * quad4(0, Cart::yy, k, m) -
1916 PmC(2) * quad4(0, Cart::yy, k, m + 1) +
1917 fac2[k] * dip4(0, Cart::yy, ind2[k], m + 1);
1918 quad4(0, Cart::yzz, k, m) =
1919 PmB(1) * quad4(0, Cart::zz, k, m) -
1920 PmC(1) * quad4(0, Cart::zz, k, m + 1) +
1921 fac1[k] * dip4(0, Cart::zz, ind1[k], m + 1);
1922 quad4(0, Cart::zzz, k, m) =
1923 PmB(2) * quad4(0, Cart::zz, k, m) -
1924 PmC(2) * quad4(0, Cart::zz, k, m + 1) +
1925 fac2[k] * dip4(0, Cart::zz, ind2[k], m + 1) +
1926 2 * fak *
1927 (quad4(0, Cart::z, k, m) -
1928 quad4(0, Cart::z, k, m + 1));
1929 }
1930 }
1931 //------------------------------------------------------
1932
1933 // Integrals p - f d - f f - f g - f
1934 for (Index m = 0; m < lmax_col - 2; m++) {
1935 for (Index i = 1; i < n_orbitals[lmax_row]; i++) {
1936 int nx_i = nx[i];
1937 int ny_i = ny[i];
1938 int nz_i = nz[i];
1939 int ilx_i = i_less_x[i];
1940 int ily_i = i_less_y[i];
1941 int ilz_i = i_less_z[i];
1942 for (Index k = 0; k < 5; k++) {
1943 double term_x =
1944 2 * fak *
1945 (quad4(i, Cart::x, k, m) - quad4(i, Cart::x, k, m + 1));
1946 double term_y =
1947 2 * fak *
1948 (quad4(i, Cart::y, k, m) - quad4(i, Cart::y, k, m + 1));
1949 double term_z =
1950 2 * fak *
1951 (quad4(i, Cart::z, k, m) - quad4(i, Cart::z, k, m + 1));
1952 quad4(i, Cart::xxx, k, m) =
1953 PmB(0) * quad4(i, Cart::xx, k, m) -
1954 PmC(0) * quad4(i, Cart::xx, k, m + 1) +
1955 fac0[k] * dip4(i, Cart::xx, ind0[k], m + 1) +
1956 nx_i * fak *
1957 (quad4(ilx_i, Cart::xx, k, m) -
1958 quad4(ilx_i, Cart::xx, k, m + 1)) +
1959 term_x;
1960 quad4(i, Cart::xxy, k, m) =
1961 PmB(1) * quad4(i, Cart::xx, k, m) -
1962 PmC(1) * quad4(i, Cart::xx, k, m + 1) +
1963 fac1[k] * dip4(i, Cart::xx, ind1[k], m + 1) +
1964 ny_i * fak *
1965 (quad4(ily_i, Cart::xx, k, m) -
1966 quad4(ily_i, Cart::xx, k, m + 1));
1967 quad4(i, Cart::xxz, k, m) =
1968 PmB(2) * quad4(i, Cart::xx, k, m) -
1969 PmC(2) * quad4(i, Cart::xx, k, m + 1) +
1970 fac2[k] * dip4(i, Cart::xx, ind2[k], m + 1) +
1971 nz_i * fak *
1972 (quad4(ilz_i, Cart::xx, k, m) -
1973 quad4(ilz_i, Cart::xx, k, m + 1));
1974 quad4(i, Cart::xyy, k, m) =
1975 PmB(0) * quad4(i, Cart::yy, k, m) -
1976 PmC(0) * quad4(i, Cart::yy, k, m + 1) +
1977 fac0[k] * dip4(i, Cart::yy, ind0[k], m + 1) +
1978 nx_i * fak *
1979 (quad4(ilx_i, Cart::yy, k, m) -
1980 quad4(ilx_i, Cart::yy, k, m + 1));
1981 quad4(i, Cart::xyz, k, m) =
1982 PmB(0) * quad4(i, Cart::yz, k, m) -
1983 PmC(0) * quad4(i, Cart::yz, k, m + 1) +
1984 fac0[k] * dip4(i, Cart::yz, ind0[k], m + 1) +
1985 nx_i * fak *
1986 (quad4(ilx_i, Cart::yz, k, m) -
1987 quad4(ilx_i, Cart::yz, k, m + 1));
1988 quad4(i, Cart::xzz, k, m) =
1989 PmB(0) * quad4(i, Cart::zz, k, m) -
1990 PmC(0) * quad4(i, Cart::zz, k, m + 1) +
1991 fac0[k] * dip4(i, Cart::zz, ind0[k], m + 1) +
1992 nx_i * fak *
1993 (quad4(ilx_i, Cart::zz, k, m) -
1994 quad4(ilx_i, Cart::zz, k, m + 1));
1995 quad4(i, Cart::yyy, k, m) =
1996 PmB(1) * quad4(i, Cart::yy, k, m) -
1997 PmC(1) * quad4(i, Cart::yy, k, m + 1) +
1998 fac1[k] * dip4(i, Cart::yy, ind1[k], m + 1) +
1999 ny_i * fak *
2000 (quad4(ily_i, Cart::yy, k, m) -
2001 quad4(ily_i, Cart::yy, k, m + 1)) +
2002 term_y;
2003 quad4(i, Cart::yyz, k, m) =
2004 PmB(2) * quad4(i, Cart::yy, k, m) -
2005 PmC(2) * quad4(i, Cart::yy, k, m + 1) +
2006 fac2[k] * dip4(i, Cart::yy, ind2[k], m + 1) +
2007 nz_i * fak *
2008 (quad4(ilz_i, Cart::yy, k, m) -
2009 quad4(ilz_i, Cart::yy, k, m + 1));
2010 quad4(i, Cart::yzz, k, m) =
2011 PmB(1) * quad4(i, Cart::zz, k, m) -
2012 PmC(1) * quad4(i, Cart::zz, k, m + 1) +
2013 fac1[k] * dip4(i, Cart::zz, ind1[k], m + 1) +
2014 ny_i * fak *
2015 (quad4(ily_i, Cart::zz, k, m) -
2016 quad4(ily_i, Cart::zz, k, m + 1));
2017 quad4(i, Cart::zzz, k, m) =
2018 PmB(2) * quad4(i, Cart::zz, k, m) -
2019 PmC(2) * quad4(i, Cart::zz, k, m + 1) +
2020 fac2[k] * dip4(i, Cart::zz, ind2[k], m + 1) +
2021 nz_i * fak *
2022 (quad4(ilz_i, Cart::zz, k, m) -
2023 quad4(ilz_i, Cart::zz, k, m + 1)) +
2024 term_z;
2025 }
2026 }
2027 }
2028 //------------------------------------------------------
2029
2030 } // end if (lmax_col > 2)
2031
2032 if (lmax_col > 3) {
2033
2034 // Integrals s - g
2035 for (Index m = 0; m < lmax_col - 3; m++) {
2036 for (Index k = 0; k < 5; k++) {
2037 double term_xx = fak * (quad4(0, Cart::xx, k, m) -
2038 quad4(0, Cart::xx, k, m + 1));
2039 double term_yy = fak * (quad4(0, Cart::yy, k, m) -
2040 quad4(0, Cart::yy, k, m + 1));
2041 double term_zz = fak * (quad4(0, Cart::zz, k, m) -
2042 quad4(0, Cart::zz, k, m + 1));
2043 quad4(0, Cart::xxxx, k, m) =
2044 PmB(0) * quad4(0, Cart::xxx, k, m) -
2045 PmC(0) * quad4(0, Cart::xxx, k, m + 1) +
2046 fac0[k] * dip4(0, Cart::xxx, ind0[k], m + 1) +
2047 3 * term_xx;
2048 quad4(0, Cart::xxxy, k, m) =
2049 PmB(1) * quad4(0, Cart::xxx, k, m) -
2050 PmC(1) * quad4(0, Cart::xxx, k, m + 1) +
2051 fac1[k] * dip4(0, Cart::xxx, ind1[k], m + 1);
2052 quad4(0, Cart::xxxz, k, m) =
2053 PmB(2) * quad4(0, Cart::xxx, k, m) -
2054 PmC(2) * quad4(0, Cart::xxx, k, m + 1) +
2055 fac2[k] * dip4(0, Cart::xxx, ind2[k], m + 1);
2056 quad4(0, Cart::xxyy, k, m) =
2057 PmB(0) * quad4(0, Cart::xyy, k, m) -
2058 PmC(0) * quad4(0, Cart::xyy, k, m + 1) +
2059 fac0[k] * dip4(0, Cart::xyy, ind0[k], m + 1) + term_yy;
2060 quad4(0, Cart::xxyz, k, m) =
2061 PmB(1) * quad4(0, Cart::xxz, k, m) -
2062 PmC(1) * quad4(0, Cart::xxz, k, m + 1) +
2063 fac1[k] * dip4(0, Cart::xxz, ind1[k], m + 1);
2064 quad4(0, Cart::xxzz, k, m) =
2065 PmB(0) * quad4(0, Cart::xzz, k, m) -
2066 PmC(0) * quad4(0, Cart::xzz, k, m + 1) +
2067 fac0[k] * dip4(0, Cart::xzz, ind0[k], m + 1) + term_zz;
2068 quad4(0, Cart::xyyy, k, m) =
2069 PmB(0) * quad4(0, Cart::yyy, k, m) -
2070 PmC(0) * quad4(0, Cart::yyy, k, m + 1) +
2071 fac0[k] * dip4(0, Cart::yyy, ind0[k], m + 1);
2072 quad4(0, Cart::xyyz, k, m) =
2073 PmB(0) * quad4(0, Cart::yyz, k, m) -
2074 PmC(0) * quad4(0, Cart::yyz, k, m + 1) +
2075 fac0[k] * dip4(0, Cart::yyz, ind0[k], m + 1);
2076 quad4(0, Cart::xyzz, k, m) =
2077 PmB(0) * quad4(0, Cart::yzz, k, m) -
2078 PmC(0) * quad4(0, Cart::yzz, k, m + 1) +
2079 fac0[k] * dip4(0, Cart::yzz, ind0[k], m + 1);
2080 quad4(0, Cart::xzzz, k, m) =
2081 PmB(0) * quad4(0, Cart::zzz, k, m) -
2082 PmC(0) * quad4(0, Cart::zzz, k, m + 1) +
2083 fac0[k] * dip4(0, Cart::zzz, ind0[k], m + 1);
2084 quad4(0, Cart::yyyy, k, m) =
2085 PmB(1) * quad4(0, Cart::yyy, k, m) -
2086 PmC(1) * quad4(0, Cart::yyy, k, m + 1) +
2087 fac1[k] * dip4(0, Cart::yyy, ind1[k], m + 1) +
2088 3 * term_yy;
2089 quad4(0, Cart::yyyz, k, m) =
2090 PmB(2) * quad4(0, Cart::yyy, k, m) -
2091 PmC(2) * quad4(0, Cart::yyy, k, m + 1) +
2092 fac2[k] * dip4(0, Cart::yyy, ind2[k], m + 1);
2093 quad4(0, Cart::yyzz, k, m) =
2094 PmB(1) * quad4(0, Cart::yzz, k, m) -
2095 PmC(1) * quad4(0, Cart::yzz, k, m + 1) +
2096 fac1[k] * dip4(0, Cart::yzz, ind1[k], m + 1) + term_zz;
2097 quad4(0, Cart::yzzz, k, m) =
2098 PmB(1) * quad4(0, Cart::zzz, k, m) -
2099 PmC(1) * quad4(0, Cart::zzz, k, m + 1) +
2100 fac1[k] * dip4(0, Cart::zzz, ind1[k], m + 1);
2101 quad4(0, Cart::zzzz, k, m) =
2102 PmB(2) * quad4(0, Cart::zzz, k, m) -
2103 PmC(2) * quad4(0, Cart::zzz, k, m + 1) +
2104 fac2[k] * dip4(0, Cart::zzz, ind2[k], m + 1) +
2105 3 * term_zz;
2106 }
2107 }
2108 //------------------------------------------------------
2109
2110 // Integrals p - g d - g f - g g - g
2111 for (Index m = 0; m < lmax_col - 3; m++) {
2112 for (Index i = 1; i < n_orbitals[lmax_row]; i++) {
2113 int nx_i = nx[i];
2114 int ny_i = ny[i];
2115 int nz_i = nz[i];
2116 int ilx_i = i_less_x[i];
2117 int ily_i = i_less_y[i];
2118 int ilz_i = i_less_z[i];
2119 for (Index k = 0; k < 5; k++) {
2120 double term_xx = fak * (quad4(i, Cart::xx, k, m) -
2121 quad4(i, Cart::xx, k, m + 1));
2122 double term_yy = fak * (quad4(i, Cart::yy, k, m) -
2123 quad4(i, Cart::yy, k, m + 1));
2124 double term_zz = fak * (quad4(i, Cart::zz, k, m) -
2125 quad4(i, Cart::zz, k, m + 1));
2126 quad4(i, Cart::xxxx, k, m) =
2127 PmB(0) * quad4(i, Cart::xxx, k, m) -
2128 PmC(0) * quad4(i, Cart::xxx, k, m + 1) +
2129 fac0[k] * dip4(i, Cart::xxx, ind0[k], m + 1) +
2130 nx_i * fak *
2131 (quad4(ilx_i, Cart::xxx, k, m) -
2132 quad4(ilx_i, Cart::xxx, k, m + 1)) +
2133 3 * term_xx;
2134 quad4(i, Cart::xxxy, k, m) =
2135 PmB(1) * quad4(i, Cart::xxx, k, m) -
2136 PmC(1) * quad4(i, Cart::xxx, k, m + 1) +
2137 fac1[k] * dip4(i, Cart::xxx, ind1[k], m + 1) +
2138 ny_i * fak *
2139 (quad4(ily_i, Cart::xxx, k, m) -
2140 quad4(ily_i, Cart::xxx, k, m + 1));
2141 quad4(i, Cart::xxxz, k, m) =
2142 PmB(2) * quad4(i, Cart::xxx, k, m) -
2143 PmC(2) * quad4(i, Cart::xxx, k, m + 1) +
2144 fac2[k] * dip4(i, Cart::xxx, ind2[k], m + 1) +
2145 nz_i * fak *
2146 (quad4(ilz_i, Cart::xxx, k, m) -
2147 quad4(ilz_i, Cart::xxx, k, m + 1));
2148 quad4(i, Cart::xxyy, k, m) =
2149 PmB(0) * quad4(i, Cart::xyy, k, m) -
2150 PmC(0) * quad4(i, Cart::xyy, k, m + 1) +
2151 fac0[k] * dip4(i, Cart::xyy, ind0[k], m + 1) +
2152 nx_i * fak *
2153 (quad4(ilx_i, Cart::xyy, k, m) -
2154 quad4(ilx_i, Cart::xyy, k, m + 1)) +
2155 term_yy;
2156 quad4(i, Cart::xxyz, k, m) =
2157 PmB(1) * quad4(i, Cart::xxz, k, m) -
2158 PmC(1) * quad4(i, Cart::xxz, k, m + 1) +
2159 fac1[k] * dip4(i, Cart::xxz, ind1[k], m + 1) +
2160 ny_i * fak *
2161 (quad4(ily_i, Cart::xxz, k, m) -
2162 quad4(ily_i, Cart::xxz, k, m + 1));
2163 quad4(i, Cart::xxzz, k, m) =
2164 PmB(0) * quad4(i, Cart::xzz, k, m) -
2165 PmC(0) * quad4(i, Cart::xzz, k, m + 1) +
2166 fac0[k] * dip4(i, Cart::xzz, ind0[k], m + 1) +
2167 nx_i * fak *
2168 (quad4(ilx_i, Cart::xzz, k, m) -
2169 quad4(ilx_i, Cart::xzz, k, m + 1)) +
2170 term_zz;
2171 quad4(i, Cart::xyyy, k, m) =
2172 PmB(0) * quad4(i, Cart::yyy, k, m) -
2173 PmC(0) * quad4(i, Cart::yyy, k, m + 1) +
2174 fac0[k] * dip4(i, Cart::yyy, ind0[k], m + 1) +
2175 nx_i * fak *
2176 (quad4(ilx_i, Cart::yyy, k, m) -
2177 quad4(ilx_i, Cart::yyy, k, m + 1));
2178 quad4(i, Cart::xyyz, k, m) =
2179 PmB(0) * quad4(i, Cart::yyz, k, m) -
2180 PmC(0) * quad4(i, Cart::yyz, k, m + 1) +
2181 fac0[k] * dip4(i, Cart::yyz, ind0[k], m + 1) +
2182 nx_i * fak *
2183 (quad4(ilx_i, Cart::yyz, k, m) -
2184 quad4(ilx_i, Cart::yyz, k, m + 1));
2185 quad4(i, Cart::xyzz, k, m) =
2186 PmB(0) * quad4(i, Cart::yzz, k, m) -
2187 PmC(0) * quad4(i, Cart::yzz, k, m + 1) +
2188 fac0[k] * dip4(i, Cart::yzz, ind0[k], m + 1) +
2189 nx_i * fak *
2190 (quad4(ilx_i, Cart::yzz, k, m) -
2191 quad4(ilx_i, Cart::yzz, k, m + 1));
2192 quad4(i, Cart::xzzz, k, m) =
2193 PmB(0) * quad4(i, Cart::zzz, k, m) -
2194 PmC(0) * quad4(i, Cart::zzz, k, m + 1) +
2195 fac0[k] * dip4(i, Cart::zzz, ind0[k], m + 1) +
2196 nx_i * fak *
2197 (quad4(ilx_i, Cart::zzz, k, m) -
2198 quad4(ilx_i, Cart::zzz, k, m + 1));
2199 quad4(i, Cart::yyyy, k, m) =
2200 PmB(1) * quad4(i, Cart::yyy, k, m) -
2201 PmC(1) * quad4(i, Cart::yyy, k, m + 1) +
2202 fac1[k] * dip4(i, Cart::yyy, ind1[k], m + 1) +
2203 ny_i * fak *
2204 (quad4(ily_i, Cart::yyy, k, m) -
2205 quad4(ily_i, Cart::yyy, k, m + 1)) +
2206 3 * term_yy;
2207 quad4(i, Cart::yyyz, k, m) =
2208 PmB(2) * quad4(i, Cart::yyy, k, m) -
2209 PmC(2) * quad4(i, Cart::yyy, k, m + 1) +
2210 fac2[k] * dip4(i, Cart::yyy, ind2[k], m + 1) +
2211 nz_i * fak *
2212 (quad4(ilz_i, Cart::yyy, k, m) -
2213 quad4(ilz_i, Cart::yyy, k, m + 1));
2214 quad4(i, Cart::yyzz, k, m) =
2215 PmB(1) * quad4(i, Cart::yzz, k, m) -
2216 PmC(1) * quad4(i, Cart::yzz, k, m + 1) +
2217 fac1[k] * dip4(i, Cart::yzz, ind1[k], m + 1) +
2218 ny_i * fak *
2219 (quad4(ily_i, Cart::yzz, k, m) -
2220 quad4(ily_i, Cart::yzz, k, m + 1)) +
2221 term_zz;
2222 quad4(i, Cart::yzzz, k, m) =
2223 PmB(1) * quad4(i, Cart::zzz, k, m) -
2224 PmC(1) * quad4(i, Cart::zzz, k, m + 1) +
2225 fac1[k] * dip4(i, Cart::zzz, ind1[k], m + 1) +
2226 ny_i * fak *
2227 (quad4(ily_i, Cart::zzz, k, m) -
2228 quad4(ily_i, Cart::zzz, k, m + 1));
2229 quad4(i, Cart::zzzz, k, m) =
2230 PmB(2) * quad4(i, Cart::zzz, k, m) -
2231 PmC(2) * quad4(i, Cart::zzz, k, m + 1) +
2232 fac2[k] * dip4(i, Cart::zzz, ind2[k], m + 1) +
2233 nz_i * fak *
2234 (quad4(ilz_i, Cart::zzz, k, m) -
2235 quad4(ilz_i, Cart::zzz, k, m + 1)) +
2236 3 * term_zz;
2237 }
2238 }
2239 }
2240 //------------------------------------------------------
2241
2242 } // end if (lmax_col > 3)
2243
2244 multipole += quadrupole(0, 1) * Eigen::Map<Eigen::MatrixXd>(
2245 quad4.data(), nrows, ncols);
2246 multipole +=
2247 quadrupole(0, 2) * Eigen::Map<Eigen::MatrixXd>(
2248 quad4.data() + offset, nrows, ncols);
2249 multipole +=
2250 quadrupole(1, 2) * Eigen::Map<Eigen::MatrixXd>(
2251 quad4.data() + 2 * offset, nrows, ncols);
2252 multipole += 0.5 * quadrupole(0, 0) *
2253 Eigen::Map<Eigen::MatrixXd>(quad4.data() + 3 * offset,
2254 nrows, ncols);
2255 multipole += 0.5 * quadrupole(1, 1) *
2256 Eigen::Map<Eigen::MatrixXd>(quad4.data() + 4 * offset,
2257 nrows, ncols);
2258 }
2259 }
2260
2261 } // sites
2262
2263 // save to matrix
2264 cartesian += AOTransform::getNorm(shell_row.getL(), gaussian_row) *
2265 AOTransform::getNorm(shell_col.getL(), gaussian_col) *
2266 multipole.bottomRightCorner(shell_row.getCartesianNumFunc(),
2267 shell_col.getCartesianNumFunc());
2268
2269 } // shell_col Gaussians
2270 } // shell_row Gaussians
2271
2272 matrix = AOTransform::tform(shell_row.getL(), shell_col.getL(), cartesian);
2273}
2274
2275void AOMultipole::setSites(const std::vector<const StaticSite*>& sites,
2276 Moments moments) {
2277 sites_.clear();
2278 sites_.reserve(sites.size());
2279 for (const StaticSite* site : sites) {
2280 SiteData d;
2281 d.pos = site->getPos();
2282 if (moments == Moments::Induced) {
2283 // A pure point dipole: no charge, no quadrupole.
2284 d.dipole = site->getInducedDipole();
2285 if (d.dipole.norm() <= 1e-12) {
2286 continue;
2287 }
2288 d.rank = 1;
2289 d.charge = 0.0;
2290 d.quadrupole = Eigen::Matrix3d::Zero();
2291 sites_.push_back(d);
2292 continue;
2293 }
2294 // getDipole() is virtual and includes the induced dipole of a PolarSite;
2295 // Q() holds the permanent moments only.
2296 d.dipole = (moments == Moments::All)
2297 ? site->getDipole()
2298 : Eigen::Vector3d(site->Q().segment<3>(1));
2299 d.rank = site->getRank();
2300 if (d.rank < 1 && d.dipole.norm() > 1e-12) {
2301 d.rank = 1;
2302 }
2303 d.charge = site->getCharge();
2304 // factor 1.5 I am not sure about but then 6 monopoles and this tensor
2305 // agree
2306 d.quadrupole = 1.5 * site->CalculateCartesianMultipole();
2307 sites_.push_back(d);
2308 }
2309}
2310
2311void AOMultipole::FillPotential(const AOBasis& aobasis,
2312 const Eigen::Vector3d& r) {
2313 StaticSite s = StaticSite(0, "", r);
2314 s.setCharge(1.0);
2315 setSites({&s});
2316 aopotential_ = Fill(aobasis);
2317}
2318
2319void AOMultipole::FillPotential(const AOBasis& aobasis,
2320 const QMMolecule& atoms) {
2321 std::vector<StaticSite> nuclei;
2322 nuclei.reserve(atoms.size());
2323 for (const auto& atom : atoms) {
2324 nuclei.emplace_back(atom, double(atom.getNuccharge()));
2325 }
2326 std::vector<const StaticSite*> sites;
2327 for (const StaticSite& s : nuclei) {
2328 sites.push_back(&s);
2329 }
2330 setSites(sites);
2331 aopotential_ = -Fill(aobasis);
2332}
2333
2334void AOMultipole::FillPotential(
2335 const AOBasis& aobasis,
2336 const std::vector<std::unique_ptr<StaticSite>>& externalsites,
2337 Moments moments) {
2338 std::vector<const StaticSite*> sites;
2339 sites.reserve(externalsites.size());
2340 for (const std::unique_ptr<StaticSite>& site : externalsites) {
2341 sites.push_back(site.get());
2342 }
2343 setSites(sites, moments);
2344 aopotential_ = -Fill(aobasis);
2345}
2346
2347} // namespace xtp
2348} // namespace votca
votca::xtp::AOBasis
Container to hold Basisfunctions for all atoms.
Definition aobasis.h:42
votca::xtp::AOMultipole::setSites
void setSites(const std::vector< const StaticSite * > &sites, Moments moments=Moments::All)
Definition aomultipole.cc:2275
votca::xtp::AOMultipole::FillBlock
void FillBlock(Eigen::Block< Eigen::MatrixXd > &matrix, const AOShell &shell_row, const AOShell &shell_col) const override
Definition aomultipole.cc:148
votca::xtp::AOMultipole::FillPotential
void FillPotential(const AOBasis &aobasis, const QMMolecule &atoms)
Definition aomultipole.cc:2319
votca::xtp::AOMultipole::Moments::All
@ All
permanent moments plus induced dipoles
Definition aopotential.h:71
votca::xtp::AOMultipole::Moments::Induced
@ Induced
induced dipoles only; sites without one are skipped
Definition aopotential.h:73
votca::xtp::AOMultipole::sites_
std::vector< SiteData > sites_
Definition aopotential.h:100
votca::xtp::AOPotential< double >::Fill
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > Fill(const AOBasis &aobasis) const
Definition aopotential.cc:27
votca::xtp::AOPotential< double >::aopotential_
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > aopotential_
Definition aopotential.h:49
votca::xtp::AOShell::getCartesianNumFunc
Index getCartesianNumFunc() const
Definition aoshell.h:104
votca::xtp::AOShell::getPos
const Eigen::Vector3d & getPos() const
Definition aoshell.h:113
votca::xtp::AOTransform::XIntegrate
static Eigen::VectorXd XIntegrate(Index size, double U)
Definition aotransform.cc:116
votca::xtp::AOTransform::i_less_y
static std::array< int, 165 > i_less_y()
Definition aotransform.cc:228
votca::xtp::AOTransform::i_less_z
static std::array< int, 165 > i_less_z()
Definition aotransform.cc:243
votca::xtp::AOTransform::getBlockSize
static Index getBlockSize(Index lmax)
Definition aotransform.cc:168
votca::xtp::AOTransform::tform
static Matrix tform(L l_row, L l_col, const Matrix &cartesian)
transforms a cartesian shell to a spherical cartesian shell
Definition aotransform.cc:68
votca::xtp::AOTransform::i_less_x
static std::array< int, 165 > i_less_x()
Definition aotransform.cc:214
votca::xtp::AOTransform::getNorm
static double getNorm(L l, const AOGaussianPrimitive &gaussian)
Definition aotransform.cc:40
votca::xtp::AOTransform::ny
static std::array< int, 165 > ny()
Definition aotransform.cc:192
votca::xtp::AOTransform::nz
static std::array< int, 165 > nz()
Definition aotransform.cc:203
votca::xtp::AOTransform::n_orbitals
static std::array< int, 9 > n_orbitals()
Definition aotransform.cc:177
votca::xtp::AOTransform::nx
static std::array< int, 165 > nx()
Definition aotransform.cc:181
votca::xtp::StaticSite
Class to represent Atom/Site in electrostatic.
Definition staticsite.h:37
votca::xtp
Charge transport classes.
Definition ERIs.h:28
votca
Provides a means for comparing floating point numbers.
Definition basebead.h:33
votca::Index
Eigen::Index Index
Definition types.h:26