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