votca 2024.1-dev
Loading...
Searching...
No Matches
aoplanewave.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
23
24namespace votca {
25namespace xtp {
26
27void AOPlanewave::FillBlock(Eigen::Block<Eigen::MatrixXcd>& matrix,
28 const AOShell& shell_row,
29 const AOShell& shell_col) const {
30
31 // shell info, only lmax tells how far to go
32 Index lmax_row = Index(shell_row.getL());
33 Index lmax_col = Index(shell_col.getL());
34 // set size of internal block for recursion
35 Index nrows = AOTransform::getBlockSize(lmax_row);
36 Index ncols = AOTransform::getBlockSize(lmax_col);
37 if (lmax_col > 6 || lmax_row > 6) {
38 throw std::runtime_error(
39 "Orbitals higher than i are not yet implemented. This should not have "
40 "happened!");
41 }
42 // get shell positions
43 const Eigen::Vector3d& pos_row = shell_row.getPos(); // get position R_{i}
44 const Eigen::Vector3d& pos_col = shell_col.getPos(); // get position R_{j}
45 const Eigen::Vector3d diff = pos_row - pos_col; // get difference r_{ij}
46 const double distsq = diff.squaredNorm(); // get |R_{ij}|^2
47 // get kvector modulus
48 const double kmodulus = k_.squaredNorm(); // get |k|^2
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 Eigen::MatrixXcd cartesian = Eigen::MatrixXcd::Zero(
59 shell_row.getCartesianNumFunc(), shell_col.getCartesianNumFunc());
60
61 // iterate over Gaussians in this shell_row
62 for (const auto& gaussian_row : shell_row) {
63 // iterate over Gaussians in this shell_col
64 // get decay constant
65 const double decay_row = gaussian_row.getDecay();
66
67 for (const auto& gaussian_col : shell_col) {
68 // get decay constant
69 const double decay_col = gaussian_col.getDecay();
70
71 // some helpers
72
73 const double fak = 0.5 / (decay_row + decay_col);
74 const double fak2 = 2.0 * fak;
75 double exparg = fak2 * decay_row * decay_col * distsq;
76
77 // initialize local matrix block for unnormalized cartesians
78 Eigen::MatrixXcd olk = Eigen::MatrixXcd::Zero(nrows, ncols);
79
80 using COMPLEX = std::complex<double>; // Define an abbreviation for
81 // complex numbers
82 Eigen::Vector3cd PmA;
83 PmA.real() = fak2 * (decay_row * pos_row + decay_col * pos_col) - pos_row;
84 PmA.imag() = fak * k_;
85 Eigen::Vector3cd PmB;
86 PmB.real() = fak2 * (decay_row * pos_row + decay_col * pos_col) - pos_col;
87 PmB.imag() = fak * k_;
88
89 const COMPLEX cfak(fak, 0.0);
90 const COMPLEX cfak2(fak2, 0.0);
91
92 // calculate s-s- overlap matrix element
93 COMPLEX ssol(std::pow(4.0 * decay_row * decay_col, 0.75) *
94 std::pow(fak2, 1.5) * std::exp(-exparg),
95 0.0); // s-s element
96
97 // calculate s-W-s matrix element
98 double kdotr_row = k_.dot(pos_row);
99 double kdotr_col = k_.dot(pos_col);
100 COMPLEX kexparg(
101 fak2 * (-0.25) * (kmodulus),
102 fak2 * (decay_row) * (kdotr_row) + fak2 * (decay_col) * (kdotr_col));
103
104 olk(0, 0) = ssol * (std::exp(kexparg)); // s-W-s element
105
106 // Integral p-W-s
107 if (lmax_row > 0) {
108 olk(Cart::x, 0) = PmA(0) * olk(0, 0);
109 olk(Cart::y, 0) = PmA(1) * olk(0, 0);
110 olk(Cart::z, 0) = PmA(2) * olk(0, 0);
111 }
112 //------------------------------------------------------
113
114 // Integrals d - W - s
115 if (lmax_row > 1) {
116 COMPLEX term = (cfak) * (olk(0, 0));
117 olk(Cart::xx, 0) = PmA(0) * olk(Cart::x, 0) + term;
118 olk(Cart::xy, 0) = PmA(0) * olk(Cart::y, 0);
119 olk(Cart::xz, 0) = PmA(0) * olk(Cart::z, 0);
120 olk(Cart::yy, 0) = PmA(1) * olk(Cart::y, 0) + term;
121 olk(Cart::yz, 0) = PmA(1) * olk(Cart::z, 0);
122 olk(Cart::zz, 0) = PmA(2) * olk(Cart::z, 0) + term;
123 }
124 //------------------------------------------------------
125 // Integrals f - W - s
126 if (lmax_row > 2) {
127 olk(Cart::xxx, 0) = PmA(0) * olk(Cart::xx, 0) + cfak2 * olk(Cart::x, 0);
128 olk(Cart::xxy, 0) = PmA(1) * olk(Cart::xx, 0);
129 olk(Cart::xxz, 0) = PmA(2) * olk(Cart::xx, 0);
130 olk(Cart::xyy, 0) = PmA(0) * olk(Cart::yy, 0);
131 olk(Cart::xyz, 0) = PmA(0) * olk(Cart::yz, 0);
132 olk(Cart::xzz, 0) = PmA(0) * olk(Cart::zz, 0);
133 olk(Cart::yyy, 0) = PmA(1) * olk(Cart::yy, 0) + cfak2 * olk(Cart::y, 0);
134 olk(Cart::yyz, 0) = PmA(2) * olk(Cart::yy, 0);
135 olk(Cart::yzz, 0) = PmA(1) * olk(Cart::zz, 0);
136 olk(Cart::zzz, 0) = PmA(2) * olk(Cart::zz, 0) + cfak2 * olk(Cart::z, 0);
137 }
138 //------------------------------------------------------
139 // Integrals g - W - s
140 if (lmax_row > 3) {
141 COMPLEX term_xx = (cfak) * (olk(Cart::xx, 0));
142 COMPLEX term_yy = (cfak) * (olk(Cart::yy, 0));
143 COMPLEX term_zz = (cfak) * (olk(Cart::zz, 0));
144 olk(Cart::xxxx, 0) = PmA(0) * olk(Cart::xxx, 0) + 3.0 * term_xx;
145 olk(Cart::xxxy, 0) = PmA(1) * olk(Cart::xxx, 0);
146 olk(Cart::xxxz, 0) = PmA(2) * olk(Cart::xxx, 0);
147 olk(Cart::xxyy, 0) = PmA(0) * olk(Cart::xyy, 0) + term_yy;
148 olk(Cart::xxyz, 0) = PmA(1) * olk(Cart::xxz, 0);
149 olk(Cart::xxzz, 0) = PmA(0) * olk(Cart::xzz, 0) + term_zz;
150 olk(Cart::xyyy, 0) = PmA(0) * olk(Cart::yyy, 0);
151 olk(Cart::xyyz, 0) = PmA(0) * olk(Cart::yyz, 0);
152 olk(Cart::xyzz, 0) = PmA(0) * olk(Cart::yzz, 0);
153 olk(Cart::xzzz, 0) = PmA(0) * olk(Cart::zzz, 0);
154 olk(Cart::yyyy, 0) = PmA(1) * olk(Cart::yyy, 0) + 3.0 * term_yy;
155 olk(Cart::yyyz, 0) = PmA(2) * olk(Cart::yyy, 0);
156 olk(Cart::yyzz, 0) = PmA(1) * olk(Cart::yzz, 0) + term_zz;
157 olk(Cart::yzzz, 0) = PmA(1) * olk(Cart::zzz, 0);
158 olk(Cart::zzzz, 0) = PmA(2) * olk(Cart::zzz, 0) + 3.0 * term_zz;
159 }
160 //------------------------------------------------------
161 // Integrals h - W - s
162 if (lmax_row > 4) {
163 COMPLEX term_xxx = (cfak) * (olk(Cart::xxx, 0));
164 COMPLEX term_yyy = (cfak) * (olk(Cart::yyy, 0));
165 COMPLEX term_zzz = (cfak) * (olk(Cart::zzz, 0));
166 olk(Cart::xxxxx, 0) = PmA(0) * olk(Cart::xxxx, 0) + 4.0 * term_xxx;
167 olk(Cart::xxxxy, 0) = PmA(1) * olk(Cart::xxxx, 0);
168 olk(Cart::xxxxz, 0) = PmA(2) * olk(Cart::xxxx, 0);
169 olk(Cart::xxxyy, 0) = PmA(1) * olk(Cart::xxxy, 0) + term_xxx;
170 olk(Cart::xxxyz, 0) = PmA(1) * olk(Cart::xxxz, 0);
171 olk(Cart::xxxzz, 0) = PmA(2) * olk(Cart::xxxz, 0) + term_xxx;
172 olk(Cart::xxyyy, 0) = PmA(0) * olk(Cart::xyyy, 0) + term_yyy;
173 olk(Cart::xxyyz, 0) = PmA(2) * olk(Cart::xxyy, 0);
174 olk(Cart::xxyzz, 0) = PmA(1) * olk(Cart::xxzz, 0);
175 olk(Cart::xxzzz, 0) = PmA(0) * olk(Cart::xzzz, 0) + term_zzz;
176 olk(Cart::xyyyy, 0) = PmA(0) * olk(Cart::yyyy, 0);
177 olk(Cart::xyyyz, 0) = PmA(0) * olk(Cart::yyyz, 0);
178 olk(Cart::xyyzz, 0) = PmA(0) * olk(Cart::yyzz, 0);
179 olk(Cart::xyzzz, 0) = PmA(0) * olk(Cart::yzzz, 0);
180 olk(Cart::xzzzz, 0) = PmA(0) * olk(Cart::zzzz, 0);
181 olk(Cart::yyyyy, 0) = PmA(1) * olk(Cart::yyyy, 0) + 4.0 * term_yyy;
182 olk(Cart::yyyyz, 0) = PmA(2) * olk(Cart::yyyy, 0);
183 olk(Cart::yyyzz, 0) = PmA(2) * olk(Cart::yyyz, 0) + term_yyy;
184 olk(Cart::yyzzz, 0) = PmA(1) * olk(Cart::yzzz, 0) + term_zzz;
185 olk(Cart::yzzzz, 0) = PmA(1) * olk(Cart::zzzz, 0);
186 olk(Cart::zzzzz, 0) = PmA(2) * olk(Cart::zzzz, 0) + 4.0 * term_zzz;
187 }
188 //------------------------------------------------------
189 // Integrals i -W - s
190 if (lmax_row > 5) {
191 COMPLEX term_xxxx = (cfak) * (olk(Cart::xxxx, 0));
192 COMPLEX term_xyyy = (cfak) * (olk(Cart::xyyy, 0));
193 COMPLEX term_xzzz = (cfak) * (olk(Cart::xzzz, 0));
194 COMPLEX term_yyyy = (cfak) * (olk(Cart::yyyy, 0));
195 COMPLEX term_yyzz = (cfak) * (olk(Cart::yyzz, 0));
196 COMPLEX term_yzzz = (cfak) * (olk(Cart::yzzz, 0));
197 COMPLEX term_zzzz = (cfak) * (olk(Cart::zzzz, 0));
198 olk(Cart::xxxxxx, 0) = PmA(0) * olk(Cart::xxxxx, 0) + 5.0 * term_xxxx;
199 olk(Cart::xxxxxy, 0) = PmA(1) * olk(Cart::xxxxx, 0);
200 olk(Cart::xxxxxz, 0) = PmA(2) * olk(Cart::xxxxx, 0);
201 olk(Cart::xxxxyy, 0) = PmA(1) * olk(Cart::xxxxy, 0) + term_xxxx;
202 olk(Cart::xxxxyz, 0) = PmA(1) * olk(Cart::xxxxz, 0);
203 olk(Cart::xxxxzz, 0) = PmA(2) * olk(Cart::xxxxz, 0) + term_xxxx;
204 olk(Cart::xxxyyy, 0) = PmA(0) * olk(Cart::xxyyy, 0) + 2.0 * term_xyyy;
205 olk(Cart::xxxyyz, 0) = PmA(2) * olk(Cart::xxxyy, 0);
206 olk(Cart::xxxyzz, 0) = PmA(1) * olk(Cart::xxxzz, 0);
207 olk(Cart::xxxzzz, 0) = PmA(0) * olk(Cart::xxzzz, 0) + 2.0 * term_xzzz;
208 olk(Cart::xxyyyy, 0) = PmA(0) * olk(Cart::xyyyy, 0) + term_yyyy;
209 olk(Cart::xxyyyz, 0) = PmA(2) * olk(Cart::xxyyy, 0);
210 olk(Cart::xxyyzz, 0) = PmA(0) * olk(Cart::xyyzz, 0) + term_yyzz;
211 olk(Cart::xxyzzz, 0) = PmA(1) * olk(Cart::xxzzz, 0);
212 olk(Cart::xxzzzz, 0) = PmA(0) * olk(Cart::xzzzz, 0) + term_zzzz;
213 olk(Cart::xyyyyy, 0) = PmA(0) * olk(Cart::yyyyy, 0);
214 olk(Cart::xyyyyz, 0) = PmA(0) * olk(Cart::yyyyz, 0);
215 olk(Cart::xyyyzz, 0) = PmA(0) * olk(Cart::yyyzz, 0);
216 olk(Cart::xyyzzz, 0) = PmA(0) * olk(Cart::yyzzz, 0);
217 olk(Cart::xyzzzz, 0) = PmA(0) * olk(Cart::yzzzz, 0);
218 olk(Cart::xzzzzz, 0) = PmA(0) * olk(Cart::zzzzz, 0);
219 olk(Cart::yyyyyy, 0) = PmA(1) * olk(Cart::yyyyy, 0) + 5.0 * term_yyyy;
220 olk(Cart::yyyyyz, 0) = PmA(2) * olk(Cart::yyyyy, 0);
221 olk(Cart::yyyyzz, 0) = PmA(2) * olk(Cart::yyyyz, 0) + term_yyyy;
222 olk(Cart::yyyzzz, 0) = PmA(1) * olk(Cart::yyzzz, 0) + 2.0 * term_yzzz;
223 olk(Cart::yyzzzz, 0) = PmA(1) * olk(Cart::yzzzz, 0) + term_zzzz;
224 olk(Cart::yzzzzz, 0) = PmA(1) * olk(Cart::zzzzz, 0);
225 olk(Cart::zzzzzz, 0) = PmA(2) * olk(Cart::zzzzz, 0) + 5.0 * term_zzzz;
226 }
227 //------------------------------------------------------
228 if (lmax_col > 0) {
229
230 // Integrals s - W - p
231 olk(0, Cart::x) = PmB(0) * olk(0, 0);
232 olk(0, Cart::y) = PmB(1) * olk(0, 0);
233 olk(0, Cart::z) = PmB(2) * olk(0, 0);
234 //------------------------------------------------------
235
236 // Integrals p - W - p d - W - p f - W - p g - W - p h -
237 // W - p i - W - p
238 for (Index i = 1; i < n_orbitals[lmax_row]; i++) {
239 // COMPLEX cnx(nx[i] * fak, 0.0);
240 olk(i, Cart::x) =
241 PmB(0) * olk(i, 0) + double(nx[i]) * cfak * olk(i_less_x[i], 0);
242 // COMPLEX cny = (ny[i] * fak, 0.0);
243 olk(i, Cart::y) =
244 PmB(1) * olk(i, 0) + double(ny[i]) * cfak * olk(i_less_y[i], 0);
245 // COMPLEX cnz = (nz[i] * fak, 0.0);
246 olk(i, Cart::z) =
247 PmB(2) * olk(i, 0) + double(nz[i]) * cfak * olk(i_less_z[i], 0);
248 }
249 //------------------------------------------------------
250
251 } // end if (lmax_col > 0)
252 if (lmax_col > 1) {
253
254 // Integrals s - W - d
255 COMPLEX term = cfak * olk(0, 0);
256 olk(0, Cart::xx) = PmB(0) * olk(0, Cart::x) + term;
257 olk(0, Cart::xy) = PmB(0) * olk(0, Cart::y);
258 olk(0, Cart::xz) = PmB(0) * olk(0, Cart::z);
259 olk(0, Cart::yy) = PmB(1) * olk(0, Cart::y) + term;
260 olk(0, Cart::yz) = PmB(1) * olk(0, Cart::z);
261 olk(0, Cart::zz) = PmB(2) * olk(0, Cart::z) + term;
262 //------------------------------------------------------
263
264 // Integrals p - W - d d - W - d f - W - d g - W - d h -
265 // W - d i - W - d
266 for (Index i = 1; i < n_orbitals[lmax_row]; i++) {
267 COMPLEX term_loc = cfak * olk(i, 0);
268 olk(i, Cart::xx) = PmB(0) * olk(i, Cart::x) +
269 double(nx[i]) * cfak * olk(i_less_x[i], Cart::x) +
270 term_loc;
271 olk(i, Cart::xy) = PmB(0) * olk(i, Cart::y) +
272 double(nx[i]) * cfak * olk(i_less_x[i], Cart::y);
273 olk(i, Cart::xz) = PmB(0) * olk(i, Cart::z) +
274 double(nx[i]) * cfak * olk(i_less_x[i], Cart::z);
275 olk(i, Cart::yy) = PmB(1) * olk(i, Cart::y) +
276 double(ny[i]) * cfak * olk(i_less_y[i], Cart::y) +
277 term_loc;
278 olk(i, Cart::yz) = PmB(1) * olk(i, Cart::z) +
279 double(ny[i]) * cfak * olk(i_less_y[i], Cart::z);
280 olk(i, Cart::zz) = PmB(2) * olk(i, Cart::z) +
281 double(nz[i]) * cfak * olk(i_less_z[i], Cart::z) +
282 term_loc;
283 }
284 //------------------------------------------------------
285
286 } // end if (lmax_col > 1)
287
288 if (lmax_col > 2) {
289
290 // Integrals s - W - f
291 olk(0, Cart::xxx) =
292 PmB(0) * olk(0, Cart::xx) + 2.0 * cfak * olk(0, Cart::x);
293 olk(0, Cart::xxy) = PmB(1) * olk(0, Cart::xx);
294 olk(0, Cart::xxz) = PmB(2) * olk(0, Cart::xx);
295 olk(0, Cart::xyy) = PmB(0) * olk(0, Cart::yy);
296 olk(0, Cart::xyz) = PmB(0) * olk(0, Cart::yz);
297 olk(0, Cart::xzz) = PmB(0) * olk(0, Cart::zz);
298 olk(0, Cart::yyy) =
299 PmB(1) * olk(0, Cart::yy) + 2.0 * cfak * olk(0, Cart::y);
300 olk(0, Cart::yyz) = PmB(2) * olk(0, Cart::yy);
301 olk(0, Cart::yzz) = PmB(1) * olk(0, Cart::zz);
302 olk(0, Cart::zzz) =
303 PmB(2) * olk(0, Cart::zz) + 2.0 * cfak * olk(0, Cart::z);
304 //------------------------------------------------------
305
306 // Integrals p - f d - f f - f g - f h - f i - f
307 for (Index i = 1; i < n_orbitals[lmax_row]; i++) {
308 COMPLEX term_x = 2.0 * cfak * olk(i, Cart::x);
309 COMPLEX term_y = 2.0 * cfak * olk(i, Cart::y);
310 COMPLEX term_z = 2.0 * cfak * olk(i, Cart::z);
311 olk(i, Cart::xxx) =
312 PmB(0) * olk(i, Cart::xx) +
313 double(nx[i]) * cfak * olk(i_less_x[i], Cart::xx) + term_x;
314 olk(i, Cart::xxy) = PmB(1) * olk(i, Cart::xx) +
315 double(ny[i]) * cfak * olk(i_less_y[i], Cart::xx);
316 olk(i, Cart::xxz) = PmB(2) * olk(i, Cart::xx) +
317 double(nz[i]) * cfak * olk(i_less_z[i], Cart::xx);
318 olk(i, Cart::xyy) = PmB(0) * olk(i, Cart::yy) +
319 double(nx[i]) * cfak * olk(i_less_x[i], Cart::yy);
320 olk(i, Cart::xyz) = PmB(0) * olk(i, Cart::yz) +
321 double(nx[i]) * cfak * olk(i_less_x[i], Cart::yz);
322 olk(i, Cart::xzz) = PmB(0) * olk(i, Cart::zz) +
323 double(nx[i]) * cfak * olk(i_less_x[i], Cart::zz);
324 olk(i, Cart::yyy) =
325 PmB(1) * olk(i, Cart::yy) +
326 double(ny[i]) * cfak * olk(i_less_y[i], Cart::yy) + term_y;
327 olk(i, Cart::yyz) = PmB(2) * olk(i, Cart::yy) +
328 double(nz[i]) * cfak * olk(i_less_z[i], Cart::yy);
329 olk(i, Cart::yzz) = PmB(1) * olk(i, Cart::zz) +
330 double(ny[i]) * cfak * olk(i_less_y[i], Cart::zz);
331 olk(i, Cart::zzz) =
332 PmB(2) * olk(i, Cart::zz) +
333 double(nz[i]) * cfak * olk(i_less_z[i], Cart::zz) + term_z;
334 }
335 //------------------------------------------------------
336
337 } // end if (lmax_col > 2)
338
339 if (lmax_col > 3) {
340
341 // Integrals s - W - g
342 COMPLEX term_xx = cfak * olk(0, Cart::xx);
343 COMPLEX term_yy = cfak * olk(0, Cart::yy);
344 COMPLEX term_zz = cfak * olk(0, Cart::zz);
345 olk(0, Cart::xxxx) = PmB(0) * olk(0, Cart::xxx) + 3.0 * term_xx;
346 olk(0, Cart::xxxy) = PmB(1) * olk(0, Cart::xxx);
347 olk(0, Cart::xxxz) = PmB(2) * olk(0, Cart::xxx);
348 olk(0, Cart::xxyy) = PmB(0) * olk(0, Cart::xyy) + term_yy;
349 olk(0, Cart::xxyz) = PmB(1) * olk(0, Cart::xxz);
350 olk(0, Cart::xxzz) = PmB(0) * olk(0, Cart::xzz) + term_zz;
351 olk(0, Cart::xyyy) = PmB(0) * olk(0, Cart::yyy);
352 olk(0, Cart::xyyz) = PmB(0) * olk(0, Cart::yyz);
353 olk(0, Cart::xyzz) = PmB(0) * olk(0, Cart::yzz);
354 olk(0, Cart::xzzz) = PmB(0) * olk(0, Cart::zzz);
355 olk(0, Cart::yyyy) = PmB(1) * olk(0, Cart::yyy) + 3.0 * term_yy;
356 olk(0, Cart::yyyz) = PmB(2) * olk(0, Cart::yyy);
357 olk(0, Cart::yyzz) = PmB(1) * olk(0, Cart::yzz) + term_zz;
358 olk(0, Cart::yzzz) = PmB(1) * olk(0, Cart::zzz);
359 olk(0, Cart::zzzz) = PmB(2) * olk(0, Cart::zzz) + 3.0 * term_zz;
360 //------------------------------------------------------
361
362 // Integrals p - g d - g f - g g - g h - g i - g
363 for (Index i = 1; i < n_orbitals[lmax_row]; i++) {
364 COMPLEX term_xx_loc = cfak * olk(i, Cart::xx);
365 COMPLEX term_yy_loc = cfak * olk(i, Cart::yy);
366 COMPLEX term_zz_loc = cfak * olk(i, Cart::zz);
367 olk(i, Cart::xxxx) =
368 PmB(0) * olk(i, Cart::xxx) +
369 double(nx[i]) * cfak * olk(i_less_x[i], Cart::xxx) +
370 3.0 * term_xx_loc;
371 olk(i, Cart::xxxy) =
372 PmB(1) * olk(i, Cart::xxx) +
373 double(ny[i]) * cfak * olk(i_less_y[i], Cart::xxx);
374 olk(i, Cart::xxxz) =
375 PmB(2) * olk(i, Cart::xxx) +
376 double(nz[i]) * cfak * olk(i_less_z[i], Cart::xxx);
377 olk(i, Cart::xxyy) =
378 PmB(0) * olk(i, Cart::xyy) +
379 double(nx[i]) * cfak * olk(i_less_x[i], Cart::xyy) + term_yy_loc;
380 olk(i, Cart::xxyz) =
381 PmB(1) * olk(i, Cart::xxz) +
382 double(ny[i]) * cfak * olk(i_less_y[i], Cart::xxz);
383 olk(i, Cart::xxzz) =
384 PmB(0) * olk(i, Cart::xzz) +
385 double(nx[i]) * cfak * olk(i_less_x[i], Cart::xzz) + term_zz_loc;
386 olk(i, Cart::xyyy) =
387 PmB(0) * olk(i, Cart::yyy) +
388 double(nx[i]) * cfak * olk(i_less_x[i], Cart::yyy);
389 olk(i, Cart::xyyz) =
390 PmB(0) * olk(i, Cart::yyz) +
391 double(nx[i]) * cfak * olk(i_less_x[i], Cart::yyz);
392 olk(i, Cart::xyzz) =
393 PmB(0) * olk(i, Cart::yzz) +
394 double(nx[i]) * cfak * olk(i_less_x[i], Cart::yzz);
395 olk(i, Cart::xzzz) =
396 PmB(0) * olk(i, Cart::zzz) +
397 double(nx[i]) * cfak * olk(i_less_x[i], Cart::zzz);
398 olk(i, Cart::yyyy) =
399 PmB(1) * olk(i, Cart::yyy) +
400 double(ny[i]) * cfak * olk(i_less_y[i], Cart::yyy) +
401 3.0 * term_yy_loc;
402 olk(i, Cart::yyyz) =
403 PmB(2) * olk(i, Cart::yyy) +
404 double(nz[i]) * cfak * olk(i_less_z[i], Cart::yyy);
405 olk(i, Cart::yyzz) =
406 PmB(1) * olk(i, Cart::yzz) +
407 double(ny[i]) * cfak * olk(i_less_y[i], Cart::yzz) + term_zz_loc;
408 olk(i, Cart::yzzz) =
409 PmB(1) * olk(i, Cart::zzz) +
410 double(ny[i]) * cfak * olk(i_less_y[i], Cart::zzz);
411 olk(i, Cart::zzzz) =
412 PmB(2) * olk(i, Cart::zzz) +
413 double(nz[i]) * cfak * olk(i_less_z[i], Cart::zzz) +
414 3.0 * term_zz_loc;
415 }
416 //------------------------------------------------------
417
418 } // end if (lmax_col > 3)
419
420 if (lmax_col > 4) {
421
422 // Integrals s - h
423 COMPLEX term_xxx = cfak * olk(0, Cart::xxx);
424 COMPLEX term_yyy = cfak * olk(0, Cart::yyy);
425 COMPLEX term_zzz = cfak * olk(0, Cart::zzz);
426 olk(0, Cart::xxxxx) = PmB(0) * olk(0, Cart::xxxx) + 4.0 * term_xxx;
427 olk(0, Cart::xxxxy) = PmB(1) * olk(0, Cart::xxxx);
428 olk(0, Cart::xxxxz) = PmB(2) * olk(0, Cart::xxxx);
429 olk(0, Cart::xxxyy) = PmB(1) * olk(0, Cart::xxxy) + term_xxx;
430 olk(0, Cart::xxxyz) = PmB(1) * olk(0, Cart::xxxz);
431 olk(0, Cart::xxxzz) = PmB(2) * olk(0, Cart::xxxz) + term_xxx;
432 olk(0, Cart::xxyyy) = PmB(0) * olk(0, Cart::xyyy) + term_yyy;
433 olk(0, Cart::xxyyz) = PmB(2) * olk(0, Cart::xxyy);
434 olk(0, Cart::xxyzz) = PmB(1) * olk(0, Cart::xxzz);
435 olk(0, Cart::xxzzz) = PmB(0) * olk(0, Cart::xzzz) + term_zzz;
436 olk(0, Cart::xyyyy) = PmB(0) * olk(0, Cart::yyyy);
437 olk(0, Cart::xyyyz) = PmB(0) * olk(0, Cart::yyyz);
438 olk(0, Cart::xyyzz) = PmB(0) * olk(0, Cart::yyzz);
439 olk(0, Cart::xyzzz) = PmB(0) * olk(0, Cart::yzzz);
440 olk(0, Cart::xzzzz) = PmB(0) * olk(0, Cart::zzzz);
441 olk(0, Cart::yyyyy) = PmB(1) * olk(0, Cart::yyyy) + 4.0 * term_yyy;
442 olk(0, Cart::yyyyz) = PmB(2) * olk(0, Cart::yyyy);
443 olk(0, Cart::yyyzz) = PmB(2) * olk(0, Cart::yyyz) + term_yyy;
444 olk(0, Cart::yyzzz) = PmB(1) * olk(0, Cart::yzzz) + term_zzz;
445 olk(0, Cart::yzzzz) = PmB(1) * olk(0, Cart::zzzz);
446 olk(0, Cart::zzzzz) = PmB(2) * olk(0, Cart::zzzz) + 4.0 * term_zzz;
447 //------------------------------------------------------
448
449 // Integrals p - h d - h f - h g - h h - h i - h
450 for (Index i = 1; i < n_orbitals[lmax_row]; i++) {
451 COMPLEX term_xxx_loc = cfak * olk(i, Cart::xxx);
452 COMPLEX term_yyy_loc = cfak * olk(i, Cart::yyy);
453 COMPLEX term_zzz_loc = cfak * olk(i, Cart::zzz);
454 olk(i, Cart::xxxxx) =
455 PmB(0) * olk(i, Cart::xxxx) +
456 double(nx[i]) * cfak * olk(i_less_x[i], Cart::xxxx) +
457 4.0 * term_xxx_loc;
458 olk(i, Cart::xxxxy) =
459 PmB(1) * olk(i, Cart::xxxx) +
460 double(ny[i]) * cfak * olk(i_less_y[i], Cart::xxxx);
461 olk(i, Cart::xxxxz) =
462 PmB(2) * olk(i, Cart::xxxx) +
463 double(nz[i]) * cfak * olk(i_less_z[i], Cart::xxxx);
464 olk(i, Cart::xxxyy) =
465 PmB(1) * olk(i, Cart::xxxy) +
466 double(ny[i]) * cfak * olk(i_less_y[i], Cart::xxxy) +
467 term_xxx_loc;
468 olk(i, Cart::xxxyz) =
469 PmB(1) * olk(i, Cart::xxxz) +
470 double(ny[i]) * cfak * olk(i_less_y[i], Cart::xxxz);
471 olk(i, Cart::xxxzz) =
472 PmB(2) * olk(i, Cart::xxxz) +
473 double(nz[i]) * cfak * olk(i_less_z[i], Cart::xxxz) +
474 term_xxx_loc;
475 olk(i, Cart::xxyyy) =
476 PmB(0) * olk(i, Cart::xyyy) +
477 double(nx[i]) * cfak * olk(i_less_x[i], Cart::xyyy) +
478 term_yyy_loc;
479 olk(i, Cart::xxyyz) =
480 PmB(2) * olk(i, Cart::xxyy) +
481 double(nz[i]) * cfak * olk(i_less_z[i], Cart::xxyy);
482 olk(i, Cart::xxyzz) =
483 PmB(1) * olk(i, Cart::xxzz) +
484 double(ny[i]) * cfak * olk(i_less_y[i], Cart::xxzz);
485 olk(i, Cart::xxzzz) =
486 PmB(0) * olk(i, Cart::xzzz) +
487 double(nx[i]) * cfak * olk(i_less_x[i], Cart::xzzz) +
488 term_zzz_loc;
489 olk(i, Cart::xyyyy) =
490 PmB(0) * olk(i, Cart::yyyy) +
491 double(nx[i]) * cfak * olk(i_less_x[i], Cart::yyyy);
492 olk(i, Cart::xyyyz) =
493 PmB(0) * olk(i, Cart::yyyz) +
494 double(nx[i]) * cfak * olk(i_less_x[i], Cart::yyyz);
495 olk(i, Cart::xyyzz) =
496 PmB(0) * olk(i, Cart::yyzz) +
497 double(nx[i]) * cfak * olk(i_less_x[i], Cart::yyzz);
498 olk(i, Cart::xyzzz) =
499 PmB(0) * olk(i, Cart::yzzz) +
500 double(nx[i]) * cfak * olk(i_less_x[i], Cart::yzzz);
501 olk(i, Cart::xzzzz) =
502 PmB(0) * olk(i, Cart::zzzz) +
503 double(nx[i]) * cfak * olk(i_less_x[i], Cart::zzzz);
504 olk(i, Cart::yyyyy) =
505 PmB(1) * olk(i, Cart::yyyy) +
506 double(ny[i]) * cfak * olk(i_less_y[i], Cart::yyyy) +
507 4.0 * term_yyy_loc;
508 olk(i, Cart::yyyyz) =
509 PmB(2) * olk(i, Cart::yyyy) +
510 double(nz[i]) * cfak * olk(i_less_z[i], Cart::yyyy);
511 olk(i, Cart::yyyzz) =
512 PmB(2) * olk(i, Cart::yyyz) +
513 double(nz[i]) * cfak * olk(i_less_z[i], Cart::yyyz) +
514 term_yyy_loc;
515 olk(i, Cart::yyzzz) =
516 PmB(1) * olk(i, Cart::yzzz) +
517 double(ny[i]) * cfak * olk(i_less_y[i], Cart::yzzz) +
518 term_zzz_loc;
519 olk(i, Cart::yzzzz) =
520 PmB(1) * olk(i, Cart::zzzz) +
521 double(ny[i]) * cfak * olk(i_less_y[i], Cart::zzzz);
522 olk(i, Cart::zzzzz) =
523 PmB(2) * olk(i, Cart::zzzz) +
524 double(nz[i]) * cfak * olk(i_less_z[i], Cart::zzzz) +
525 4.0 * term_zzz_loc;
526 }
527 //------------------------------------------------------
528
529 } // end if (lmax_col > 4)
530
531 if (lmax_col > 5) {
532
533 // Integrals s - W -i
534 COMPLEX term_xxxx = cfak * olk(0, Cart::xxxx);
535 COMPLEX term_xyyy = cfak * olk(0, Cart::xyyy);
536 COMPLEX term_xzzz = cfak * olk(0, Cart::xzzz);
537 COMPLEX term_yyyy = cfak * olk(0, Cart::yyyy);
538 COMPLEX term_yyzz = cfak * olk(0, Cart::yyzz);
539 COMPLEX term_yzzz = cfak * olk(0, Cart::yzzz);
540 COMPLEX term_zzzz = cfak * olk(0, Cart::zzzz);
541 olk(0, Cart::xxxxxx) = PmB(0) * olk(0, Cart::xxxxx) + 5.0 * term_xxxx;
542 olk(0, Cart::xxxxxy) = PmB(1) * olk(0, Cart::xxxxx);
543 olk(0, Cart::xxxxxz) = PmB(2) * olk(0, Cart::xxxxx);
544 olk(0, Cart::xxxxyy) = PmB(1) * olk(0, Cart::xxxxy) + term_xxxx;
545 olk(0, Cart::xxxxyz) = PmB(1) * olk(0, Cart::xxxxz);
546 olk(0, Cart::xxxxzz) = PmB(2) * olk(0, Cart::xxxxz) + term_xxxx;
547 olk(0, Cart::xxxyyy) = PmB(0) * olk(0, Cart::xxyyy) + 2.0 * term_xyyy;
548 olk(0, Cart::xxxyyz) = PmB(2) * olk(0, Cart::xxxyy);
549 olk(0, Cart::xxxyzz) = PmB(1) * olk(0, Cart::xxxzz);
550 olk(0, Cart::xxxzzz) = PmB(0) * olk(0, Cart::xxzzz) + 2.0 * term_xzzz;
551 olk(0, Cart::xxyyyy) = PmB(0) * olk(0, Cart::xyyyy) + term_yyyy;
552 olk(0, Cart::xxyyyz) = PmB(2) * olk(0, Cart::xxyyy);
553 olk(0, Cart::xxyyzz) = PmB(0) * olk(0, Cart::xyyzz) + term_yyzz;
554 olk(0, Cart::xxyzzz) = PmB(1) * olk(0, Cart::xxzzz);
555 olk(0, Cart::xxzzzz) = PmB(0) * olk(0, Cart::xzzzz) + term_zzzz;
556 olk(0, Cart::xyyyyy) = PmB(0) * olk(0, Cart::yyyyy);
557 olk(0, Cart::xyyyyz) = PmB(0) * olk(0, Cart::yyyyz);
558 olk(0, Cart::xyyyzz) = PmB(0) * olk(0, Cart::yyyzz);
559 olk(0, Cart::xyyzzz) = PmB(0) * olk(0, Cart::yyzzz);
560 olk(0, Cart::xyzzzz) = PmB(0) * olk(0, Cart::yzzzz);
561 olk(0, Cart::xzzzzz) = PmB(0) * olk(0, Cart::zzzzz);
562 olk(0, Cart::yyyyyy) = PmB(1) * olk(0, Cart::yyyyy) + 5.0 * term_yyyy;
563 olk(0, Cart::yyyyyz) = PmB(2) * olk(0, Cart::yyyyy);
564 olk(0, Cart::yyyyzz) = PmB(2) * olk(0, Cart::yyyyz) + term_yyyy;
565 olk(0, Cart::yyyzzz) = PmB(1) * olk(0, Cart::yyzzz) + 2.0 * term_yzzz;
566 olk(0, Cart::yyzzzz) = PmB(1) * olk(0, Cart::yzzzz) + term_zzzz;
567 olk(0, Cart::yzzzzz) = PmB(1) * olk(0, Cart::zzzzz);
568 olk(0, Cart::zzzzzz) = PmB(2) * olk(0, Cart::zzzzz) + 5.0 * term_zzzz;
569 //------------------------------------------------------
570
571 // Integrals p - W - i d - W - i f - W - i g - W -i h -
572 // W - i i - W - i
573 for (Index i = 1; i < n_orbitals[lmax_row]; i++) {
574 COMPLEX term_xxxx_loc = cfak * olk(i, Cart::xxxx);
575 COMPLEX term_xyyy_loc = cfak * olk(i, Cart::xyyy);
576 COMPLEX term_xzzz_loc = cfak * olk(i, Cart::xzzz);
577 COMPLEX term_yyyy_loc = cfak * olk(i, Cart::yyyy);
578 COMPLEX term_yyzz_loc = cfak * olk(i, Cart::yyzz);
579 COMPLEX term_yzzz_loc = cfak * olk(i, Cart::yzzz);
580 COMPLEX term_zzzz_loc = cfak * olk(i, Cart::zzzz);
581 olk(i, Cart::xxxxxx) =
582 PmB(0) * olk(i, Cart::xxxxx) +
583 double(nx[i]) * cfak * olk(i_less_x[i], Cart::xxxxx) +
584 5.0 * term_xxxx_loc;
585 olk(i, Cart::xxxxxy) =
586 PmB(1) * olk(i, Cart::xxxxx) +
587 double(ny[i]) * cfak * olk(i_less_y[i], Cart::xxxxx);
588 olk(i, Cart::xxxxxz) =
589 PmB(2) * olk(i, Cart::xxxxx) +
590 double(nz[i]) * cfak * olk(i_less_z[i], Cart::xxxxx);
591 olk(i, Cart::xxxxyy) =
592 PmB(1) * olk(i, Cart::xxxxy) +
593 double(ny[i]) * cfak * olk(i_less_y[i], Cart::xxxxy) +
594 term_xxxx_loc;
595 olk(i, Cart::xxxxyz) =
596 PmB(1) * olk(i, Cart::xxxxz) +
597 double(ny[i]) * cfak * olk(i_less_y[i], Cart::xxxxz);
598 olk(i, Cart::xxxxzz) =
599 PmB(2) * olk(i, Cart::xxxxz) +
600 double(nz[i]) * cfak * olk(i_less_z[i], Cart::xxxxz) +
601 term_xxxx_loc;
602 olk(i, Cart::xxxyyy) =
603 PmB(0) * olk(i, Cart::xxyyy) +
604 double(nx[i]) * cfak * olk(i_less_x[i], Cart::xxyyy) +
605 2.0 * term_xyyy_loc;
606 olk(i, Cart::xxxyyz) =
607 PmB(2) * olk(i, Cart::xxxyy) +
608 double(nz[i]) * cfak * olk(i_less_z[i], Cart::xxxyy);
609 olk(i, Cart::xxxyzz) =
610 PmB(1) * olk(i, Cart::xxxzz) +
611 double(ny[i]) * cfak * olk(i_less_y[i], Cart::xxxzz);
612 olk(i, Cart::xxxzzz) =
613 PmB(0) * olk(i, Cart::xxzzz) +
614 double(nx[i]) * cfak * olk(i_less_x[i], Cart::xxzzz) +
615 2.0 * term_xzzz_loc;
616 olk(i, Cart::xxyyyy) =
617 PmB(0) * olk(i, Cart::xyyyy) +
618 double(nx[i]) * cfak * olk(i_less_x[i], Cart::xyyyy) +
619 term_yyyy_loc;
620 olk(i, Cart::xxyyyz) =
621 PmB(2) * olk(i, Cart::xxyyy) +
622 double(nz[i]) * cfak * olk(i_less_z[i], Cart::xxyyy);
623 olk(i, Cart::xxyyzz) =
624 PmB(0) * olk(i, Cart::xyyzz) +
625 double(nx[i]) * cfak * olk(i_less_x[i], Cart::xyyzz) +
626 term_yyzz_loc;
627 olk(i, Cart::xxyzzz) =
628 PmB(1) * olk(i, Cart::xxzzz) +
629 double(ny[i]) * cfak * olk(i_less_y[i], Cart::xxzzz);
630 olk(i, Cart::xxzzzz) =
631 PmB(0) * olk(i, Cart::xzzzz) +
632 double(nx[i]) * cfak * olk(i_less_x[i], Cart::xzzzz) +
633 term_zzzz_loc;
634 olk(i, Cart::xyyyyy) =
635 PmB(0) * olk(i, Cart::yyyyy) +
636 double(nx[i]) * cfak * olk(i_less_x[i], Cart::yyyyy);
637 olk(i, Cart::xyyyyz) =
638 PmB(0) * olk(i, Cart::yyyyz) +
639 double(nx[i]) * cfak * olk(i_less_x[i], Cart::yyyyz);
640 olk(i, Cart::xyyyzz) =
641 PmB(0) * olk(i, Cart::yyyzz) +
642 double(nx[i]) * cfak * olk(i_less_x[i], Cart::yyyzz);
643 olk(i, Cart::xyyzzz) =
644 PmB(0) * olk(i, Cart::yyzzz) +
645 double(nx[i]) * cfak * olk(i_less_x[i], Cart::yyzzz);
646 olk(i, Cart::xyzzzz) =
647 PmB(0) * olk(i, Cart::yzzzz) +
648 double(nx[i]) * cfak * olk(i_less_x[i], Cart::yzzzz);
649 olk(i, Cart::xzzzzz) =
650 PmB(0) * olk(i, Cart::zzzzz) +
651 double(nx[i]) * cfak * olk(i_less_x[i], Cart::zzzzz);
652 olk(i, Cart::yyyyyy) =
653 PmB(1) * olk(i, Cart::yyyyy) +
654 double(ny[i]) * cfak * olk(i_less_y[i], Cart::yyyyy) +
655 5.0 * term_yyyy_loc;
656 olk(i, Cart::yyyyyz) =
657 PmB(2) * olk(i, Cart::yyyyy) +
658 double(nz[i]) * cfak * olk(i_less_z[i], Cart::yyyyy);
659 olk(i, Cart::yyyyzz) =
660 PmB(2) * olk(i, Cart::yyyyz) +
661 double(nz[i]) * cfak * olk(i_less_z[i], Cart::yyyyz) +
662 term_yyyy_loc;
663 olk(i, Cart::yyyzzz) =
664 PmB(1) * olk(i, Cart::yyzzz) +
665 double(ny[i]) * cfak * olk(i_less_y[i], Cart::yyzzz) +
666 2.0 * term_yzzz_loc;
667 olk(i, Cart::yyzzzz) =
668 PmB(1) * olk(i, Cart::yzzzz) +
669 double(ny[i]) * cfak * olk(i_less_y[i], Cart::yzzzz) +
670 term_zzzz_loc;
671 olk(i, Cart::yzzzzz) =
672 PmB(1) * olk(i, Cart::zzzzz) +
673 double(ny[i]) * cfak * olk(i_less_y[i], Cart::zzzzz);
674 olk(i, Cart::zzzzzz) =
675 PmB(2) * olk(i, Cart::zzzzz) +
676 double(nz[i]) * cfak * olk(i_less_z[i], Cart::zzzzz) +
677 5.0 * term_zzzz_loc;
678 }
679 //------------------------------------------------------
680
681 } // end if (lmax_col > 5)
682
683 // save to matrix
684 cartesian += AOTransform::getNorm(shell_row.getL(), gaussian_row) *
685 AOTransform::getNorm(shell_col.getL(), gaussian_col) *
686 olk.bottomRightCorner(shell_row.getCartesianNumFunc(),
687 shell_col.getCartesianNumFunc());
688
689 } // close Gaussian shell_col
690
691 } // close Gaussian shell_row
692
693 matrix = AOTransform::tform(shell_row.getL(), shell_col.getL(), cartesian);
694
695} // End AOPlanewave
696
698 const std::vector<Eigen::Vector3d>& kpoints) {
699
701 Eigen::MatrixXcd::Zero(aobasis.AOBasisSize(), aobasis.AOBasisSize());
702
703 for (const auto& kpoint : kpoints) {
704 setkVector(kpoint);
705 aopotential_ += Fill(aobasis);
706 }
707
708 return;
709}
710
711} // namespace xtp
712} // namespace votca
Container to hold Basisfunctions for all atoms.
Definition aobasis.h:42
Index AOBasisSize() const
Definition aobasis.h:46
void FillPotential(const AOBasis &aobasis, const std::vector< Eigen::Vector3d > &kpoints)
void FillBlock(Eigen::Block< Eigen::MatrixXcd > &matrix, const AOShell &shell_row, const AOShell &shell_col) const override
void setkVector(const Eigen::Vector3d &k)
Definition aopotential.h:94
Eigen::Vector3d k_
Definition aopotential.h:95
Eigen::Matrix< std::complex< double >, Eigen::Dynamic, Eigen::Dynamic > Fill(const AOBasis &aobasis) const
Eigen::Matrix< std::complex< double >, Eigen::Dynamic, Eigen::Dynamic > aopotential_
Definition aopotential.h:49
Index getCartesianNumFunc() const
Definition aoshell.h:104
const Eigen::Vector3d & getPos() const
Definition aoshell.h:113
static std::array< int, 165 > i_less_y()
static std::array< int, 165 > i_less_z()
static Index getBlockSize(Index lmax)
static Matrix tform(L l_row, L l_col, const Matrix &cartesian)
transforms a cartesian shell to a spherical cartesian shell
static std::array< int, 165 > i_less_x()
static double getNorm(L l, const AOGaussianPrimitive &gaussian)
static std::array< int, 165 > ny()
static std::array< int, 165 > nz()
static std::array< int, 9 > n_orbitals()
static std::array< int, 165 > nx()
base class for all analysis tools
Definition basebead.h:33
Eigen::Index Index
Definition types.h:26