votca 2024.2-dev
Loading...
Searching...
No Matches
libint2_calls.cc
Go to the documentation of this file.
1/*
2 * Copyright 2009-2024 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/ERIs.h"
22#include "votca/xtp/aobasis.h"
23#include "votca/xtp/aomatrix.h"
26
27// include libint last otherwise it overrides eigen
29#define LIBINT2_CONSTEXPR_STATICS 0
30#if defined(__clang__)
31#pragma clang diagnostic push
32#pragma clang diagnostic ignored "-W#warnings"
33#elif defined(__GNUC__)
34#pragma GCC diagnostic push
35#pragma GCC diagnostic ignored "-Warray-bounds"
36#pragma GCC diagnostic ignored "-Wcpp"
37#endif
38#include <libint2.hpp>
39#include <libint2/statics_definition.h>
40#if defined(__clang__)
41#pragma clang diagnostic pop
42#elif defined(__GNUC__)
43#pragma GCC diagnostic pop
44#endif
45
46namespace votca {
47namespace xtp {
48
49std::vector<std::vector<Index>> AOBasis::ComputeShellPairs(
50 double threshold) const {
51
52 Index nthreads = OPENMP::getMaxThreads();
53
54 std::vector<libint2::Shell> shells = GenerateLibintBasis();
55
56 // construct the 2-electron repulsion integrals engine
57 std::vector<libint2::Engine> engines;
58 engines.reserve(nthreads);
59 engines.emplace_back(libint2::Operator::overlap, getMaxNprim(), getMaxL(), 0);
60 for (Index i = 1; i != nthreads; ++i) {
61 engines.push_back(engines[0]);
62 }
63
64 std::vector<std::vector<Index>> pairs(shells.size());
65
66#pragma omp parallel for schedule(dynamic)
67 for (Index s1 = 0; s1 < Index(shells.size()); ++s1) {
68 Index thread_id = OPENMP::getThreadId();
69
70 libint2::Engine& engine = engines[thread_id];
71 const libint2::Engine::target_ptr_vec& buf = engine.results();
72 Index n1 = shells[s1].size();
73
74 for (Index s2 = 0; s2 <= s1; ++s2) {
75 bool on_same_center = (shells[s1].O == shells[s2].O);
76 bool significant = on_same_center;
77 if (!on_same_center) {
78 Index n2 = shells[s2].size();
79 engine.compute(shells[s1], shells[s2]);
80 Eigen::Map<const Eigen::MatrixXd> buf_mat(buf[0], n1, n2);
81 significant = (buf_mat.norm() >= threshold);
82 }
83 if (significant) {
84 pairs[s1].push_back(s2);
85 }
86 }
87 std::sort(pairs[s1].begin(), pairs[s1].end());
88 }
89 return pairs;
90}
91
92using MatrixLibInt =
93 Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
94
95template <libint2::Operator obtype,
96 typename OperatorParams =
97 typename libint2::operator_traits<obtype>::oper_params_type>
98std::array<MatrixLibInt, libint2::operator_traits<obtype>::nopers>
100 OperatorParams oparams = OperatorParams()) {
101
102 Index nthreads = OPENMP::getMaxThreads();
103 std::vector<libint2::Shell> shells = aobasis.GenerateLibintBasis();
104
105 std::vector<std::vector<Index>> shellpair_list = aobasis.ComputeShellPairs();
106
107 Index nopers = static_cast<Index>(libint2::operator_traits<obtype>::nopers);
108 std::array<MatrixLibInt, libint2::operator_traits<obtype>::nopers> result;
109 for (MatrixLibInt& r : result) {
110 r = MatrixLibInt::Zero(aobasis.AOBasisSize(), aobasis.AOBasisSize());
111 }
112
113 std::vector<libint2::Engine> engines(nthreads);
114 engines[0] = libint2::Engine(obtype, aobasis.getMaxNprim(),
115 static_cast<int>(aobasis.getMaxL()), 0);
116 engines[0].set_params(oparams);
117 for (Index i = 1; i < nthreads; ++i) {
118 engines[i] = engines[0];
119 }
120
121 std::vector<Index> shell2bf = aobasis.getMapToBasisFunctions();
122
123#pragma omp parallel for schedule(dynamic)
124 for (Index s1 = 0; s1 < aobasis.getNumofShells(); ++s1) {
125 Index thread_id = OPENMP::getThreadId();
126 libint2::Engine& engine = engines[thread_id];
127 const libint2::Engine::target_ptr_vec& buf = engine.results();
128
129 Index bf1 = shell2bf[s1];
130 Index n1 = shells[s1].size();
131
132 for (Index s2 : shellpair_list[s1]) {
133
134 engine.compute(shells[s1], shells[s2]);
135 if (buf[0] == nullptr) {
136 continue; // if all integrals screened out, skip to next shell set
137 }
138 Index bf2 = shell2bf[s2];
139 Index n2 = shells[s2].size();
140 for (unsigned int op = 0; op != nopers; ++op) {
141 Eigen::Map<const MatrixLibInt> buf_mat(buf[op], n1, n2);
142 result[op].block(bf1, bf2, n1, n2) = buf_mat;
143 if (s1 != s2) { // if s1 >= s2, copy {s1,s2} to the corresponding
144 // {s2,s1} block, note the transpose!
145 result[op].block(bf2, bf1, n2, n1) = buf_mat.transpose();
146 }
147 }
148 }
149 }
150 return result;
151}
152
153/***********************************
154 * KINETIC
155 ***********************************/
159
160/***********************************
161 * OVERLAP
162 ***********************************/
166
167Eigen::MatrixXd AOOverlap::singleShellOverlap(const AOShell& shell) const {
168 libint2::Shell::do_enforce_unit_normalization(false);
169 libint2::Operator obtype = libint2::Operator::overlap;
170 try {
171 libint2::Engine engine(obtype, shell.getSize(),
172 static_cast<int>(shell.getL()), 0);
173
174 const libint2::Engine::target_ptr_vec& buf = engine.results();
175
176 libint2::Shell s = shell.LibintShell();
177 engine.compute(s, s);
178
179 Eigen::Map<const MatrixLibInt> buf_mat(buf[0], shell.getNumFunc(),
180 shell.getNumFunc());
181 return buf_mat;
182 } catch (const libint2::Engine::lmax_exceeded& error) {
183 std::ostringstream oss;
184 oss << "\nA libint error occured:\n"
185 << error.what() << "\n"
186 << "\nYou requested a computation for a shell with angular momentum "
187 << error.lmax_requested()
188 << ",\nbut your libint package only supports angular momenta upto "
189 << error.lmax_limit() - 1 << ".\n"
190 << "\nTo fix this error you will need to reinstall libint with "
191 "support\n"
192 "for higher angular momenta. If you installed your own libint it\n"
193 "should be reconfigured and installed with the option "
194 "--with-max-am=<maxAngularMomentum>.\n"
195 "If you installed libint with your OS package manager, you will\n"
196 "need to setup you own libint installation with the \n"
197 "--with-max-am=<maxAngularMomentum> option set."
198 << std::endl;
199 throw std::runtime_error(oss.str());
200 }
201}
202
203/***********************************
204 * DIPOLE
205 ***********************************/
206void AODipole::Fill(const AOBasis& aobasis) {
207 auto results = computeOneBodyIntegrals<libint2::Operator::emultipole1,
208 std::array<libint2::Shell::real_t, 3>>(
209 aobasis, r_);
210
211 for (Index i = 0; i < 3; i++) {
212 aomatrix_[i] = results[1 + i]; // emultipole1 returns: overlap, x-dipole,
213 // y-dipole, z-dipole
214 }
215}
216
217/***********************************
218 * COULOMB
219 ***********************************/
220void AOCoulomb::Fill(const AOBasis& aobasis) {
222}
223
225 Index nthreads = OPENMP::getMaxThreads();
226 std::vector<libint2::Shell> shells = aobasis.GenerateLibintBasis();
227 std::vector<Index> shell2bf = aobasis.getMapToBasisFunctions();
228
229 aomatrix_ =
230 Eigen::MatrixXd::Zero(aobasis.AOBasisSize(), aobasis.AOBasisSize());
231
232 // build engines for each thread
233 std::vector<libint2::Engine> engines(nthreads);
234 engines[0] =
235 libint2::Engine(libint2::Operator::coulomb, aobasis.getMaxNprim(),
236 static_cast<int>(aobasis.getMaxL()), 0);
237 engines[0].set(libint2::BraKet::xs_xs);
238 for (Index i = 1; i != nthreads; ++i) {
239 engines[i] = engines[0];
240 }
241
242#pragma omp parallel for schedule(dynamic)
243 for (Index s1 = 0; s1 < aobasis.getNumofShells(); ++s1) {
244 libint2::Engine& engine = engines[OPENMP::getThreadId()];
245 const libint2::Engine::target_ptr_vec& buf = engine.results();
246
247 Index bf1 = shell2bf[s1];
248 Index n1 = shells[s1].size();
249 // cannot use shellpairs because this is still a two-center integral and
250 // overlap screening would give wrong result
251 for (Index s2 = 0; s2 <= s1; ++s2) {
252
253 engine.compute2<libint2::Operator::coulomb, libint2::BraKet::xs_xs, 0>(
254 shells[s1], libint2::Shell::unit(), shells[s2],
255 libint2::Shell::unit());
256
257 if (buf[0] == nullptr) {
258 continue; // if all integrals screened out, skip to next shell set
259 }
260 Index bf2 = shell2bf[s2];
261 Index n2 = shells[s2].size();
262
263 Eigen::Map<const MatrixLibInt> buf_mat(buf[0], n1, n2);
264 aomatrix_.block(bf1, bf2, n1, n2) = buf_mat;
265 if (s1 != s2) { // if s1 >= s2, copy {s1,s2} to the corresponding
266 // {s2,s1} block, note the transpose!
267 aomatrix_.block(bf2, bf1, n2, n1) = buf_mat.transpose();
268 }
269 }
270 }
271}
272
273Eigen::MatrixXd ERIs::ComputeSchwarzShells(const AOBasis& basis) const {
274
275 Index noshells = basis.getNumofShells();
276
277 Eigen::MatrixXd result = Eigen::MatrixXd::Zero(noshells, noshells);
278 Index nthreads = OPENMP::getMaxThreads();
279 std::vector<libint2::Engine> engines(nthreads);
280 double epsilon = 0.0;
281 engines[0] = libint2::Engine(libint2::Operator::coulomb, basis.getMaxNprim(),
282 static_cast<int>(basis.getMaxL()), 0, epsilon);
283
284 for (Index i = 1; i < nthreads; ++i) {
285 engines[i] = engines[0];
286 }
287
288 std::vector<libint2::Shell> shells = basis.GenerateLibintBasis();
289
290#pragma omp parallel for schedule(dynamic)
291 for (Index s1 = 0l; s1 < basis.getNumofShells(); ++s1) {
292 Index thread_id = OPENMP::getThreadId();
293 libint2::Engine& engine = engines[thread_id];
294 const libint2::Engine::target_ptr_vec& buf = engine.results();
295 Index n1 = shells[s1].size();
296
297 for (Index s2 = 0l; s2 <= s1; ++s2) {
298 Index n2 = shells[s2].size();
299 Index n12 = n1 * n2;
300
301 engines[thread_id]
302 .compute2<libint2::Operator::coulomb, libint2::BraKet::xx_xx, 0>(
303 shells[s1], shells[s2], shells[s1], shells[s2]);
304
305 Eigen::Map<const MatrixLibInt> buf_mat(buf[0], n12, n12);
306
307 result(s2, s1) = std::sqrt(buf_mat.cwiseAbs().maxCoeff());
308 }
309 }
310 return result.selfadjointView<Eigen::Upper>();
311}
312
313template <bool with_exchange>
314std::array<Eigen::MatrixXd, 2> ERIs::Compute4c(const Eigen::MatrixXd& dmat,
315 double error) const {
316 assert(schwarzscreen_.rows() > 0 && schwarzscreen_.cols() > 0 &&
317 "Please call Initialize_4c before running this");
318 Index nthreads = OPENMP::getMaxThreads();
319
320 Eigen::MatrixXd hartree = Eigen::MatrixXd::Zero(dmat.rows(), dmat.cols());
321 Eigen::MatrixXd exchange;
322 if (with_exchange) {
323 exchange = Eigen::MatrixXd::Zero(dmat.rows(), dmat.cols());
324 }
325 Eigen::MatrixXd dnorm_block = ComputeShellBlockNorm(dmat);
326 double fock_precision = error;
327 // engine precision controls primitive truncation, assume worst-case scenario
328 // (all primitive combinations add up constructively)
329 Index max_nprim4 = maxnprim_ * maxnprim_ * maxnprim_ * maxnprim_;
330 double engine_precision = std::min(fock_precision / dnorm_block.maxCoeff(),
331 std::numeric_limits<double>::epsilon()) /
332 double(max_nprim4);
333 std::vector<libint2::Engine> engines(nthreads);
334 engines[0] = libint2::Engine(libint2::Operator::coulomb, int(maxnprim_),
335 int(maxL_), 0);
336 engines[0].set_precision(engine_precision); // shellset-dependent precision
337 // control will likely break
338 // positive definiteness
339 // stick with this simple recipe
340 for (Index i = 1; i < nthreads; ++i) {
341 engines[i] = engines[0];
342 }
343 Index nshells = basis_.size();
344
345#pragma omp parallel for schedule(dynamic) reduction(+ : hartree) \
346 reduction(+ : exchange)
347 for (Index s1 = 0; s1 < nshells; ++s1) {
348 Index thread_id = OPENMP::getThreadId();
349 libint2::Engine& engine = engines[thread_id];
350 const auto& buf = engine.results();
351 Index start_1 = starts_[s1];
352 const libint2::Shell& shell1 = basis_[s1];
353 Index n1 = shell1.size();
354
355 auto sp12_iter = shellpairdata_[s1].begin();
356 for (Index s2 : shellpairs_[s1]) {
357 Index start_2 = starts_[s2];
358 const libint2::Shell& shell2 = basis_[s2];
359 Index n2 = shell2.size();
360 double dnorm_12 = dnorm_block(s1, s2);
361 const libint2::ShellPair* sp12 = &(*sp12_iter);
362 ++sp12_iter;
363
364 for (Index s3 = 0; s3 <= s1; ++s3) {
365
366 Index start_3 = starts_[s3];
367 const libint2::Shell& shell3 = basis_[s3];
368 Index n3 = shell3.size();
369 auto sp34_iter = shellpairdata_[s3].begin();
370 double dnorm_123 = std::max(dnorm_block(s1, s3),
371 std::max(dnorm_block(s2, s3), dnorm_12));
372 Index s4max = (s1 == s3) ? s2 : s3;
373 for (Index s4 : shellpairs_[s3]) {
374 if (s4 > s4max) {
375 break;
376 } // for each s3, s4 are stored in monotonically increasing
377 // order
378
379 const libint2::ShellPair* sp34 = &(*sp34_iter);
380 // must update the iter even if going to skip s4
381 ++sp34_iter;
382 double dnorm_1234 =
383 std::max(dnorm_block(s1, s4),
384 std::max(dnorm_block(s2, s4),
385 std::max(dnorm_block(s3, s4), dnorm_123)));
386
387 if (dnorm_1234 * schwarzscreen_(s1, s2) * schwarzscreen_(s3, s4) <
388 fock_precision) {
389 continue;
390 }
391
392 const libint2::Shell& shell4 = basis_[s4];
393 engine
394 .compute2<libint2::Operator::coulomb, libint2::BraKet::xx_xx, 0>(
395 shell1, shell2, shell3, shell4, sp12, sp34);
396 const auto* buf_1234 = buf[0];
397 if (buf_1234 == nullptr) {
398 continue; // if all integrals screened out, skip to next quartet
399 }
400 Index start_4 = starts_[s4];
401 Index n4 = shell4.size();
402 Index s12_deg = (s1 == s2) ? 1 : 2;
403 Index s34_deg = (s3 == s4) ? 1 : 2;
404 Index s12_34_deg = (s1 == s3) ? (s2 == s4 ? 1 : 2) : 2;
405 Index s1234_deg = s12_deg * s34_deg * s12_34_deg;
406
407 for (Index f1 = 0, f1234 = 0; f1 != n1; ++f1) {
408 const Index bf1 = f1 + start_1;
409 for (Index f2 = 0; f2 != n2; ++f2) {
410 const Index bf2 = f2 + start_2;
411 for (Index f3 = 0; f3 != n3; ++f3) {
412 const Index bf3 = f3 + start_3;
413 for (Index f4 = 0; f4 != n4; ++f4, ++f1234) {
414 const Index bf4 = f4 + start_4;
415
416 const double value = buf_1234[f1234];
417
418 const double value_scal_by_deg = value * double(s1234_deg);
419
420 hartree(bf1, bf2) += dmat(bf3, bf4) * value_scal_by_deg;
421 hartree(bf3, bf4) += dmat(bf1, bf2) * value_scal_by_deg;
422 if (with_exchange) {
423 exchange(bf1, bf3) -= dmat(bf2, bf4) * value_scal_by_deg;
424 exchange(bf2, bf3) -= dmat(bf1, bf4) * value_scal_by_deg;
425 exchange(bf2, bf4) -= dmat(bf1, bf3) * value_scal_by_deg;
426 exchange(bf1, bf4) -= dmat(bf2, bf3) * value_scal_by_deg;
427 }
428 }
429 }
430 }
431 }
432 }
433 }
434 }
435 }
436 std::array<Eigen::MatrixXd, 2> result2;
437 // 0.25=0.5(symmetrisation)*0.5(our dmat has a factor 2)
438 result2[0] = 0.25 * (hartree + hartree.transpose());
439 if (with_exchange) {
440 // prefactor
441 result2[1] = 0.125 * (exchange + exchange.transpose());
442 }
443 return result2;
444}
445
446template std::array<Eigen::MatrixXd, 2> ERIs::Compute4c<true>(
447 const Eigen::MatrixXd& dmat, double error) const;
448template std::array<Eigen::MatrixXd, 2> ERIs::Compute4c<false>(
449 const Eigen::MatrixXd& dmat, double error) const;
450
451void TCMatrix_dft::Fill(const AOBasis& auxbasis, const AOBasis& dftbasis) {
452 {
453 AOCoulomb auxAOcoulomb;
454 auxAOcoulomb.Fill(auxbasis);
455 inv_sqrt_ = auxAOcoulomb.Pseudo_InvSqrt(1e-8);
456 removedfunctions_ = auxAOcoulomb.Removedfunctions();
457 }
458 matrix_ = std::vector<Symmetric_Matrix>(auxbasis.AOBasisSize());
459
460#pragma omp parallel for schedule(dynamic, 4)
461 for (Index i = 0; i < auxbasis.AOBasisSize(); i++) {
462 matrix_[i] = Symmetric_Matrix(dftbasis.AOBasisSize());
463 }
464
465 Index nthreads = OPENMP::getMaxThreads();
466 std::vector<libint2::Shell> dftshells = dftbasis.GenerateLibintBasis();
467 std::vector<libint2::Shell> auxshells = auxbasis.GenerateLibintBasis();
468 std::vector<libint2::Engine> engines(nthreads);
469 engines[0] = libint2::Engine(
470 libint2::Operator::coulomb,
471 std::max(dftbasis.getMaxNprim(), auxbasis.getMaxNprim()),
472 static_cast<int>(std::max(dftbasis.getMaxL(), auxbasis.getMaxL())), 0);
473 engines[0].set(libint2::BraKet::xs_xx);
474 for (Index i = 1; i < nthreads; ++i) {
475 engines[i] = engines[0];
476 }
477
478 std::vector<Index> shell2bf = dftbasis.getMapToBasisFunctions();
479 std::vector<Index> auxshell2bf = auxbasis.getMapToBasisFunctions();
480
481#pragma omp parallel for schedule(dynamic)
482 for (Index is = dftbasis.getNumofShells() - 1; is >= 0; is--) {
483
484 libint2::Engine& engine = engines[OPENMP::getThreadId()];
485 const libint2::Engine::target_ptr_vec& buf = engine.results();
486 const libint2::Shell& dftshell = dftshells[is];
487 Index start = shell2bf[is];
488 std::vector<Eigen::MatrixXd> block(dftshell.size());
489 for (Index i = 0; i < Index(dftshell.size()); i++) {
490 Index size = start + i + 1;
491 block[i] = Eigen::MatrixXd::Zero(auxbasis.AOBasisSize(), size);
492 }
493
494 for (Index aux = 0; aux < auxbasis.getNumofShells(); aux++) {
495 const libint2::Shell& auxshell = auxshells[aux];
496 Index aux_start = auxshell2bf[aux];
497
498 for (Index dis = 0; dis <= is; dis++) {
499
500 const libint2::Shell& shell_col = dftshells[dis];
501 Index col_start = shell2bf[dis];
502 engine.compute2<libint2::Operator::coulomb, libint2::BraKet::xs_xx, 0>(
503 auxshell, libint2::Shell::unit(), dftshell, shell_col);
504
505 if (buf[0] == nullptr) {
506 continue;
507 }
508 Eigen::TensorMap<Eigen::Tensor<const double, 3, Eigen::RowMajor> const>
509 result(buf[0], auxshell.size(), dftshell.size(), shell_col.size());
510
511 for (size_t left = 0; left < dftshell.size(); left++) {
512 for (size_t auxf = 0; auxf < auxshell.size(); auxf++) {
513 for (size_t col = 0; col < shell_col.size(); col++) {
514 // symmetry
515 if ((col_start + col) > (start + left)) {
516 break;
517 }
518 block[left](aux_start + auxf, col_start + col) =
519 result(auxf, left, col);
520 }
521 }
522 }
523 }
524 }
525
526 for (Index i = 0; i < Index(block.size()); ++i) {
527 Eigen::MatrixXd temp = inv_sqrt_ * block[i];
528 for (Index mu = 0; mu < temp.rows(); ++mu) {
529 for (Index j = 0; j < temp.cols(); ++j) {
530 matrix_[mu](i + start, j) = temp(mu, j);
531 }
532 }
533 }
534 }
535
536 return;
537}
538
539/*
540 * Determines the 3-center integrals for a given shell in the aux basis
541 * by calculating the 3-center repulsion integral of the functions in the
542 * aux shell with ALL functions in the DFT basis set
543 */
544std::vector<Eigen::MatrixXd> ComputeAO3cBlock(const libint2::Shell& auxshell,
545 const AOBasis& dftbasis,
546 libint2::Engine& engine) {
547 std::vector<Eigen::MatrixXd> ao3c = std::vector<Eigen::MatrixXd>(
548 auxshell.size(),
549 Eigen::MatrixXd::Zero(dftbasis.AOBasisSize(), dftbasis.AOBasisSize()));
550
551 std::vector<libint2::Shell> dftshells = dftbasis.GenerateLibintBasis();
552 std::vector<Index> shell2bf = dftbasis.getMapToBasisFunctions();
553
554 const libint2::Engine::target_ptr_vec& buf = engine.results();
555 // alpha-loop over the "left" DFT basis function
556 for (Index row = 0; row < Index(dftshells.size()); row++) {
557
558 const libint2::Shell& shell_row = dftshells[row];
559 const Index row_start = shell2bf[row];
560 // ThreecMatrix is symmetric, restrict explicit calculation to triangular
561 // matrix
562 for (Index col = 0; col <= row; col++) {
563 const libint2::Shell& shell_col = dftshells[col];
564 const Index col_start = shell2bf[col];
565
566 engine.compute2<libint2::Operator::coulomb, libint2::BraKet::xs_xx, 0>(
567 auxshell, libint2::Shell::unit(), shell_col, shell_row);
568
569 if (buf[0] == nullptr) {
570 continue;
571 }
572 Eigen::TensorMap<Eigen::Tensor<const double, 3, Eigen::RowMajor> const>
573 result(buf[0], auxshell.size(), shell_col.size(), shell_row.size());
574
575 for (size_t aux_c = 0; aux_c < auxshell.size(); aux_c++) {
576 for (size_t col_c = 0; col_c < shell_col.size(); col_c++) {
577 for (size_t row_c = 0; row_c < shell_row.size(); row_c++) {
578 // ao3c is col major result is row-major so this is ideal
579 ao3c[aux_c](row_start + row_c, col_start + col_c) =
580 result(aux_c, col_c, row_c);
581 } // ROW copy
582 } // COL copy
583 } // AUX copy
584
585 } // gamma-loop
586 } // alpha-loop
587
588 for (Eigen::MatrixXd& mat : ao3c) {
589 mat.triangularView<Eigen::Upper>() =
590 mat.triangularView<Eigen::Lower>().transpose();
591 }
592 return ao3c;
593}
594
595void TCMatrix_gwbse::Fill3cMO(const AOBasis& auxbasis, const AOBasis& dftbasis,
596 const Eigen::MatrixXd& dft_orbitals) {
597
598 const Eigen::MatrixXd dftm = dft_orbitals.middleCols(mmin_, mtotal_);
599 const Eigen::MatrixXd dftn =
600 dft_orbitals.middleCols(nmin_, ntotal_).transpose();
601
602 OpenMP_CUDA transform;
603 transform.setOperators(dftn, dftm);
604 Index nthreads = OPENMP::getMaxThreads();
605
606 std::vector<libint2::Shell> auxshells = auxbasis.GenerateLibintBasis();
607 std::vector<libint2::Engine> engines(nthreads);
608 engines[0] = libint2::Engine(
609 libint2::Operator::coulomb,
610 std::max(dftbasis.getMaxNprim(), auxbasis.getMaxNprim()),
611 static_cast<int>(std::max(dftbasis.getMaxL(), auxbasis.getMaxL())), 0);
612 engines[0].set(libint2::BraKet::xs_xx);
613 for (Index i = 1; i < nthreads; ++i) {
614 engines[i] = engines[0];
615 }
616 std::vector<Index> auxshell2bf = auxbasis.getMapToBasisFunctions();
617
618#pragma omp parallel
619 {
620 Index threadid = OPENMP::getThreadId();
621#pragma omp for schedule(dynamic)
622 for (Index aux = 0; aux < Index(auxshells.size()); aux++) {
623 const libint2::Shell& auxshell = auxshells[aux];
624
625 std::vector<Eigen::MatrixXd> ao3c =
626 ComputeAO3cBlock(auxshell, dftbasis, engines[threadid]);
627
628 // this is basically a transpose of AO3c and at the same time the ao->mo
629 // transformation
630 // we do not want to put it into matrix_ straight away is because,
631 // matrix_ is shared between all threads and we want a nice clean access
632 // pattern to it
633 std::vector<Eigen::MatrixXd> block = std::vector<Eigen::MatrixXd>(
634 mtotal_, Eigen::MatrixXd::Zero(ntotal_, ao3c.size()));
635
636 Index dim = static_cast<Index>(ao3c.size());
637 for (Index k = 0; k < dim; ++k) {
638 transform.MultiplyLeftRight(ao3c[k], threadid);
639 for (Index i = 0; i < ao3c[k].cols(); ++i) {
640 block[i].col(k) = ao3c[k].col(i);
641 }
642 }
643
644 // put into correct position
645 for (Index m_level = 0; m_level < mtotal_; m_level++) {
646 matrix_[m_level].middleCols(auxshell2bf[aux], auxshell.size()) =
647 block[m_level];
648 } // m-th DFT orbital
649 } // shells of GW basis set
650 }
651}
652
653} // namespace xtp
654} // namespace votca
Container to hold Basisfunctions for all atoms.
Definition aobasis.h:42
Index getMaxL() const
Definition aobasis.cc:29
AOShellIterator begin() const
Definition aobasis.h:49
Index AOBasisSize() const
Definition aobasis.h:46
AOShellIterator end() const
Definition aobasis.h:50
std::vector< Index > getMapToBasisFunctions() const
Definition aobasis.cc:45
Index getNumofShells() const
Definition aobasis.h:56
Index getMaxNprim() const
Definition aobasis.cc:37
std::vector< std::vector< Index > > ComputeShellPairs(double threshold=1e-20) const
std::vector< libint2::Shell > GenerateLibintBasis() const
Definition aobasis.cc:119
void Fill(const AOBasis &aobasis) final
Eigen::MatrixXd aomatrix_
Definition aomatrix.h:82
void computeCoulombIntegrals(const AOBasis &aobasis)
Eigen::MatrixXd Pseudo_InvSqrt(double etol)
Definition aomatrix.cc:88
Index Removedfunctions() const
Definition aomatrix.h:77
std::array< Eigen::MatrixXd, 3 > aomatrix_
Definition aomatrix.h:101
std::array< libint2::Shell::real_t, 3 > r_
Definition aomatrix.h:102
void Fill(const AOBasis &aobasis) final
void Fill(const AOBasis &aobasis) final
Eigen::MatrixXd aomatrix_
Definition aomatrix.h:44
void Fill(const AOBasis &aobasis) final
Eigen::MatrixXd singleShellOverlap(const AOShell &shell) const
Eigen::MatrixXd aomatrix_
Definition aomatrix.h:64
libint2::Shell LibintShell() const
Definition aoshell.cc:65
Index getSize() const
Definition aoshell.h:109
Index getNumFunc() const
Definition aoshell.h:103
std::vector< std::vector< libint2::ShellPair > > shellpairdata_
Definition ERIs.h:78
std::vector< libint2::Shell > basis_
Definition ERIs.h:74
Eigen::MatrixXd ComputeSchwarzShells(const AOBasis &dftbasis) const
std::vector< std::vector< Index > > shellpairs_
Definition ERIs.h:77
std::array< Eigen::MatrixXd, 2 > Compute4c(const Eigen::MatrixXd &dmat, double error) const
Index maxL_
Definition ERIs.h:80
Eigen::MatrixXd schwarzscreen_
Definition ERIs.h:98
Eigen::MatrixXd ComputeShellBlockNorm(const Eigen::MatrixXd &dmat) const
Definition ERIs.cc:63
std::vector< Index > starts_
Definition ERIs.h:75
Index maxnprim_
Definition ERIs.h:79
void MultiplyLeftRight(Eigen::MatrixXd &matrix, Index OpenmpThread)
void setOperators(const std::vector< Eigen::MatrixXd > &tensor, const Eigen::MatrixXd &rightoperator)
void Fill(const AOBasis &auxbasis, const AOBasis &dftbasis)
std::vector< Symmetric_Matrix > matrix_
Definition threecenter.h:63
void Fill3cMO(const AOBasis &auxbasis, const AOBasis &dftbasis, const Eigen::MatrixXd &dft_orbitals)
std::vector< Eigen::MatrixXd > matrix_
Eigen::MatrixXd inv_sqrt_
Definition threecenter.h:49
Index getMaxThreads()
Definition eigen.h:128
Index getThreadId()
Definition eigen.h:143
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > MatrixLibInt
Definition aoecp.cc:56
std::vector< Eigen::MatrixXd > ComputeAO3cBlock(const libint2::Shell &auxshell, const AOBasis &dftbasis, libint2::Engine &engine)
std::array< MatrixLibInt, libint2::operator_traits< obtype >::nopers > computeOneBodyIntegrals(const AOBasis &aobasis, OperatorParams oparams=OperatorParams())
base class for all analysis tools
Definition basebead.h:33
Eigen::Index Index
Definition types.h:26