50 double threshold)
const {
57 std::vector<libint2::Engine> engines;
58 engines.reserve(nthreads);
60 for (
Index i = 1; i != nthreads; ++i) {
61 engines.push_back(engines[0]);
64 std::vector<std::vector<Index>> pairs(shells.size());
66#pragma omp parallel for schedule(dynamic)
67 for (
Index s1 = 0; s1 <
Index(shells.size()); ++s1) {
70 libint2::Engine& engine = engines[thread_id];
71 const libint2::Engine::target_ptr_vec& buf = engine.results();
72 Index n1 = shells[s1].size();
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);
84 pairs[s1].push_back(s2);
87 std::sort(pairs[s1].
begin(), pairs[s1].
end());
100 OperatorParams oparams = OperatorParams()) {
107 Index nopers =
static_cast<Index>(libint2::operator_traits<obtype>::nopers);
108 std::array<MatrixLibInt, libint2::operator_traits<obtype>::nopers> result;
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];
123#pragma omp parallel for schedule(dynamic)
126 libint2::Engine& engine = engines[thread_id];
127 const libint2::Engine::target_ptr_vec& buf = engine.results();
129 Index bf1 = shell2bf[s1];
130 Index n1 = shells[s1].size();
132 for (
Index s2 : shellpair_list[s1]) {
134 engine.compute(shells[s1], shells[s2]);
135 if (buf[0] ==
nullptr) {
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;
145 result[op].block(bf2, bf1, n2, n1) = buf_mat.transpose();
168 libint2::Shell::do_enforce_unit_normalization(
false);
169 libint2::Operator obtype = libint2::Operator::overlap;
171 libint2::Engine engine(obtype, shell.
getSize(),
172 static_cast<int>(shell.
getL()), 0);
174 const libint2::Engine::target_ptr_vec& buf = engine.results();
177 engine.compute(s, s);
179 Eigen::Map<const MatrixLibInt> buf_mat(buf[0], shell.
getNumFunc(),
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 "
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."
199 throw std::runtime_error(oss.str());
233 std::vector<libint2::Engine> engines(nthreads);
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];
242#pragma omp parallel for schedule(dynamic)
245 const libint2::Engine::target_ptr_vec& buf = engine.results();
247 Index bf1 = shell2bf[s1];
248 Index n1 = shells[s1].size();
251 for (
Index s2 = 0; s2 <= s1; ++s2) {
253 engine.compute2<libint2::Operator::coulomb, libint2::BraKet::xs_xs, 0>(
254 shells[s1], libint2::Shell::unit(), shells[s2],
255 libint2::Shell::unit());
257 if (buf[0] ==
nullptr) {
260 Index bf2 = shell2bf[s2];
261 Index n2 = shells[s2].size();
263 Eigen::Map<const MatrixLibInt> buf_mat(buf[0], n1, n2);
264 aomatrix_.block(bf1, bf2, n1, n2) = buf_mat;
267 aomatrix_.block(bf2, bf1, n2, n1) = buf_mat.transpose();
277 Eigen::MatrixXd result = Eigen::MatrixXd::Zero(noshells, noshells);
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);
284 for (
Index i = 1; i < nthreads; ++i) {
285 engines[i] = engines[0];
290#pragma omp parallel for schedule(dynamic)
293 libint2::Engine& engine = engines[thread_id];
294 const libint2::Engine::target_ptr_vec& buf = engine.results();
295 Index n1 = shells[s1].size();
297 for (
Index s2 = 0l; s2 <= s1; ++s2) {
298 Index n2 = shells[s2].size();
302 .compute2<libint2::Operator::coulomb, libint2::BraKet::xx_xx, 0>(
303 shells[s1], shells[s2], shells[s1], shells[s2]);
305 Eigen::Map<const MatrixLibInt> buf_mat(buf[0], n12, n12);
307 result(s2, s1) = std::sqrt(buf_mat.cwiseAbs().maxCoeff());
310 return result.selfadjointView<Eigen::Upper>();
315 double error)
const {
317 "Please call Initialize_4c before running this");
320 Eigen::MatrixXd hartree = Eigen::MatrixXd::Zero(dmat.rows(), dmat.cols());
321 Eigen::MatrixXd exchange;
323 exchange = Eigen::MatrixXd::Zero(dmat.rows(), dmat.cols());
326 double fock_precision = error;
330 double engine_precision = std::min(fock_precision / dnorm_block.maxCoeff(),
331 std::numeric_limits<double>::epsilon()) /
333 std::vector<libint2::Engine> engines(nthreads);
334 engines[0] = libint2::Engine(libint2::Operator::coulomb,
int(
maxnprim_),
336 engines[0].set_precision(engine_precision);
340 for (
Index i = 1; i < nthreads; ++i) {
341 engines[i] = engines[0];
345#pragma omp parallel for schedule(dynamic) reduction(+ : hartree) \
346 reduction(+ : exchange)
347 for (
Index s1 = 0; s1 < nshells; ++s1) {
349 libint2::Engine& engine = engines[thread_id];
350 const auto& buf = engine.results();
352 const libint2::Shell& shell1 =
basis_[s1];
353 Index n1 = shell1.size();
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);
364 for (
Index s3 = 0; s3 <= s1; ++s3) {
367 const libint2::Shell& shell3 =
basis_[s3];
368 Index n3 = shell3.size();
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;
379 const libint2::ShellPair* sp34 = &(*sp34_iter);
383 std::max(dnorm_block(s1, s4),
384 std::max(dnorm_block(s2, s4),
385 std::max(dnorm_block(s3, s4), dnorm_123)));
392 const libint2::Shell& shell4 =
basis_[s4];
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) {
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;
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;
416 const double value = buf_1234[f1234];
418 const double value_scal_by_deg = value * double(s1234_deg);
420 hartree(bf1, bf2) += dmat(bf3, bf4) * value_scal_by_deg;
421 hartree(bf3, bf4) += dmat(bf1, bf2) * value_scal_by_deg;
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;
436 std::array<Eigen::MatrixXd, 2> result2;
438 result2[0] = 0.25 * (hartree + hartree.transpose());
441 result2[1] = 0.125 * (exchange + exchange.transpose());
454 auxAOcoulomb.
Fill(auxbasis);
460#pragma omp parallel for schedule(dynamic, 4)
468 std::vector<libint2::Engine> engines(nthreads);
469 engines[0] = libint2::Engine(
470 libint2::Operator::coulomb,
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];
481#pragma omp parallel for schedule(dynamic)
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++) {
495 const libint2::Shell& auxshell = auxshells[aux];
496 Index aux_start = auxshell2bf[aux];
498 for (
Index dis = 0; dis <= is; dis++) {
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);
505 if (buf[0] ==
nullptr) {
508 Eigen::TensorMap<Eigen::Tensor<const double, 3, Eigen::RowMajor>
const>
509 result(buf[0], auxshell.size(), dftshell.size(), shell_col.size());
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++) {
515 if ((col_start + col) > (start + left)) {
518 block[left](aux_start + auxf, col_start + col) =
519 result(auxf, left, col);
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);
546 libint2::Engine& engine) {
547 std::vector<Eigen::MatrixXd> ao3c = std::vector<Eigen::MatrixXd>(
554 const libint2::Engine::target_ptr_vec& buf = engine.results();
556 for (
Index row = 0; row <
Index(dftshells.size()); row++) {
558 const libint2::Shell& shell_row = dftshells[row];
559 const Index row_start = shell2bf[row];
562 for (
Index col = 0; col <= row; col++) {
563 const libint2::Shell& shell_col = dftshells[col];
564 const Index col_start = shell2bf[col];
566 engine.compute2<libint2::Operator::coulomb, libint2::BraKet::xs_xx, 0>(
567 auxshell, libint2::Shell::unit(), shell_col, shell_row);
569 if (buf[0] ==
nullptr) {
572 Eigen::TensorMap<Eigen::Tensor<const double, 3, Eigen::RowMajor>
const>
573 result(buf[0], auxshell.size(), shell_col.size(), shell_row.size());
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++) {
579 ao3c[aux_c](row_start + row_c, col_start + col_c) =
580 result(aux_c, col_c, row_c);
588 for (Eigen::MatrixXd& mat : ao3c) {
589 mat.triangularView<Eigen::Upper>() =
590 mat.triangularView<Eigen::Lower>().transpose();
596 const Eigen::MatrixXd& dft_orbitals) {
598 const Eigen::MatrixXd dftm = dft_orbitals.middleCols(
mmin_,
mtotal_);
599 const Eigen::MatrixXd dftn =
607 std::vector<libint2::Engine> engines(nthreads);
608 engines[0] = libint2::Engine(
609 libint2::Operator::coulomb,
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];
621#pragma omp for schedule(dynamic)
622 for (
Index aux = 0; aux <
Index(auxshells.size()); aux++) {
623 const libint2::Shell& auxshell = auxshells[aux];
625 std::vector<Eigen::MatrixXd> ao3c =
633 std::vector<Eigen::MatrixXd> block = std::vector<Eigen::MatrixXd>(
637 for (
Index k = 0; k < dim; ++k) {
639 for (
Index i = 0; i < ao3c[k].cols(); ++i) {
640 block[i].col(k) = ao3c[k].col(i);
646 matrix_[m_level].middleCols(auxshell2bf[aux], auxshell.size()) =