29#define LIBINT2_CONSTEXPR_STATICS 0
31#include <libint2/statics_definition.h>
37 double threshold)
const {
44 std::vector<libint2::Engine> engines;
45 engines.reserve(nthreads);
47 for (
Index i = 1; i != nthreads; ++i) {
48 engines.push_back(engines[0]);
51 std::vector<std::vector<Index>> pairs(shells.size());
53#pragma omp parallel for schedule(dynamic)
54 for (
Index s1 = 0; s1 <
Index(shells.size()); ++s1) {
57 libint2::Engine& engine = engines[thread_id];
58 const libint2::Engine::target_ptr_vec& buf = engine.results();
59 Index n1 = shells[s1].size();
61 for (
Index s2 = 0; s2 <= s1; ++s2) {
62 bool on_same_center = (shells[s1].O == shells[s2].O);
63 bool significant = on_same_center;
64 if (!on_same_center) {
65 Index n2 = shells[s2].size();
66 engine.compute(shells[s1], shells[s2]);
67 Eigen::Map<const Eigen::MatrixXd> buf_mat(buf[0], n1, n2);
68 significant = (buf_mat.norm() >= threshold);
71 pairs[s1].push_back(s2);
74 std::sort(pairs[s1].
begin(), pairs[s1].
end());
80 Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
82template <libint2::Operator obtype,
83 typename OperatorParams =
84 typename libint2::operator_traits<obtype>::oper_params_type>
85std::array<MatrixLibInt, libint2::operator_traits<obtype>::nopers>
87 OperatorParams oparams = OperatorParams()) {
94 Index nopers =
static_cast<Index>(libint2::operator_traits<obtype>::nopers);
95 std::array<MatrixLibInt, libint2::operator_traits<obtype>::nopers> result;
100 std::vector<libint2::Engine> engines(nthreads);
101 engines[0] = libint2::Engine(obtype, aobasis.
getMaxNprim(),
102 static_cast<int>(aobasis.
getMaxL()), 0);
103 engines[0].set_params(oparams);
104 for (
Index i = 1; i < nthreads; ++i) {
105 engines[i] = engines[0];
110#pragma omp parallel for schedule(dynamic)
113 libint2::Engine& engine = engines[thread_id];
114 const libint2::Engine::target_ptr_vec& buf = engine.results();
116 Index bf1 = shell2bf[s1];
117 Index n1 = shells[s1].size();
119 for (
Index s2 : shellpair_list[s1]) {
121 engine.compute(shells[s1], shells[s2]);
122 if (buf[0] ==
nullptr) {
125 Index bf2 = shell2bf[s2];
126 Index n2 = shells[s2].size();
127 for (
unsigned int op = 0; op != nopers; ++op) {
128 Eigen::Map<const MatrixLibInt> buf_mat(buf[op], n1, n2);
129 result[op].block(bf1, bf2, n1, n2) = buf_mat;
132 result[op].block(bf2, bf1, n2, n1) = buf_mat.transpose();
144 aomatrix_ = computeOneBodyIntegrals<libint2::Operator::kinetic>(aobasis)[0];
151 aomatrix_ = computeOneBodyIntegrals<libint2::Operator::overlap>(aobasis)[0];
155 libint2::Shell::do_enforce_unit_normalization(
false);
156 libint2::Operator obtype = libint2::Operator::overlap;
158 libint2::Engine engine(obtype, shell.
getSize(),
159 static_cast<int>(shell.
getL()), 0);
161 const libint2::Engine::target_ptr_vec& buf = engine.results();
164 engine.compute(s, s);
166 Eigen::Map<const MatrixLibInt> buf_mat(buf[0], shell.
getNumFunc(),
169 }
catch (
const libint2::Engine::lmax_exceeded& error) {
170 std::ostringstream oss;
171 oss <<
"\nA libint error occured:\n"
172 << error.what() <<
"\n"
173 <<
"\nYou requested a computation for a shell with angular momentum "
174 << error.lmax_requested()
175 <<
",\nbut your libint package only supports angular momenta upto "
176 << error.lmax_limit() - 1 <<
".\n"
177 <<
"\nTo fix this error you will need to reinstall libint with "
179 "for higher angular momenta. If you installed your own libint it\n"
180 "should be reconfigured and installed with the option "
181 "--with-max-am=<maxAngularMomentum>.\n"
182 "If you installed libint with your OS package manager, you will\n"
183 "need to setup you own libint installation with the \n"
184 "--with-max-am=<maxAngularMomentum> option set."
186 throw std::runtime_error(oss.str());
195 std::array<libint2::Shell::real_t, 3>>(
198 for (
Index i = 0; i < 3; i++) {
220 std::vector<libint2::Engine> engines(nthreads);
222 libint2::Engine(libint2::Operator::coulomb, aobasis.
getMaxNprim(),
223 static_cast<int>(aobasis.
getMaxL()), 0);
224 engines[0].set(libint2::BraKet::xs_xs);
225 for (
Index i = 1; i != nthreads; ++i) {
226 engines[i] = engines[0];
229#pragma omp parallel for schedule(dynamic)
232 const libint2::Engine::target_ptr_vec& buf = engine.results();
234 Index bf1 = shell2bf[s1];
235 Index n1 = shells[s1].size();
238 for (
Index s2 = 0; s2 <= s1; ++s2) {
240 engine.compute2<libint2::Operator::coulomb, libint2::BraKet::xs_xs, 0>(
241 shells[s1], libint2::Shell::unit(), shells[s2],
242 libint2::Shell::unit());
244 if (buf[0] ==
nullptr) {
247 Index bf2 = shell2bf[s2];
248 Index n2 = shells[s2].size();
250 Eigen::Map<const MatrixLibInt> buf_mat(buf[0], n1, n2);
251 aomatrix_.block(bf1, bf2, n1, n2) = buf_mat;
254 aomatrix_.block(bf2, bf1, n2, n1) = buf_mat.transpose();
264 Eigen::MatrixXd result = Eigen::MatrixXd::Zero(noshells, noshells);
266 std::vector<libint2::Engine> engines(nthreads);
267 double epsilon = 0.0;
268 engines[0] = libint2::Engine(libint2::Operator::coulomb, basis.
getMaxNprim(),
269 static_cast<int>(basis.
getMaxL()), 0, epsilon);
271 for (
Index i = 1; i < nthreads; ++i) {
272 engines[i] = engines[0];
277#pragma omp parallel for schedule(dynamic)
280 libint2::Engine& engine = engines[thread_id];
281 const libint2::Engine::target_ptr_vec& buf = engine.results();
282 Index n1 = shells[s1].size();
284 for (
Index s2 = 0l; s2 <= s1; ++s2) {
285 Index n2 = shells[s2].size();
289 .compute2<libint2::Operator::coulomb, libint2::BraKet::xx_xx, 0>(
290 shells[s1], shells[s2], shells[s1], shells[s2]);
292 Eigen::Map<const MatrixLibInt> buf_mat(buf[0], n12, n12);
294 result(s2, s1) = std::sqrt(buf_mat.cwiseAbs().maxCoeff());
297 return result.selfadjointView<Eigen::Upper>();
300template <
bool with_exchange>
302 double error)
const {
304 "Please call Initialize_4c before running this");
307 Eigen::MatrixXd hartree = Eigen::MatrixXd::Zero(dmat.rows(), dmat.cols());
308 Eigen::MatrixXd exchange;
310 exchange = Eigen::MatrixXd::Zero(dmat.rows(), dmat.cols());
313 double fock_precision = error;
317 double engine_precision = std::min(fock_precision / dnorm_block.maxCoeff(),
318 std::numeric_limits<double>::epsilon()) /
320 std::vector<libint2::Engine> engines(nthreads);
321 engines[0] = libint2::Engine(libint2::Operator::coulomb,
int(
maxnprim_),
323 engines[0].set_precision(engine_precision);
327 for (
Index i = 1; i < nthreads; ++i) {
328 engines[i] = engines[0];
332#pragma omp parallel for schedule(dynamic) reduction(+ : hartree) \
333 reduction(+ : exchange)
334 for (
Index s1 = 0; s1 < nshells; ++s1) {
336 libint2::Engine& engine = engines[thread_id];
337 const auto& buf = engine.results();
339 const libint2::Shell& shell1 =
basis_[s1];
340 Index n1 = shell1.size();
345 const libint2::Shell& shell2 =
basis_[s2];
346 Index n2 = shell2.size();
347 double dnorm_12 = dnorm_block(s1, s2);
348 const libint2::ShellPair* sp12 = &(*sp12_iter);
351 for (
Index s3 = 0; s3 <= s1; ++s3) {
354 const libint2::Shell& shell3 =
basis_[s3];
355 Index n3 = shell3.size();
357 double dnorm_123 = std::max(dnorm_block(s1, s3),
358 std::max(dnorm_block(s2, s3), dnorm_12));
359 Index s4max = (s1 == s3) ? s2 : s3;
366 const libint2::ShellPair* sp34 = &(*sp34_iter);
370 std::max(dnorm_block(s1, s4),
371 std::max(dnorm_block(s2, s4),
372 std::max(dnorm_block(s3, s4), dnorm_123)));
379 const libint2::Shell& shell4 =
basis_[s4];
381 .compute2<libint2::Operator::coulomb, libint2::BraKet::xx_xx, 0>(
382 shell1, shell2, shell3, shell4, sp12, sp34);
383 const auto* buf_1234 = buf[0];
384 if (buf_1234 ==
nullptr) {
388 Index n4 = shell4.size();
389 Index s12_deg = (s1 == s2) ? 1 : 2;
390 Index s34_deg = (s3 == s4) ? 1 : 2;
391 Index s12_34_deg = (s1 == s3) ? (s2 == s4 ? 1 : 2) : 2;
392 Index s1234_deg = s12_deg * s34_deg * s12_34_deg;
394 for (
Index f1 = 0, f1234 = 0; f1 != n1; ++f1) {
395 const Index bf1 = f1 + start_1;
396 for (
Index f2 = 0; f2 != n2; ++f2) {
397 const Index bf2 = f2 + start_2;
398 for (
Index f3 = 0; f3 != n3; ++f3) {
399 const Index bf3 = f3 + start_3;
400 for (
Index f4 = 0; f4 != n4; ++f4, ++f1234) {
401 const Index bf4 = f4 + start_4;
403 const double value = buf_1234[f1234];
405 const double value_scal_by_deg = value * double(s1234_deg);
407 hartree(bf1, bf2) += dmat(bf3, bf4) * value_scal_by_deg;
408 hartree(bf3, bf4) += dmat(bf1, bf2) * value_scal_by_deg;
410 exchange(bf1, bf3) -= dmat(bf2, bf4) * value_scal_by_deg;
411 exchange(bf2, bf3) -= dmat(bf1, bf4) * value_scal_by_deg;
412 exchange(bf2, bf4) -= dmat(bf1, bf3) * value_scal_by_deg;
413 exchange(bf1, bf4) -= dmat(bf2, bf3) * value_scal_by_deg;
423 std::array<Eigen::MatrixXd, 2> result2;
425 result2[0] = 0.25 * (hartree + hartree.transpose());
428 result2[1] = 0.125 * (exchange + exchange.transpose());
433template std::array<Eigen::MatrixXd, 2> ERIs::Compute4c<true>(
434 const Eigen::MatrixXd& dmat,
double error)
const;
435template std::array<Eigen::MatrixXd, 2> ERIs::Compute4c<false>(
436 const Eigen::MatrixXd& dmat,
double error)
const;
441 auxAOcoulomb.
Fill(auxbasis);
447#pragma omp parallel for schedule(dynamic, 4)
455 std::vector<libint2::Engine> engines(nthreads);
456 engines[0] = libint2::Engine(
457 libint2::Operator::coulomb,
459 static_cast<int>(std::max(dftbasis.
getMaxL(), auxbasis.
getMaxL())), 0);
460 engines[0].set(libint2::BraKet::xs_xx);
461 for (
Index i = 1; i < nthreads; ++i) {
462 engines[i] = engines[0];
468#pragma omp parallel for schedule(dynamic)
472 const libint2::Engine::target_ptr_vec& buf = engine.results();
473 const libint2::Shell& dftshell = dftshells[is];
474 Index start = shell2bf[is];
475 std::vector<Eigen::MatrixXd> block(dftshell.size());
476 for (
Index i = 0; i <
Index(dftshell.size()); i++) {
482 const libint2::Shell& auxshell = auxshells[aux];
483 Index aux_start = auxshell2bf[aux];
485 for (
Index dis = 0; dis <= is; dis++) {
487 const libint2::Shell& shell_col = dftshells[dis];
488 Index col_start = shell2bf[dis];
489 engine.compute2<libint2::Operator::coulomb, libint2::BraKet::xs_xx, 0>(
490 auxshell, libint2::Shell::unit(), dftshell, shell_col);
492 if (buf[0] ==
nullptr) {
495 Eigen::TensorMap<Eigen::Tensor<const double, 3, Eigen::RowMajor>
const>
496 result(buf[0], auxshell.size(), dftshell.size(), shell_col.size());
498 for (
size_t left = 0; left < dftshell.size(); left++) {
499 for (
size_t auxf = 0; auxf < auxshell.size(); auxf++) {
500 for (
size_t col = 0; col < shell_col.size(); col++) {
502 if ((col_start + col) > (start + left)) {
505 block[left](aux_start + auxf, col_start + col) =
506 result(auxf, left, col);
513 for (
Index i = 0; i <
Index(block.size()); ++i) {
514 Eigen::MatrixXd temp =
inv_sqrt_ * block[i];
515 for (
Index mu = 0; mu < temp.rows(); ++mu) {
516 for (
Index j = 0; j < temp.cols(); ++j) {
517 matrix_[mu](i + start, j) = temp(mu, j);
533 libint2::Engine& engine) {
534 std::vector<Eigen::MatrixXd> ao3c = std::vector<Eigen::MatrixXd>(
541 const libint2::Engine::target_ptr_vec& buf = engine.results();
543 for (
Index row = 0; row <
Index(dftshells.size()); row++) {
545 const libint2::Shell& shell_row = dftshells[row];
546 const Index row_start = shell2bf[row];
549 for (
Index col = 0; col <= row; col++) {
550 const libint2::Shell& shell_col = dftshells[col];
551 const Index col_start = shell2bf[col];
553 engine.compute2<libint2::Operator::coulomb, libint2::BraKet::xs_xx, 0>(
554 auxshell, libint2::Shell::unit(), shell_col, shell_row);
556 if (buf[0] ==
nullptr) {
559 Eigen::TensorMap<Eigen::Tensor<const double, 3, Eigen::RowMajor>
const>
560 result(buf[0], auxshell.size(), shell_col.size(), shell_row.size());
562 for (
size_t aux_c = 0; aux_c < auxshell.size(); aux_c++) {
563 for (
size_t col_c = 0; col_c < shell_col.size(); col_c++) {
564 for (
size_t row_c = 0; row_c < shell_row.size(); row_c++) {
566 ao3c[aux_c](row_start + row_c, col_start + col_c) =
567 result(aux_c, col_c, row_c);
575 for (Eigen::MatrixXd& mat : ao3c) {
576 mat.triangularView<Eigen::Upper>() =
577 mat.triangularView<Eigen::Lower>().transpose();
583 const Eigen::MatrixXd& dft_orbitals) {
585 const Eigen::MatrixXd dftm = dft_orbitals.middleCols(
mmin_,
mtotal_);
586 const Eigen::MatrixXd dftn =
594 std::vector<libint2::Engine> engines(nthreads);
595 engines[0] = libint2::Engine(
596 libint2::Operator::coulomb,
598 static_cast<int>(std::max(dftbasis.
getMaxL(), auxbasis.
getMaxL())), 0);
599 engines[0].set(libint2::BraKet::xs_xx);
600 for (
Index i = 1; i < nthreads; ++i) {
601 engines[i] = engines[0];
608#pragma omp for schedule(dynamic)
609 for (
Index aux = 0; aux <
Index(auxshells.size()); aux++) {
610 const libint2::Shell& auxshell = auxshells[aux];
612 std::vector<Eigen::MatrixXd> ao3c =
620 std::vector<Eigen::MatrixXd> block = std::vector<Eigen::MatrixXd>(
624 for (
Index k = 0; k < dim; ++k) {
626 for (
Index i = 0; i < ao3c[k].cols(); ++i) {
627 block[i].col(k) = ao3c[k].col(i);
633 matrix_[m_level].middleCols(auxshell2bf[aux], auxshell.size()) =
Container to hold Basisfunctions for all atoms.
AOShellIterator begin() const
Index AOBasisSize() const
AOShellIterator end() const
std::vector< Index > getMapToBasisFunctions() const
Index getNumofShells() const
Index getMaxNprim() const
std::vector< std::vector< Index > > ComputeShellPairs(double threshold=1e-20) const
std::vector< libint2::Shell > GenerateLibintBasis() const
void Fill(const AOBasis &aobasis) final
Eigen::MatrixXd aomatrix_
void computeCoulombIntegrals(const AOBasis &aobasis)
Eigen::MatrixXd Pseudo_InvSqrt(double etol)
Index Removedfunctions() const
std::array< Eigen::MatrixXd, 3 > aomatrix_
std::array< libint2::Shell::real_t, 3 > r_
void Fill(const AOBasis &aobasis) final
void Fill(const AOBasis &aobasis) final
Eigen::MatrixXd aomatrix_
void Fill(const AOBasis &aobasis) final
Eigen::MatrixXd singleShellOverlap(const AOShell &shell) const
Eigen::MatrixXd aomatrix_
libint2::Shell LibintShell() const
std::vector< std::vector< libint2::ShellPair > > shellpairdata_
std::vector< libint2::Shell > basis_
Eigen::MatrixXd ComputeSchwarzShells(const AOBasis &dftbasis) const
std::vector< std::vector< Index > > shellpairs_
std::array< Eigen::MatrixXd, 2 > Compute4c(const Eigen::MatrixXd &dmat, double error) const
Eigen::MatrixXd schwarzscreen_
Eigen::MatrixXd ComputeShellBlockNorm(const Eigen::MatrixXd &dmat) const
std::vector< Index > starts_
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_
void Fill3cMO(const AOBasis &auxbasis, const AOBasis &dftbasis, const Eigen::MatrixXd &dft_orbitals)
std::vector< Eigen::MatrixXd > matrix_
Eigen::MatrixXd inv_sqrt_
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > MatrixLibInt
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