33 #include <unordered_map> 37 #ifdef LIBINT2_HAVE_BTAS 38 # include <btas/btas.h> 39 #else // LIBINT2_HAVE_BTAS 40 # error "libint2::lcao requires BTAS" 44 #include <libint2/diis.h> 45 #include <libint2/util/intpart_iter.h> 46 #include <libint2/chemistry/sto3g_atomic_density.h> 47 #include <libint2.hpp> 53 typedef btas::RangeNd<CblasRowMajor, std::array<long, 2>> Range2;
54 typedef btas::RangeNd<CblasRowMajor, std::array<long, 3>> Range3;
55 typedef btas::RangeNd<CblasRowMajor, std::array<long, 4>> Range4;
56 typedef btas::Tensor<double, Range2> Tensor2d;
57 typedef btas::Tensor<double, Range3> Tensor3d;
58 typedef btas::Tensor<double, Range3> Tensor4d;
60 typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
65 typedef Eigen::DiagonalMatrix<double, Eigen::Dynamic, Eigen::Dynamic>
69 using libint2::libint2::Atom;
70 using libint2::BasisSet;
71 using libint2::Operator;
72 using libint2::BraKet;
74 std::vector<libint2::Atom> read_geometry(
const std::string& filename);
75 Matrix compute_soad(
const std::vector<libint2::Atom>& atoms);
77 Matrix compute_shellblock_norm(
const BasisSet& obs,
const Matrix& A);
79 template <Operator obtype>
81 const BasisSet& obs,
const std::vector<libint2::Atom>& atoms = std::vector<libint2::Atom>());
83 #if LIBINT2_DERIV_ONEBODY_ORDER 84 template <Operator obtype>
85 std::vector<Matrix> compute_1body_ints_deriv(
unsigned deriv_order,
87 const std::vector<libint2::Atom>& atoms);
88 #endif // LIBINT2_DERIV_ONEBODY_ORDER 90 template <lib
int2::Operator Kernel = lib
int2::Operator::coulomb>
91 Matrix compute_schwarz_ints(
92 const BasisSet& bs1,
const BasisSet& bs2 = BasisSet(),
93 bool use_2norm =
false,
94 typename libint2::operator_traits<Kernel>::oper_params_type params =
95 libint2::operator_traits<Kernel>::default_params());
96 Matrix compute_do_ints(
const BasisSet& bs1,
const BasisSet& bs2 = BasisSet(),
97 bool use_2norm =
false 100 using shellpair_list_t = std::unordered_map<size_t, std::vector<size_t>>;
101 shellpair_list_t obs_shellpair_list;
107 shellpair_list_t compute_shellpair_list(
const BasisSet& bs1,
108 const BasisSet& bs2 = BasisSet(),
109 double threshold = 1e-12);
111 Matrix compute_2body_fock(
112 const BasisSet& obs,
const Matrix& D,
113 double precision = std::numeric_limits<
115 const Matrix& Schwarz = Matrix()
119 Matrix compute_2body_fock_general(
120 const BasisSet& obs,
const Matrix& D,
const BasisSet& D_bs,
121 bool D_is_sheldiagonal =
false,
122 double precision = std::numeric_limits<
126 #if LIBINT2_DERIV_ERI_ORDER 127 template <
unsigned deriv_order>
128 std::vector<Matrix> compute_2body_fock_deriv(
129 const BasisSet& obs,
const std::vector<libint2::Atom>& atoms,
131 double precision = std::numeric_limits<
133 const Matrix& Schwarz = Matrix()
136 #endif // LIBINT2_DERIV_ERI_ORDER 143 std::tuple<Matrix, Matrix, double> conditioning_orthogonalizer(
144 const Matrix& S,
double S_condition_number_threshold);
146 #ifdef LIBINT2_HAVE_BTAS 147 #define HAVE_DENSITY_FITTING 1 150 const BasisSet& dfbs;
151 DFFockEngine(
const BasisSet& _obs,
const BasisSet& _dfbs)
152 : obs(_obs), dfbs(_dfbs) {}
154 typedef btas::RangeNd<CblasRowMajor, std::array<long, 3>> Range3d;
155 typedef btas::Tensor<double, Range3d> Tensor3d;
159 Matrix compute_2body_fock_dfC(
const Matrix& Cocc);
161 #endif // HAVE_DENSITY_FITTING 167 template <
typename Lambda>
172 auto thread_id = omp_get_thread_num();
175 #else // use C++11 threads 176 std::vector<std::thread> threads;
177 for (
int thread_id = 0; thread_id != libint2::nthreads; ++thread_id) {
178 if (thread_id != nthreads - 1)
179 threads.push_back(std::thread(lambda, thread_id));
183 for (
int thread_id = 0; thread_id < nthreads - 1; ++thread_id)
184 threads[thread_id].join();
189 int main(
int argc,
char* argv[]) {
201 const auto filename = (argc > 1) ? argv[1] :
"h2o.xyz";
202 const auto basisname = (argc > 2) ? argv[2] :
"aug-cc-pVDZ";
203 bool do_density_fitting =
false;
204 #ifdef HAVE_DENSITY_FITTING 205 do_density_fitting = (argc > 3);
206 const auto dfbasisname = do_density_fitting ? argv[3] :
"";
208 std::vector<libint2::Atom> atoms = read_geometry(filename);
212 using libint2::nthreads;
213 auto nthreads_cstr = getenv(
"LIBINT_NUM_THREADS");
215 if (nthreads_cstr && strcmp(nthreads_cstr,
"")) {
216 std::istringstream iss(nthreads_cstr);
218 if (nthreads > 1 << 16 || nthreads <= 0) nthreads = 1;
221 omp_set_num_threads(nthreads);
223 std::cout <<
"Will scale over " << nthreads
229 <<
" threads" << std::endl;
234 for (
auto i = 0; i < atoms.size(); ++i) nelectron += atoms[i].atomic_number;
235 const auto ndocc = nelectron / 2;
236 cout <<
"# of electrons = " << nelectron << endl;
240 for (
auto i = 0; i < atoms.size(); i++)
241 for (
auto j = i + 1; j < atoms.size(); j++) {
242 auto xij = atoms[i].x - atoms[j].x;
243 auto yij = atoms[i].y - atoms[j].y;
244 auto zij = atoms[i].z - atoms[j].z;
245 auto r2 = xij * xij + yij * yij + zij * zij;
247 enuc += atoms[i].atomic_number * atoms[j].atomic_number / r;
249 cout <<
"Nuclear repulsion energy = " << std::setprecision(15) << enuc
254 cout <<
"libint2::Atomic Cartesian coordinates (a.u.):" << endl;
255 for (
const auto& a : atoms)
256 std::cout << a.atomic_number <<
" " << a.x <<
" " << a.y <<
" " << a.z
259 BasisSet obs(basisname, atoms);
260 cout <<
"orbital basis set rank = " << obs.nbf() << endl;
262 #ifdef HAVE_DENSITY_FITTING 264 if (do_density_fitting) {
265 dfbs = BasisSet(dfbasisname, atoms);
266 cout <<
"density-fitting basis set rank = " << dfbs.nbf() << endl;
268 #endif // HAVE_DENSITY_FITTING 279 obs_shellpair_list = compute_shellpair_list(obs);
281 for (
auto& sp : obs_shellpair_list) {
282 nsp += sp.second.size();
284 std::cout <<
"# of {all,non-negligible} shell-pairs = {" 285 << obs.size() * (obs.size() + 1) / 2 <<
"," << nsp <<
"}" 290 auto S = compute_1body_ints<Operator::overlap>(obs)[0];
291 auto T = compute_1body_ints<Operator::kinetic>(obs)[0];
299 double XtX_condition_number;
308 double S_condition_number_threshold =
309 1.0 / std::numeric_limits<double>::epsilon();
310 std::tie(X, Xinv, XtX_condition_number) =
311 conditioning_orthogonalizer(S, S_condition_number_threshold);
317 const auto tstart = std::chrono::high_resolution_clock::now();
319 auto D_minbs = compute_soad(atoms);
320 BasisSet minbs(
"STO-3G", atoms);
326 std::cout <<
"projecting SOAD into AO basis ... ";
328 F += compute_2body_fock_general(
329 obs, D_minbs, minbs,
true ,
330 std::numeric_limits<double>::epsilon()
337 Eigen::SelfAdjointEigenSolver<Matrix> eig_solver(X.transpose() * F * X);
338 auto C = X * eig_solver.eigenvectors();
341 C_occ = C.leftCols(ndocc);
342 D = C_occ * C_occ.transpose();
344 const auto tstop = std::chrono::high_resolution_clock::now();
345 const std::chrono::duration<double> time_elapsed = tstop - tstart;
346 std::cout <<
"done (" << time_elapsed.count() <<
" s)" << std::endl;
351 auto K = compute_schwarz_ints<>(obs, BasisSet(),
true);
354 #ifdef HAVE_DENSITY_FITTING 355 std::unique_ptr<DFFockEngine> dffockengine(
356 do_density_fitting ?
new DFFockEngine(obs, dfbs) :
nullptr);
357 #endif // HAVE_DENSITY_FITTING 363 const auto maxiter = 100;
364 const auto conv = 1e-12;
366 auto rms_error = 1.0;
367 auto ediff_rel = 0.0;
369 auto n2 = D.cols() * D.rows();
375 bool reset_incremental_fock_formation =
false;
376 bool incremental_Fbuild_started =
false;
377 double start_incremental_F_threshold = 1e-5;
378 double next_reset_threshold = 0.0;
379 size_t last_reset_iteration = 0;
381 if (do_density_fitting) start_incremental_F_threshold = 0.0;
384 const auto tstart = std::chrono::high_resolution_clock::now();
391 if (not incremental_Fbuild_started &&
392 rms_error < start_incremental_F_threshold) {
393 incremental_Fbuild_started =
true;
394 reset_incremental_fock_formation =
false;
395 last_reset_iteration = iter - 1;
396 next_reset_threshold = rms_error / 1e1;
397 std::cout <<
"== started incremental fock build" << std::endl;
399 if (reset_incremental_fock_formation || not incremental_Fbuild_started) {
403 if (reset_incremental_fock_formation && incremental_Fbuild_started) {
404 reset_incremental_fock_formation =
false;
405 last_reset_iteration = iter;
406 next_reset_threshold = rms_error / 1e1;
407 std::cout <<
"== reset incremental fock build" << std::endl;
411 if (not do_density_fitting) {
413 const auto precision_F = std::min(
414 std::min(1e-3 / XtX_condition_number, 1e-7),
415 std::max(rms_error / 1e4, std::numeric_limits<double>::epsilon()));
416 F += compute_2body_fock(obs, D_diff, precision_F, K);
418 #if HAVE_DENSITY_FITTING 420 F = H + dffockengine->compute_2body_fock_dfC(C_occ);
427 #endif // HAVE_DENSITY_FITTING 430 ehf = D.cwiseProduct(H + F).sum();
431 ediff_rel = std::abs((ehf - ehf_last) / ehf);
434 Matrix FD_comm = F * D * S - S * D * F;
435 rms_error = FD_comm.norm() / n2;
436 if (rms_error < next_reset_threshold || iter - last_reset_iteration >= 8)
437 reset_incremental_fock_formation =
true;
443 diis.extrapolate(F_diis, FD_comm);
448 Eigen::SelfAdjointEigenSolver<Matrix> eig_solver(X.transpose() * F_diis *
450 evals = eig_solver.eigenvalues();
451 auto C = X * eig_solver.eigenvectors();
454 C_occ = C.leftCols(ndocc);
455 D = C_occ * C_occ.transpose();
458 const auto tstop = std::chrono::high_resolution_clock::now();
459 const std::chrono::duration<double> time_elapsed = tstop - tstart;
462 std::cout <<
"\n\nIter E(HF) D(E)/E " 463 "RMS([F,D])/nn Time(s)\n";
464 printf(
" %02d %20.12f %20.12e %20.12e %10.5lf\n", iter, ehf + enuc,
465 ediff_rel, rms_error, time_elapsed.count());
467 }
while (((ediff_rel > conv) || (rms_error > conv)) && (iter < maxiter));
469 printf(
"** Hartree-Fock energy = %20.12f\n", ehf + enuc);
471 auto Mu = compute_1body_ints<Operator::emultipole2>(obs);
474 for (
int xyz = 0; xyz != 3; ++xyz)
476 D.cwiseProduct(Mu[xyz + 1])
478 for (
int k = 0; k != 6; ++k)
480 D.cwiseProduct(Mu[k + 4])
482 std::cout <<
"** edipole = ";
483 std::copy(mu.begin(), mu.end(),
484 std::ostream_iterator<double>(std::cout,
" "));
485 std::cout << std::endl;
486 std::cout <<
"** equadrupole = ";
487 std::copy(qu.begin(), qu.end(),
488 std::ostream_iterator<double>(std::cout,
" "));
489 std::cout << std::endl;
492 #if LIBINT2_DERIV_ONEBODY_ORDER 494 Matrix F1 = Matrix::Zero(atoms.size(), 3);
495 Matrix F_Pulay = Matrix::Zero(atoms.size(), 3);
499 auto T1 = compute_1body_ints_deriv<Operator::kinetic>(1, obs, atoms);
500 auto V1 = compute_1body_ints_deriv<Operator::nuclear>(1, obs, atoms);
501 for (
auto atom = 0, i = 0; atom != atoms.size(); ++atom) {
502 for (
auto xyz = 0; xyz != 3; ++xyz, ++i) {
503 auto force = 2 * (T1[i] + V1[i]).cwiseProduct(D).sum();
504 F1(atom, xyz) += force;
512 DiagonalMatrix evals_occ(evals.topRows(ndocc));
513 Matrix W = C_occ * evals_occ * C_occ.transpose();
514 auto S1 = compute_1body_ints_deriv<Operator::overlap>(1, obs, atoms);
515 for (
auto atom = 0, i = 0; atom != atoms.size(); ++atom) {
516 for (
auto xyz = 0; xyz != 3; ++xyz, ++i) {
517 auto force = 2 * S1[i].cwiseProduct(W).sum();
518 F_Pulay(atom, xyz) -= force;
522 std::cout <<
"** 1-body forces = ";
523 for (
int atom = 0; atom != atoms.size(); ++atom)
524 for (
int xyz = 0; xyz != 3; ++xyz) std::cout << F1(atom, xyz) <<
" ";
525 std::cout << std::endl;
526 std::cout <<
"** Pulay forces = ";
527 for (
int atom = 0; atom != atoms.size(); ++atom)
528 for (
int xyz = 0; xyz != 3; ++xyz)
529 std::cout << F_Pulay(atom, xyz) <<
" ";
530 std::cout << std::endl;
531 #endif // LIBINT2_DERIV_ONEBODY_ORDER 533 #if LIBINT2_DERIV_ERI_ORDER 535 Matrix F2 = Matrix::Zero(atoms.size(), 3);
540 auto G1 = compute_2body_fock_deriv<1>(obs, atoms, D);
541 for (
auto atom = 0, i = 0; atom != atoms.size(); ++atom) {
542 for (
auto xyz = 0; xyz != 3; ++xyz, ++i) {
544 auto force = G1[i].cwiseProduct(D).sum();
545 F2(atom, xyz) += force;
549 std::cout <<
"** 2-body forces = ";
550 for (
int atom = 0; atom != atoms.size(); ++atom)
551 for (
int xyz = 0; xyz != 3; ++xyz) std::cout << F2(atom, xyz) <<
" ";
552 std::cout << std::endl;
557 #if LIBINT2_DERIV_ONEBODY_ORDER && LIBINT2_DERIV_ERI_ORDER 559 Matrix FN = Matrix::Zero(atoms.size(), 3);
563 for (
auto a1 = 1; a1 != atoms.size(); ++a1) {
564 const auto& atom1 = atoms[a1];
565 for (
auto a2 = 0; a2 < a1; ++a2) {
566 const auto& atom2 = atoms[a2];
568 auto x12 = atom1.x - atom2.x;
569 auto y12 = atom1.y - atom2.y;
570 auto z12 = atom1.z - atom2.z;
571 auto r12_2 = x12 * x12 + y12 * y12 + z12 * z12;
572 auto r12 = sqrt(r12_2);
573 auto r12_3 = r12 * r12_2;
575 auto z1z2_over_r12_3 =
576 atom1.atomic_number * atom2.atomic_number / r12_3;
578 auto fx = -x12 * z1z2_over_r12_3;
579 auto fy = -y12 * z1z2_over_r12_3;
580 auto fz = -z12 * z1z2_over_r12_3;
590 std::cout <<
"** nuclear repulsion forces = ";
591 for (
int atom = 0; atom != atoms.size(); ++atom)
592 for (
int xyz = 0; xyz != 3; ++xyz) std::cout << FN(atom, xyz) <<
" ";
593 std::cout << std::endl;
595 auto F = F1 + F_Pulay + F2 + FN;
596 std::cout <<
"** Hartree-Fock forces = ";
597 for (
int atom = 0; atom != atoms.size(); ++atom)
598 for (
int xyz = 0; xyz != 3; ++xyz) std::cout << F(atom, xyz) <<
" ";
599 std::cout << std::endl;
604 const auto ncoords = 3 * atoms.size();
606 const auto nelem = ncoords * (ncoords+1) / 2;
607 #if LIBINT2_DERIV_ONEBODY_ORDER > 1 609 Matrix H1 = Matrix::Zero(ncoords, ncoords);
610 Matrix H_Pulay = Matrix::Zero(ncoords, ncoords);
614 auto T2 = compute_1body_ints_deriv<Operator::kinetic>(2, obs, atoms);
615 auto V2 = compute_1body_ints_deriv<Operator::nuclear>(2, obs, atoms);
616 for (
auto row = 0, i = 0; row != ncoords; ++row) {
617 for (
auto col = row; col != ncoords; ++col, ++i) {
618 auto hess = 2 * (T2[i] + V2[i]).cwiseProduct(D).sum();
619 H1(row, col) += hess;
627 DiagonalMatrix evals_occ(evals.topRows(ndocc));
628 Matrix W = C_occ * evals_occ * C_occ.transpose();
629 auto S2 = compute_1body_ints_deriv<Operator::overlap>(2, obs, atoms);
630 for (
auto row = 0, i = 0; row != ncoords; ++row) {
631 for (
auto col = row; col != ncoords; ++col, ++i) {
632 auto hess = 2 * S2[i].cwiseProduct(W).sum();
633 H_Pulay(row, col) -= hess;
637 std::cout <<
"** 1-body hessian = ";
638 for (
auto row = 0, i = 0; row != ncoords; ++row) {
639 for (
auto col = row; col != ncoords; ++col) {
640 std::cout << H1(row, col) <<
" ";
643 std::cout << std::endl;
645 std::cout <<
"** Pulay hessian = ";
646 for (
auto row = 0, i = 0; row != ncoords; ++row) {
647 for (
auto col = row; col != ncoords; ++col) {
648 std::cout << H_Pulay(row, col) <<
" ";
651 std::cout << std::endl;
652 #endif // LIBINT2_DERIV_ONEBODY_ORDER > 1 654 #if LIBINT2_DERIV_ERI_ORDER > 1 656 Matrix H2 = Matrix::Zero(ncoords, ncoords);
661 auto G2 = compute_2body_fock_deriv<2>(obs, atoms, D);
662 for (
auto row = 0, i = 0; row != ncoords; ++row) {
663 for (
auto col = row; col != ncoords; ++col, ++i) {
665 auto hess = G2[i].cwiseProduct(D).sum();
666 H2(row, col) += hess;
670 std::cout <<
"** 2-body hessian = ";
671 for (
auto row = 0, i = 0; row != ncoords; ++row) {
672 for (
auto col = row; col != ncoords; ++col) {
673 std::cout << H2(row, col) <<
" ";
676 std::cout << std::endl;
681 #if LIBINT2_DERIV_ONEBODY_ORDER > 1 && LIBINT2_DERIV_ERI_ORDER > 1 684 Matrix HN = Matrix::Zero(ncoords, ncoords);
689 for (
auto a1 = 1; a1 != atoms.size(); ++a1) {
690 const auto& atom1 = atoms[a1];
691 for (
auto a2 = 0; a2 < a1; ++a2) {
692 const auto& atom2 = atoms[a2];
694 auto x12 = atom1.x - atom2.x;
695 auto y12 = atom1.y - atom2.y;
696 auto z12 = atom1.z - atom2.z;
697 auto x12_2 = x12 * x12;
698 auto y12_2 = y12 * y12;
699 auto z12_2 = z12 * z12;
700 auto r12_2 = x12 * x12 + y12 * y12 + z12 * z12;
701 auto r12 = sqrt(r12_2);
702 auto r12_5 = r12 * r12_2 * r12_2;
704 auto z1z2_over_r12_5 =
705 atom1.atomic_number * atom2.atomic_number / r12_5;
707 HN(3*a1 + 0, 3*a1 + 0) += z1z2_over_r12_5 * (3*x12_2 - r12_2);
708 HN(3*a1 + 1, 3*a1 + 1) += z1z2_over_r12_5 * (3*y12_2 - r12_2);
709 HN(3*a1 + 2, 3*a1 + 2) += z1z2_over_r12_5 * (3*z12_2 - r12_2);
710 HN(3*a1 + 0, 3*a1 + 1) += z1z2_over_r12_5 * (3*x12*y12);
711 HN(3*a1 + 0, 3*a1 + 2) += z1z2_over_r12_5 * (3*x12*z12);
712 HN(3*a1 + 1, 3*a1 + 2) += z1z2_over_r12_5 * (3*y12*z12);
714 HN(3*a2 + 0, 3*a2 + 0) += z1z2_over_r12_5 * (3*x12_2 - r12_2);
715 HN(3*a2 + 1, 3*a2 + 1) += z1z2_over_r12_5 * (3*y12_2 - r12_2);
716 HN(3*a2 + 2, 3*a2 + 2) += z1z2_over_r12_5 * (3*z12_2 - r12_2);
717 HN(3*a2 + 0, 3*a2 + 1) += z1z2_over_r12_5 * (3*x12*y12);
718 HN(3*a2 + 0, 3*a2 + 2) += z1z2_over_r12_5 * (3*x12*z12);
719 HN(3*a2 + 1, 3*a2 + 2) += z1z2_over_r12_5 * (3*y12*z12);
721 HN(3*a2 + 0, 3*a1 + 0) -= z1z2_over_r12_5 * (3*x12_2 - r12_2);
722 HN(3*a2 + 1, 3*a1 + 1) -= z1z2_over_r12_5 * (3*y12_2 - r12_2);
723 HN(3*a2 + 2, 3*a1 + 2) -= z1z2_over_r12_5 * (3*z12_2 - r12_2);
724 HN(3*a2 + 1, 3*a1 + 0) -= z1z2_over_r12_5 * (3*y12*x12);
725 HN(3*a2 + 2, 3*a1 + 0) -= z1z2_over_r12_5 * (3*z12*x12);
726 HN(3*a2 + 2, 3*a1 + 1) -= z1z2_over_r12_5 * (3*z12*y12);
727 HN(3*a2 + 0, 3*a1 + 1) -= z1z2_over_r12_5 * (3*x12*y12);
728 HN(3*a2 + 0, 3*a1 + 2) -= z1z2_over_r12_5 * (3*x12*z12);
729 HN(3*a2 + 1, 3*a1 + 2) -= z1z2_over_r12_5 * (3*y12*z12);
733 std::cout <<
"** nuclear repulsion hessian = ";
734 for (
auto row = 0, i = 0; row != ncoords; ++row) {
735 for (
auto col = row; col != ncoords; ++col) {
736 std::cout << HN(row, col) <<
" ";
739 std::cout << std::endl;
741 auto H = H1 + H_Pulay + H2 + HN;
742 std::cout <<
"** Hartree-Fock hessian = ";
743 for (
auto row = 0, i = 0; row != ncoords; ++row) {
744 for (
auto col = row; col != ncoords; ++col) {
745 std::cout << H(row, col) <<
" ";
748 std::cout << std::endl;
757 catch (
const char* ex) {
758 cerr <<
"caught exception: " << ex << endl;
760 }
catch (std::string& ex) {
761 cerr <<
"caught exception: " << ex << endl;
763 }
catch (std::exception& ex) {
764 cerr << ex.what() << endl;
767 cerr <<
"caught unknown exception\n";
774 std::vector<libint2::Atom> read_geometry(
const std::string& filename) {
775 std::cout <<
"Will read geometry from " << filename << std::endl;
776 std::ifstream is(filename);
778 char errmsg[256] =
"Could not open file ";
779 strncpy(errmsg + 20, filename.c_str(), 235);
781 throw std::ios_base::failure(errmsg);
788 std::ostringstream oss;
793 std::istringstream iss(oss.str());
797 if (filename.rfind(
".xyz") != std::string::npos)
800 throw "only .xyz files are accepted";
807 Matrix compute_soad(
const std::vector<libint2::Atom>& atoms) {
810 for (
const auto& atom : atoms) {
811 const auto Z = atom.atomic_number;
816 Matrix D = Matrix::Zero(nao, nao);
817 size_t ao_offset = 0;
818 for (
const auto& atom : atoms) {
819 const auto Z = atom.atomic_number;
821 for(
const auto& occ: occvec) {
822 D(ao_offset, ao_offset) = occ;
830 Matrix compute_shellblock_norm(
const BasisSet& obs,
const Matrix& A) {
831 const auto nsh = obs.size();
832 Matrix Ash(nsh, nsh);
834 auto shell2bf = obs.shell2bf();
835 for (
size_t s1 = 0; s1 != nsh; ++s1) {
836 const auto& s1_first = shell2bf[s1];
837 const auto& s1_size = obs[s1].size();
838 for (
size_t s2 = 0; s2 != nsh; ++s2) {
839 const auto& s2_first = shell2bf[s2];
840 const auto& s2_size = obs[s2].size();
842 Ash(s1, s2) = A.block(s1_first, s2_first, s1_size, s2_size)
843 .lpNorm<Eigen::Infinity>();
850 template <Operator obtype,
typename OperatorParams>
852 const BasisSet& obs, OperatorParams params) {
853 const auto n = obs.nbf();
854 const auto nshells = obs.size();
855 using libint2::nthreads;
858 const unsigned int nopers = libint2::operator_traits<obtype>::nopers;
860 for (
auto& r : result) r = Matrix::Zero(n, n);
863 std::vector<libint2::Engine> engines(nthreads);
864 engines[0] = libint2::Engine(obtype, obs.max_nprim(), obs.max_l(), 0);
868 if (obtype == Operator::nuclear || obtype == Operator::erf_nuclear ||
869 obtype == Operator::erfc_nuclear) {
870 engines[0].set_params(params);
872 for (
size_t i = 1; i != nthreads; ++i) {
873 engines[i] = engines[0];
876 auto shell2bf = obs.shell2bf();
878 auto compute = [&](
int thread_id) {
880 const auto& buf = engines[thread_id].results();
885 for (
auto s1 = 0l, s12 = 0l; s1 != nshells; ++s1) {
886 auto bf1 = shell2bf[s1];
887 auto n1 = obs[s1].size();
888 auto s1_offset = s1 * (s1+1) / 2;
890 for(
auto s2: obs_shellpair_list[s1]) {
891 auto s12 = s1_offset + s2;
892 if (s12 % nthreads != thread_id)
continue;
894 auto bf2 = shell2bf[s2];
895 auto n2 = obs[s2].size();
900 engines[thread_id].compute(obs[s1], obs[s2]);
902 for (
unsigned int op = 0; op != nopers; ++op) {
905 Eigen::Map<const Matrix> buf_mat(buf[op], n1, n2);
906 result[op].block(bf1, bf2, n1, n2) = buf_mat;
909 result[op].block(bf2, bf1, n2, n1) = buf_mat.transpose();
920 #if LIBINT2_DERIV_ONEBODY_ORDER 921 template <Operator obtype>
922 std::vector<Matrix> compute_1body_ints_deriv(
unsigned deriv_order,
924 const std::vector<libint2::Atom>& atoms) {
925 using libint2::nthreads;
926 const auto n = obs.nbf();
927 const auto nshells = obs.size();
928 constexpr
auto nopers = libint2::operator_traits<obtype>::nopers;
929 const auto nresults =
930 nopers * libint2::num_geometrical_derivatives(atoms.size(), deriv_order);
931 typedef std::vector<Matrix> result_type;
932 result_type result(nresults);
933 for (
auto& r : result) r = Matrix::Zero(n, n);
936 std::vector<libint2::Engine> engines(nthreads);
938 libint2::Engine(obtype, obs.max_nprim(), obs.max_l(), deriv_order);
942 if (obtype == Operator::nuclear) {
943 std::vector<std::pair<double, std::array<double, 3>>> q;
944 for (
const auto& atom : atoms) {
945 q.push_back({static_cast<double>(atom.atomic_number),
946 {{atom.x, atom.y, atom.z}}});
948 engines[0].set_params(q);
950 for (
size_t i = 1; i != nthreads; ++i) {
951 engines[i] = engines[0];
954 auto shell2bf = obs.shell2bf();
955 auto shell2atom = obs.shell2atom(atoms);
957 const auto natoms = atoms.size();
958 const auto two_times_ncoords = 6*natoms;
959 const auto nderivcenters_shset =
960 2 + ((obtype == Operator::nuclear) ? natoms : 0);
962 auto compute = [&](
int thread_id) {
964 const auto& buf = engines[thread_id].results();
969 for (
auto s1 = 0l, s12 = 0l; s1 != nshells; ++s1) {
970 auto bf1 = shell2bf[s1];
971 auto n1 = obs[s1].size();
972 auto atom1 = shell2atom[s1];
975 auto s1_offset = s1 * (s1+1) / 2;
976 for(
auto s2: obs_shellpair_list[s1]) {
977 auto s12 = s1_offset + s2;
978 if (s12 % nthreads != thread_id)
continue;
980 auto bf2 = shell2bf[s2];
981 auto n2 = obs[s2].size();
982 auto atom2 = shell2atom[s2];
987 engines[thread_id].compute(obs[s1], obs[s2]);
991 auto add_shellset_to_dest = [&](std::size_t op, std::size_t idx,
992 double scale = 1.0) {
995 Eigen::Map<const Matrix> buf_mat(buf[idx], n1, n2);
997 result[op].block(bf1, bf2, n1, n2) += buf_mat;
999 result[op].block(bf1, bf2, n1, n2) += scale * buf_mat;
1003 result[op].block(bf2, bf1, n2, n1) += buf_mat.transpose();
1005 result[op].block(bf2, bf1, n2, n1) += scale * buf_mat.transpose();
1009 switch (deriv_order) {
1011 for (std::size_t op = 0; op != nopers; ++op) {
1012 add_shellset_to_dest(op, op);
1034 std::size_t shellset_idx = 0;
1035 for (
auto c = 0; c != nderivcenters_shset; ++c) {
1036 auto atom = (c == 0) ? atom1 : ((c == 1) ? atom2 : c - 2);
1037 auto op_start = 3 * atom * nopers;
1038 auto op_fence = op_start + nopers;
1039 for (
auto xyz = 0; xyz != 3;
1040 ++xyz, op_start += nopers, op_fence += nopers) {
1041 for (
unsigned int op = op_start; op != op_fence;
1042 ++op, ++shellset_idx) {
1043 add_shellset_to_dest(op, shellset_idx);
1068 #define upper_triangle_index(n2, i, j) \ 1069 (std::min((i), (j))) * ((n2) - (std::min((i), (j))) - 1) / 2 + \ 1070 (std::max((i), (j))) 1073 std::size_t shellset_idx = 0;
1074 for (
auto c1 = 0; c1 != nderivcenters_shset; ++c1) {
1075 auto a1 = (c1 == 0) ? atom1 : ((c1 == 1) ? atom2 : c1 - 2);
1076 auto coord1 = 3 * a1;
1077 for (
auto xyz1 = 0; xyz1 != 3; ++xyz1, ++coord1) {
1079 for (
auto c2 = c1; c2 != nderivcenters_shset; ++c2) {
1080 auto a2 = (c2 == 0) ? atom1 : ((c2 == 1) ? atom2 : c2 - 2);
1081 auto xyz2_start = (c1 == c2) ? xyz1 : 0;
1082 auto coord2 = 3 * a2 + xyz2_start;
1083 for (
auto xyz2 = xyz2_start; xyz2 != 3;
1086 double scale = (coord1 == coord2 && c1 != c2) ? 2.0 : 1.0;
1088 const auto coord12 =
1089 upper_triangle_index(two_times_ncoords, coord1, coord2);
1090 auto op_start = coord12 * nopers;
1091 auto op_fence = op_start + nopers;
1092 for (
auto op = op_start; op != op_fence;
1093 ++op, ++shellset_idx) {
1094 add_shellset_to_dest(op, shellset_idx, scale);
1101 #undef upper_triangle_index 1104 assert(
false &&
"not yet implemented");
1106 using ShellSetDerivIterator =
1108 std::vector<unsigned int>>;
1109 ShellSetDerivIterator shellset_diter(deriv_order,
1110 nderivcenters_shset);
1111 while (shellset_diter) {
1112 const auto& deriv = *shellset_diter;
1127 template <lib
int2::Operator Kernel>
1128 Matrix compute_schwarz_ints(
1129 const BasisSet& bs1,
const BasisSet& _bs2,
bool use_2norm,
1130 typename libint2::operator_traits<Kernel>::oper_params_type params) {
1131 const BasisSet& bs2 = (_bs2.empty() ? bs1 : _bs2);
1132 const auto nsh1 = bs1.size();
1133 const auto nsh2 = bs2.size();
1134 const auto bs1_equiv_bs2 = (&bs1 == &bs2);
1136 Matrix K = Matrix::Zero(nsh1, nsh2);
1139 using libint2::Engine;
1140 using libint2::nthreads;
1141 std::vector<Engine> engines(nthreads);
1145 engines[0] = Engine(Kernel, std::max(bs1.max_nprim(), bs2.max_nprim()),
1146 std::max(bs1.max_l(), bs2.max_l()), 0, epsilon, params);
1147 for (
size_t i = 1; i != nthreads; ++i) {
1148 engines[i] = engines[0];
1151 std::cout <<
"computing Schwarz bound prerequisites (kernel=" << (int)Kernel
1158 auto compute = [&](
int thread_id) {
1160 const auto& buf = engines[thread_id].results();
1163 for (
auto s1 = 0l, s12 = 0l; s1 != nsh1; ++s1) {
1164 auto n1 = bs1[s1].size();
1166 auto s2_max = bs1_equiv_bs2 ? s1 : nsh2 - 1;
1167 for (
auto s2 = 0; s2 <= s2_max; ++s2, ++s12) {
1168 if (s12 % nthreads != thread_id)
continue;
1170 auto n2 = bs2[s2].size();
1173 engines[thread_id].compute2<Kernel, BraKet::xx_xx, 0>(bs1[s1], bs2[s2],
1175 assert(buf[0] !=
nullptr &&
1176 "to compute Schwarz ints turn off primitive screening");
1179 Eigen::Map<const Matrix> buf_mat(buf[0], n12, n12);
1180 auto norm2 = use_2norm ? buf_mat.diagonal().norm()
1181 : buf_mat.diagonal().lpNorm<Eigen::Infinity>();
1182 K(s1, s2) = std::sqrt(norm2);
1183 if (bs1_equiv_bs2) K(s2, s1) = K(s1, s2);
1191 std::cout <<
"done (" << timer.read(0) <<
" s)" << std::endl;
1196 Matrix compute_do_ints(
const BasisSet& bs1,
const BasisSet& bs2,
1198 return compute_schwarz_ints<libint2::Operator::delta>(bs1, bs2, use_2norm);
1201 shellpair_list_t compute_shellpair_list(
const BasisSet& bs1,
1202 const BasisSet& _bs2,
1203 const double threshold) {
1204 const BasisSet& bs2 = (_bs2.empty() ? bs1 : _bs2);
1205 const auto nsh1 = bs1.size();
1206 const auto nsh2 = bs2.size();
1207 const auto bs1_equiv_bs2 = (&bs1 == &bs2);
1209 using libint2::nthreads;
1212 using libint2::Engine;
1213 std::vector<Engine> engines;
1214 engines.reserve(nthreads);
1215 engines.emplace_back(Operator::overlap,
1216 std::max(bs1.max_nprim(), bs2.max_nprim()),
1217 std::max(bs1.max_l(), bs2.max_l()), 0);
1218 for (
size_t i = 1; i != nthreads; ++i) {
1219 engines.push_back(engines[0]);
1222 std::cout <<
"computing non-negligible shell-pair list ... ";
1228 shellpair_list_t result;
1232 auto compute = [&](
int thread_id) {
1234 auto& engine = engines[thread_id];
1235 const auto& buf = engine.results();
1238 for (
auto s1 = 0l, s12 = 0l; s1 != nsh1; ++s1) {
1240 if (result.find(s1) == result.end())
1241 result.insert(std::make_pair(s1, std::vector<size_t>()));
1244 auto n1 = bs1[s1].size();
1246 auto s2_max = bs1_equiv_bs2 ? s1 : nsh2 - 1;
1247 for (
auto s2 = 0; s2 <= s2_max; ++s2, ++s12) {
1248 if (s12 % nthreads != thread_id)
continue;
1250 auto on_same_center = (bs1[s1].O == bs2[s2].O);
1251 bool significant = on_same_center;
1252 if (not on_same_center) {
1253 auto n2 = bs2[s2].size();
1254 engines[thread_id].compute(bs1[s1], bs2[s2]);
1255 Eigen::Map<const Matrix> buf_mat(buf[0], n1, n2);
1256 auto norm = buf_mat.norm();
1257 significant = (norm >= threshold);
1262 result[s1].emplace_back(s2);
1273 auto sort = [&](
int thread_id) {
1274 for (
auto s1 = 0l; s1 != nsh1; ++s1) {
1275 if (s1 % nthreads == thread_id) {
1276 auto& list = result[s1];
1277 std::sort(list.begin(), list.end());
1285 std::cout <<
"done (" << timer.
read(0) <<
" s)" << std::endl;
1302 std::tuple<Matrix, Matrix, size_t, double, double> gensqrtinv(
1303 const Matrix& S,
bool symmetric =
false,
1304 double max_condition_number = 1e8) {
1305 Eigen::SelfAdjointEigenSolver<Matrix> eig_solver(S);
1306 auto U = eig_solver.eigenvectors();
1307 auto s = eig_solver.eigenvalues();
1308 auto s_max = s.maxCoeff();
1309 auto condition_number = std::min(
1310 s_max / std::max(s.minCoeff(), std::numeric_limits<double>::min()),
1311 1.0 / std::numeric_limits<double>::epsilon());
1312 auto threshold = s_max / max_condition_number;
1315 for (
long i = n - 1; i >= 0; --i) {
1316 if (s(i) >= threshold) {
1322 auto sigma = s.bottomRows(n_cond);
1323 auto result_condition_number = sigma.maxCoeff() / sigma.minCoeff();
1324 auto sigma_sqrt = sigma.array().sqrt().matrix().asDiagonal();
1325 auto sigma_invsqrt = sigma.array().sqrt().inverse().matrix().asDiagonal();
1328 auto U_cond = U.block(0, n - n_cond, n, n_cond);
1329 Matrix X = U_cond * sigma_invsqrt;
1330 Matrix Xinv = U_cond * sigma_sqrt;
1333 X = X * U_cond.transpose();
1334 Xinv = Xinv * U_cond.transpose();
1336 return std::make_tuple(X, Xinv,
size_t(n_cond), condition_number,
1337 result_condition_number);
1340 std::tuple<Matrix, Matrix, double> conditioning_orthogonalizer(
1341 const Matrix& S,
double S_condition_number_threshold) {
1343 double S_condition_number;
1344 double XtX_condition_number;
1347 assert(S.rows() == S.cols());
1349 std::tie(X, Xinv, obs_rank, S_condition_number, XtX_condition_number) =
1350 gensqrtinv(S,
false, S_condition_number_threshold);
1351 auto obs_nbf_omitted = (long)S.rows() - (long)obs_rank;
1352 std::cout <<
"overlap condition number = " << S_condition_number;
1353 if (obs_nbf_omitted > 0)
1354 std::cout <<
" (dropped " << obs_nbf_omitted <<
" " 1355 << (obs_nbf_omitted > 1 ?
"fns" :
"fn") <<
" to reduce to " 1356 << XtX_condition_number <<
")";
1357 std::cout << std::endl;
1359 if (obs_nbf_omitted > 0) {
1360 Matrix should_be_I = X.transpose() * S * X;
1361 Matrix I = Matrix::Identity(should_be_I.rows(), should_be_I.cols());
1362 std::cout <<
"||X^t * S * X - I||_2 = " << (should_be_I - I).norm()
1363 <<
" (should be 0)" << std::endl;
1366 return std::make_tuple(X, Xinv, XtX_condition_number);
1369 Matrix compute_2body_2index_ints(
const BasisSet& bs) {
1370 using libint2::nthreads;
1371 const auto n = bs.nbf();
1372 const auto nshells = bs.size();
1373 Matrix result = Matrix::Zero(n, n);
1376 using libint2::Engine;
1377 std::vector<Engine> engines(nthreads);
1379 Engine(libint2::Operator::coulomb, bs.max_nprim(), bs.max_l(), 0);
1380 engines[0].set(BraKet::xs_xs);
1381 for (
size_t i = 1; i != nthreads; ++i) {
1382 engines[i] = engines[0];
1385 auto shell2bf = bs.shell2bf();
1386 auto unitshell = Shell::unit();
1388 auto compute = [&](
int thread_id) {
1390 auto& engine = engines[thread_id];
1391 const auto& buf = engine.results();
1396 for (
auto s1 = 0l, s12 = 0l; s1 != nshells; ++s1) {
1397 auto bf1 = shell2bf[s1];
1398 auto n1 = bs[s1].size();
1400 for (
auto s2 = 0; s2 <= s1; ++s2, ++s12) {
1401 if (s12 % nthreads != thread_id)
continue;
1403 auto bf2 = shell2bf[s2];
1404 auto n2 = bs[s2].size();
1407 engine.compute(bs[s1], bs[s2]);
1408 if (buf[0] ==
nullptr)
1413 Eigen::Map<const Matrix> buf_mat(buf[0], n1, n2);
1414 result.block(bf1, bf2, n1, n2) = buf_mat;
1417 result.block(bf2, bf1, n2, n1) = buf_mat.transpose();
1427 Matrix compute_2body_fock(
const BasisSet& obs,
const Matrix& D,
1428 double precision,
const Matrix& Schwarz) {
1429 const auto n = obs.nbf();
1430 const auto nshells = obs.size();
1431 using libint2::nthreads;
1432 std::vector<Matrix> G(nthreads, Matrix::Zero(n, n));
1434 const auto do_schwarz_screen = Schwarz.cols() != 0 && Schwarz.rows() != 0;
1435 Matrix D_shblk_norm =
1436 compute_shellblock_norm(obs, D);
1438 auto fock_precision = precision;
1441 auto max_nprim = obs.max_nprim();
1442 auto max_nprim4 = max_nprim * max_nprim * max_nprim * max_nprim;
1443 auto engine_precision = std::min(fock_precision / D_shblk_norm.maxCoeff(),
1444 std::numeric_limits<double>::epsilon()) /
1448 using libint2::Engine;
1449 std::vector<Engine> engines(nthreads);
1450 engines[0] = Engine(Operator::coulomb, obs.max_nprim(), obs.max_l(), 0);
1451 engines[0].set_precision(engine_precision);
1455 std::cout <<
"compute_2body_fock:precision = " << precision << std::endl;
1456 std::cout <<
"Engine::precision = " << engines[0].precision() << std::endl;
1457 for (
size_t i = 1; i != nthreads; ++i) {
1458 engines[i] = engines[0];
1460 std::atomic<size_t> num_ints_computed{0};
1462 #if defined(REPORT_INTEGRAL_TIMINGS) 1463 std::vector<libint2::Timers<1>> timers(nthreads);
1466 auto shell2bf = obs.shell2bf();
1468 auto lambda = [&](
int thread_id) {
1470 auto& engine = engines[thread_id];
1471 auto& g = G[thread_id];
1472 const auto& buf = engine.results();
1474 #if defined(REPORT_INTEGRAL_TIMINGS) 1475 auto& timer = timers[thread_id];
1481 for (
auto s1 = 0l, s1234 = 0l; s1 != nshells; ++s1) {
1482 auto bf1_first = shell2bf[s1];
1483 auto n1 = obs[s1].size();
1485 for (
const auto& s2 : obs_shellpair_list[s1]) {
1486 auto bf2_first = shell2bf[s2];
1487 auto n2 = obs[s2].size();
1489 const auto Dnorm12 = do_schwarz_screen ? D_shblk_norm(s1, s2) : 0.;
1491 for (
auto s3 = 0; s3 <= s1; ++s3) {
1492 auto bf3_first = shell2bf[s3];
1493 auto n3 = obs[s3].size();
1495 const auto Dnorm123 =
1497 ? std::max(D_shblk_norm(s1, s3),
1498 std::max(D_shblk_norm(s2, s3), Dnorm12))
1501 const auto s4_max = (s1 == s3) ? s2 : s3;
1502 for (
const auto& s4 : obs_shellpair_list[s3]) {
1507 if ((s1234++) % nthreads != thread_id)
continue;
1509 const auto Dnorm1234 =
1512 D_shblk_norm(s1, s4),
1513 std::max(D_shblk_norm(s2, s4),
1514 std::max(D_shblk_norm(s3, s4), Dnorm123)))
1517 if (do_schwarz_screen &&
1518 Dnorm1234 * Schwarz(s1, s2) * Schwarz(s3, s4) <
1522 auto bf4_first = shell2bf[s4];
1523 auto n4 = obs[s4].size();
1525 num_ints_computed += n1 * n2 * n3 * n4;
1529 auto s12_deg = (s1 == s2) ? 1.0 : 2.0;
1530 auto s34_deg = (s3 == s4) ? 1.0 : 2.0;
1531 auto s12_34_deg = (s1 == s3) ? (s2 == s4 ? 1.0 : 2.0) : 2.0;
1532 auto s1234_deg = s12_deg * s34_deg * s12_34_deg;
1534 #if defined(REPORT_INTEGRAL_TIMINGS) 1538 engine.compute2<Operator::coulomb, BraKet::xx_xx, 0>(
1539 obs[s1], obs[s2], obs[s3], obs[s4]);
1540 const auto* buf_1234 = buf[0];
1541 if (buf_1234 ==
nullptr)
1544 #if defined(REPORT_INTEGRAL_TIMINGS) 1548 Eigen::Map<MatrixXd> buf_1234_map(buf_1234, n12, n34);
1549 assert(buf_1234_map.norm() < Schwarz(s1, s2) * Schwarz(s3, s4));
1563 for (
auto f1 = 0, f1234 = 0; f1 != n1; ++f1) {
1564 const auto bf1 = f1 + bf1_first;
1565 for (
auto f2 = 0; f2 != n2; ++f2) {
1566 const auto bf2 = f2 + bf2_first;
1567 for (
auto f3 = 0; f3 != n3; ++f3) {
1568 const auto bf3 = f3 + bf3_first;
1569 for (
auto f4 = 0; f4 != n4; ++f4, ++f1234) {
1570 const auto bf4 = f4 + bf4_first;
1572 const auto value = buf_1234[f1234];
1574 const auto value_scal_by_deg = value * s1234_deg;
1576 g(bf1, bf2) += D(bf3, bf4) * value_scal_by_deg;
1577 g(bf3, bf4) += D(bf1, bf2) * value_scal_by_deg;
1578 g(bf1, bf3) -= 0.25 * D(bf2, bf4) * value_scal_by_deg;
1579 g(bf2, bf4) -= 0.25 * D(bf1, bf3) * value_scal_by_deg;
1580 g(bf1, bf4) -= 0.25 * D(bf2, bf3) * value_scal_by_deg;
1581 g(bf2, bf3) -= 0.25 * D(bf1, bf4) * value_scal_by_deg;
1596 for (
size_t i = 1; i != nthreads; ++i) {
1600 #if defined(REPORT_INTEGRAL_TIMINGS) 1601 double time_for_ints = 0.0;
1602 for (
auto& t : timers) {
1603 time_for_ints += t.read(0);
1605 std::cout <<
"time for integrals = " << time_for_ints << std::endl;
1606 for (
int t = 0; t != nthreads; ++t) engines[t].print_timers();
1609 Matrix GG = 0.5 * (G[0] + G[0].transpose());
1611 std::cout <<
"# of integrals = " << num_ints_computed << std::endl;
1617 #if LIBINT2_DERIV_ERI_ORDER 1618 template <
unsigned deriv_order>
1619 std::vector<Matrix> compute_2body_fock_deriv(
const BasisSet& obs,
1620 const std::vector<libint2::Atom>& atoms,
1621 const Matrix& D,
double precision,
1622 const Matrix& Schwarz) {
1623 const auto n = obs.nbf();
1624 const auto nshells = obs.size();
1625 const auto nderiv_shellset =
1626 libint2::num_geometrical_derivatives(4, deriv_order);
1627 const auto nderiv = libint2::num_geometrical_derivatives(
1628 atoms.size(), deriv_order);
1629 const auto ncoords_times_two = (atoms.size() * 3) * 2;
1630 using libint2::nthreads;
1631 std::vector<Matrix> G(nthreads * nderiv, Matrix::Zero(n, n));
1633 const auto do_schwarz_screen = Schwarz.cols() != 0 && Schwarz.rows() != 0;
1634 Matrix D_shblk_norm =
1635 compute_shellblock_norm(obs, D);
1637 auto fock_precision = precision;
1640 auto max_nprim = obs.max_nprim();
1641 auto max_nprim4 = max_nprim * max_nprim * max_nprim * max_nprim;
1642 auto engine_precision = std::min(fock_precision / D_shblk_norm.maxCoeff(),
1643 std::numeric_limits<double>::epsilon()) /
1647 using libint2::Engine;
1648 std::vector<Engine> engines(nthreads);
1650 Engine(Operator::coulomb, obs.max_nprim(), obs.max_l(), deriv_order);
1651 engines[0].set_precision(engine_precision);
1655 std::cout <<
"compute_2body_fock:precision = " << precision << std::endl;
1656 std::cout <<
"Engine::precision = " << engines[0].precision() << std::endl;
1657 for (
size_t i = 1; i != nthreads; ++i) {
1658 engines[i] = engines[0];
1660 std::atomic<size_t> num_ints_computed{0};
1662 #if defined(REPORT_INTEGRAL_TIMINGS) 1663 std::vector<libint2::Timers<1>> timers(nthreads);
1666 auto shell2bf = obs.shell2bf();
1667 auto shell2atom = obs.shell2atom(atoms);
1669 auto lambda = [&](
int thread_id) {
1671 auto& engine = engines[thread_id];
1672 const auto& buf = engine.results();
1674 #if defined(REPORT_INTEGRAL_TIMINGS) 1675 auto& timer = timers[thread_id];
1680 size_t shell_atoms[4];
1683 for (
auto s1 = 0l, s1234 = 0l; s1 != nshells; ++s1) {
1684 auto bf1_first = shell2bf[s1];
1685 auto n1 = obs[s1].size();
1686 shell_atoms[0] = shell2atom[s1];
1688 for (
const auto& s2 : obs_shellpair_list[s1]) {
1689 auto bf2_first = shell2bf[s2];
1690 auto n2 = obs[s2].size();
1691 shell_atoms[1] = shell2atom[s2];
1693 const auto Dnorm12 = do_schwarz_screen ? D_shblk_norm(s1, s2) : 0.;
1695 for (
auto s3 = 0; s3 <= s1; ++s3) {
1696 auto bf3_first = shell2bf[s3];
1697 auto n3 = obs[s3].size();
1698 shell_atoms[2] = shell2atom[s3];
1700 const auto Dnorm123 =
1702 ? std::max(D_shblk_norm(s1, s3),
1703 std::max(D_shblk_norm(s2, s3), Dnorm12))
1706 const auto s4_max = (s1 == s3) ? s2 : s3;
1707 for (
const auto& s4 : obs_shellpair_list[s3]) {
1712 if ((s1234++) % nthreads != thread_id)
continue;
1714 const auto Dnorm1234 =
1717 D_shblk_norm(s1, s4),
1718 std::max(D_shblk_norm(s2, s4),
1719 std::max(D_shblk_norm(s3, s4), Dnorm123)))
1722 if (do_schwarz_screen &&
1723 Dnorm1234 * Schwarz(s1, s2) * Schwarz(s3, s4) <
1727 auto bf4_first = shell2bf[s4];
1728 auto n4 = obs[s4].size();
1729 shell_atoms[3] = shell2atom[s4];
1731 const auto n1234 = n1 * n2 * n3 * n4;
1735 auto s12_deg = (s1 == s2) ? 1.0 : 2.0;
1736 auto s34_deg = (s3 == s4) ? 1.0 : 2.0;
1737 auto s12_34_deg = (s1 == s3) ? (s2 == s4 ? 1.0 : 2.0) : 2.0;
1738 auto s1234_deg = s12_deg * s34_deg * s12_34_deg;
1742 auto add_shellset_to_dest = [&](
1743 std::size_t op, std::size_t idx,
int coord1,
int coord2,
double scale = 1.0) {
1745 auto shset = buf[idx];
1746 const auto weight = scale * s1234_deg;
1748 for (
auto f1 = 0, f1234 = 0; f1 != n1; ++f1) {
1749 const auto bf1 = f1 + bf1_first;
1750 for (
auto f2 = 0; f2 != n2; ++f2) {
1751 const auto bf2 = f2 + bf2_first;
1752 for (
auto f3 = 0; f3 != n3; ++f3) {
1753 const auto bf3 = f3 + bf3_first;
1754 for (
auto f4 = 0; f4 != n4; ++f4, ++f1234) {
1755 const auto bf4 = f4 + bf4_first;
1757 const auto value = shset[f1234];
1758 const auto wvalue = value * weight;
1760 g(bf1, bf2) += D(bf3, bf4) * wvalue;
1761 g(bf3, bf4) += D(bf1, bf2) * wvalue;
1762 g(bf1, bf3) -= 0.25 * D(bf2, bf4) * wvalue;
1763 g(bf2, bf4) -= 0.25 * D(bf1, bf3) * wvalue;
1764 g(bf1, bf4) -= 0.25 * D(bf2, bf3) * wvalue;
1765 g(bf2, bf3) -= 0.25 * D(bf1, bf4) * wvalue;
1772 #if defined(REPORT_INTEGRAL_TIMINGS) 1776 engine.compute2<Operator::coulomb, BraKet::xx_xx, deriv_order>(
1777 obs[s1], obs[s2], obs[s3], obs[s4]);
1778 if (buf[0] ==
nullptr)
1780 num_ints_computed += nderiv_shellset * n1234;
1782 #if defined(REPORT_INTEGRAL_TIMINGS) 1786 switch (deriv_order) {
1788 int coord1 = 0, coord2 = 0;
1789 add_shellset_to_dest(thread_id, 0, coord1, coord2);
1793 for (
auto d = 0; d != 12; ++d) {
1794 const int a = d / 3;
1795 const int xyz = d % 3;
1797 auto coord = shell_atoms[a] * 3 + xyz;
1798 auto& g = G[thread_id * nderiv + coord];
1800 int coord1 = 0, coord2 = 0;
1802 add_shellset_to_dest(thread_id * nderiv + coord, d, coord1, coord2);
1811 #define upper_triangle_index(n2, i, j) \ 1812 (std::min((i), (j))) * ((n2) - (std::min((i), (j))) - 1) / 2 + \ 1813 (std::max((i), (j))) 1815 std::size_t shellset_idx = 0;
1816 for (
auto c1 = 0; c1 != 4; ++c1) {
1817 auto a1 = shell_atoms[c1];
1818 auto coord1 = 3 * a1;
1819 for (
auto xyz1 = 0; xyz1 != 3; ++xyz1, ++coord1) {
1820 for (
auto c2 = c1; c2 != 4; ++c2) {
1821 auto a2 = shell_atoms[c2];
1822 auto xyz2_start = (c1 == c2) ? xyz1 : 0;
1823 auto coord2 = 3 * a2 + xyz2_start;
1824 for (
auto xyz2 = xyz2_start; xyz2 != 3;
1827 (coord1 == coord2 && c1 != c2) ? 2.0 : 1.0;
1829 const auto coord12 = upper_triangle_index(
1830 ncoords_times_two, coord1, coord2);
1831 auto op = thread_id * nderiv + coord12;
1832 add_shellset_to_dest(op, shellset_idx, coord1, coord2, scale);
1839 #undef upper_triangle_index 1842 assert(deriv_order <= 2 &&
1843 "support for 3rd and higher derivatives of the Fock " 1844 "matrix not yet implemented");
1856 for (
size_t t = 1; t != nthreads; ++t) {
1857 for (
auto d = 0; d != nderiv; ++d) {
1858 G[d] += G[t * nderiv + d];
1862 #if defined(REPORT_INTEGRAL_TIMINGS) 1863 double time_for_ints = 0.0;
1864 for (
auto& t : timers) {
1865 time_for_ints += t.read(0);
1867 std::cout <<
"time for integrals = " << time_for_ints << std::endl;
1868 for (
int t = 0; t != nthreads; ++t) engines[t].print_timers();
1871 std::vector<Matrix> GG(nderiv);
1872 for (
auto d = 0; d != nderiv; ++d) {
1873 GG[d] = 0.5 * (G[d] + G[d].transpose());
1876 std::cout <<
"# of integrals = " << num_ints_computed << std::endl;
1884 Matrix compute_2body_fock_general(
const BasisSet& obs,
const Matrix& D,
1885 const BasisSet& D_bs,
bool D_is_shelldiagonal,
1887 const auto n = obs.nbf();
1888 const auto nshells = obs.size();
1889 const auto n_D = D_bs.nbf();
1890 assert(D.cols() == D.rows() && D.cols() == n_D);
1892 using libint2::nthreads;
1893 std::vector<Matrix> G(nthreads, Matrix::Zero(n, n));
1896 using libint2::Engine;
1897 std::vector<Engine> engines(nthreads);
1898 engines[0] = Engine(libint2::Operator::coulomb,
1899 std::max(obs.max_nprim(), D_bs.max_nprim()),
1900 std::max(obs.max_l(), D_bs.max_l()), 0);
1901 engines[0].set_precision(precision);
1905 for (
size_t i = 1; i != nthreads; ++i) {
1906 engines[i] = engines[0];
1908 auto shell2bf = obs.shell2bf();
1909 auto shell2bf_D = D_bs.shell2bf();
1911 auto lambda = [&](
int thread_id) {
1913 auto& engine = engines[thread_id];
1914 auto& g = G[thread_id];
1915 const auto& buf = engine.results();
1918 for (
auto s1 = 0l, s1234 = 0l; s1 != nshells; ++s1) {
1919 auto bf1_first = shell2bf[s1];
1920 auto n1 = obs[s1].size();
1922 for (
auto s2 = 0; s2 <= s1; ++s2) {
1923 auto bf2_first = shell2bf[s2];
1924 auto n2 = obs[s2].size();
1926 for (
auto s3 = 0; s3 < D_bs.size(); ++s3) {
1927 auto bf3_first = shell2bf_D[s3];
1928 auto n3 = D_bs[s3].size();
1930 auto s4_begin = D_is_shelldiagonal ? s3 : 0;
1931 auto s4_fence = D_is_shelldiagonal ? s3 + 1 : D_bs.size();
1933 for (
auto s4 = s4_begin; s4 != s4_fence; ++s4, ++s1234) {
1934 if (s1234 % nthreads != thread_id)
continue;
1936 auto bf4_first = shell2bf_D[s4];
1937 auto n4 = D_bs[s4].size();
1941 auto s12_deg = (s1 == s2) ? 1.0 : 2.0;
1944 auto s34_deg = (s3 == s4) ? 1.0 : 2.0;
1945 auto s1234_deg = s12_deg * s34_deg;
1947 engine.compute2<Operator::coulomb, BraKet::xx_xx, 0>(
1948 obs[s1], obs[s2], D_bs[s3], D_bs[s4]);
1949 const auto* buf_1234 = buf[0];
1950 if (buf_1234 !=
nullptr) {
1951 for (
auto f1 = 0, f1234 = 0; f1 != n1; ++f1) {
1952 const auto bf1 = f1 + bf1_first;
1953 for (
auto f2 = 0; f2 != n2; ++f2) {
1954 const auto bf2 = f2 + bf2_first;
1955 for (
auto f3 = 0; f3 != n3; ++f3) {
1956 const auto bf3 = f3 + bf3_first;
1957 for (
auto f4 = 0; f4 != n4; ++f4, ++f1234) {
1958 const auto bf4 = f4 + bf4_first;
1960 const auto value = buf_1234[f1234];
1961 const auto value_scal_by_deg = value * s1234_deg;
1962 g(bf1, bf2) += 2.0 * D(bf3, bf4) * value_scal_by_deg;
1970 engine.compute2<Operator::coulomb, BraKet::xx_xx, 0>(
1971 obs[s1], D_bs[s3], obs[s2], D_bs[s4]);
1972 const auto* buf_1324 = buf[0];
1973 if (buf_1324 ==
nullptr)
1976 for (
auto f1 = 0, f1324 = 0; f1 != n1; ++f1) {
1977 const auto bf1 = f1 + bf1_first;
1978 for (
auto f3 = 0; f3 != n3; ++f3) {
1979 const auto bf3 = f3 + bf3_first;
1980 for (
auto f2 = 0; f2 != n2; ++f2) {
1981 const auto bf2 = f2 + bf2_first;
1982 for (
auto f4 = 0; f4 != n4; ++f4, ++f1324) {
1983 const auto bf4 = f4 + bf4_first;
1985 const auto value = buf_1324[f1324];
1986 const auto value_scal_by_deg = value * s12_deg;
1987 g(bf1, bf2) -= D(bf3, bf4) * value_scal_by_deg;
2002 for (
size_t i = 1; i != nthreads; ++i) {
2007 return 0.5 * (G[0] + G[0].transpose());
2010 #ifdef HAVE_DENSITY_FITTING 2012 Matrix DFFockEngine::compute_2body_fock_dfC(
const Matrix& Cocc) {
2014 using libint2::nthreads;
2016 const auto n = obs.nbf();
2017 const auto ndf = dfbs.nbf();
2021 std::vector<libint2::Timers<5>> timers(nthreads);
2024 typedef btas::RangeNd<CblasRowMajor, std::array<long, 1>> Range1d;
2025 typedef btas::RangeNd<CblasRowMajor, std::array<long, 2>> Range2d;
2026 typedef btas::Tensor<double, Range1d> Tensor1d;
2027 typedef btas::Tensor<double, Range2d> Tensor2d;
2031 if (xyK.size() == 0) {
2033 wall_timer.
start(0);
2035 const auto nshells = obs.size();
2036 const auto nshells_df = dfbs.size();
2042 std::vector<libint2::Engine> engines(nthreads);
2043 engines[0] = libint2::Engine(libint2::Operator::coulomb,
2044 std::max(obs.max_nprim(), dfbs.max_nprim()),
2045 std::max(obs.max_l(), dfbs.max_l()), 0);
2046 engines[0].set(BraKet::xs_xx);
2047 for (
size_t i = 1; i != nthreads; ++i) {
2048 engines[i] = engines[0];
2051 auto shell2bf = obs.shell2bf();
2052 auto shell2bf_df = dfbs.shell2bf();
2054 Tensor3d Zxy{ndf, n, n};
2056 auto lambda = [&](
int thread_id) {
2058 auto& engine = engines[thread_id];
2059 auto& timer = timers[thread_id];
2060 const auto& results = engine.results();
2063 for (
auto s1 = 0l, s123 = 0l; s1 != nshells_df; ++s1) {
2064 auto bf1_first = shell2bf_df[s1];
2065 auto n1 = dfbs[s1].size();
2067 for (
auto s2 = 0; s2 != nshells; ++s2) {
2068 auto bf2_first = shell2bf[s2];
2069 auto n2 = obs[s2].size();
2070 const auto n12 = n1 * n2;
2072 for (
auto s3 = 0; s3 != nshells; ++s3, ++s123) {
2073 if (s123 % nthreads != thread_id)
continue;
2075 auto bf3_first = shell2bf[s3];
2076 auto n3 = obs[s3].size();
2077 const auto n123 = n12 * n3;
2081 engine.compute2<Operator::coulomb, BraKet::xs_xx, 0>(
2082 dfbs[s1], unitshell, obs[s2], obs[s3]);
2083 const auto* buf = results[0];
2090 auto lower_bound = {bf1_first, bf2_first, bf3_first};
2091 auto upper_bound = {bf1_first + n1, bf2_first + n2, bf3_first + n3};
2092 auto view = btas::make_view(
2093 Zxy.range().slice(lower_bound, upper_bound), Zxy.storage());
2094 std::copy(buf, buf + n123, view.begin());
2107 double ints_time = 0;
2108 for(
const auto& timer: timers) ints_time += timer.
read(0);
2109 std::cout <<
"time for Zxy integrals = " << ints_time <<
" (total from all threads)" << std::endl;
2110 double copy_time = 0;
2111 for(
const auto& timer: timers) copy_time += timer.
read(1);
2112 std::cout <<
"time for copying into BTAS = " << copy_time <<
" (total from all threads)"<< std::endl;
2113 std::cout <<
"wall time for Zxy integrals + copy = " << wall_timer.
read(0) << std::endl;
2117 Matrix V = compute_2body_2index_ints(dfbs);
2118 Eigen::LLT<Matrix> V_LLt(V);
2119 Matrix I = Matrix::Identity(ndf, ndf);
2120 auto L = V_LLt.matrixL();
2122 Matrix Linv = L.solve(I).transpose();
2131 Tensor2d K{ndf, ndf};
2132 std::copy(Linv.data(), Linv.data() + ndf * ndf, K.begin());
2134 xyK = Tensor3d{n, n, ndf};
2135 btas::contract(1.0, Zxy, {1, 2, 3}, K, {1, 4}, 0.0, xyK, {2, 3, 4});
2136 Zxy = Tensor3d{0, 0, 0};
2139 std::cout <<
"time for integrals metric tform = " << timers[0].read(2)
2146 const auto nocc = Cocc.cols();
2147 Tensor2d Co{n, nocc};
2148 std::copy(Cocc.data(), Cocc.data() + n * nocc, Co.begin());
2149 Tensor3d xiK{n, nocc, ndf};
2150 btas::contract(1.0, xyK, {1, 2, 3}, Co, {2, 4}, 0.0, xiK, {1, 4, 3});
2153 btas::contract(1.0, xiK, {1, 2, 3}, xiK, {4, 2, 3}, 0.0, G, {1, 4});
2156 std::cout <<
"time for exchange = " << timers[0].read(3) << std::endl;
2162 btas::contract(1.0, xiK, {1, 2, 3}, Co, {1, 2}, 0.0, Jtmp, {3});
2163 xiK = Tensor3d{0, 0, 0};
2164 btas::contract(2.0, xyK, {1, 2, 3}, Jtmp, {3}, -1.0, G, {1, 2});
2167 std::cout <<
"time for coulomb = " << timers[0].read(4) << std::endl;
2170 Matrix result(n, n);
2171 std::copy(G.cbegin(), G.cend(), result.data());
2174 #endif // HAVE_DENSITY_FITTING 2177 void api_basic_compile_test(
const BasisSet& obs) {
2179 Engine onebody_engine(
2185 auto shell2bf = obs.shell2bf();
2186 const auto& results = onebody_engine.results();
2187 for (
auto s1 = 0; s1 != obs.size(); ++s1) {
2188 for (
auto s2 = 0; s2 != obs.size(); ++s2) {
2189 std::cout <<
"compute shell set {" << s1 <<
"," << s2 <<
"} ... ";
2190 onebody_engine.compute(obs[s1], obs[s2]);
2191 const auto* ints_shellset = results[0];
2192 std::cout <<
"done" << std::endl;
2194 auto bf1 = shell2bf[s1];
2195 auto n1 = obs[s1].size();
2196 auto bf2 = shell2bf[s2];
2197 auto n2 = obs[s2].size();
2201 for (
auto f1 = 0; f1 != n1; ++f1)
2202 for (
auto f2 = 0; f2 != n2; ++f2)
2203 std::cout <<
" " << bf1 + f1 <<
" " << bf2 + f2 <<
" " 2204 << ints_shellset[f1 * n2 + f2] << std::endl;
2208 using libint2::Operator;
2210 std::vector<std::pair<double, double>> cgtg_params{
2211 {0.1, 0.2}, {0.3, 0.4}, {0.5, 0.6}};
2214 compute_schwarz_ints<Operator::cgtg>(obs, obs,
false, cgtg_params);
2215 std::cout <<
"cGTG Schwarz ints\n" << K << std::endl;
2218 auto K = compute_schwarz_ints<Operator::cgtg_x_coulomb>(obs, obs,
false,
2220 std::cout <<
"cGTG/r12 Schwarz ints\n" << K << std::endl;
2224 compute_schwarz_ints<Operator::delcgtg2>(obs, obs,
false, cgtg_params);
2225 std::cout <<
"||Del.cGTG||^2 Schwarz ints\n" << K << std::endl;
2229 Engine eri4_engine(Operator::coulomb, obs.max_nprim(), obs.max_l());
2230 Engine eri2_engine = eri4_engine;
2231 eri2_engine.set(BraKet::xs_xs);
2232 auto shell2bf = obs.shell2bf();
2233 const auto& results4 = eri4_engine.results();
2234 const auto& results2 = eri2_engine.results();
2235 for (
auto s1 = 0; s1 != obs.size(); ++s1) {
2236 for (
auto s2 = 0; s2 != obs.size(); ++s2) {
2238 eri2_engine.compute(obs[s1], obs[s2]);
2240 auto bf1 = shell2bf[s1];
2241 auto n1 = obs[s1].size();
2242 auto bf2 = shell2bf[s2];
2243 auto n2 = obs[s2].size();
2245 const auto* buf4 = results4[0];
2246 const auto* buf2 = results2[0];
2250 for (
auto f1 = 0, f12 = 0; f1 != n1; ++f1)
2251 for (
auto f2 = 0; f2 != n2; ++f2, ++f12)
2252 assert(std::abs(buf4[f12] - buf2[f12]) < 1e-12 &&
2253 "2-center ints test failed");
2258 Engine eri4_engine(Operator::coulomb, obs.max_nprim(), obs.max_l());
2259 Engine eri3_engine = eri4_engine;
2260 eri3_engine.set(BraKet::xs_xx);
2261 auto shell2bf = obs.shell2bf();
2262 const auto& results4 = eri4_engine.results();
2263 const auto& results3 = eri3_engine.results();
2264 for (
auto s1 = 0; s1 != obs.size(); ++s1) {
2265 for (
auto s2 = 0; s2 != obs.size(); ++s2) {
2266 for (
auto s3 = 0; s3 != obs.size(); ++s3) {
2267 eri4_engine.compute(obs[s1],
Shell::unit(), obs[s2], obs[s3]);
2268 eri3_engine.compute(obs[s1], obs[s2], obs[s3]);
2270 auto bf1 = shell2bf[s1];
2271 auto n1 = obs[s1].size();
2272 auto bf2 = shell2bf[s2];
2275 auto bf3 = shell2bf[s3];
2276 auto n3 = obs[s3].size();
2278 const auto* buf4 = results4[0];
2279 const auto* buf3 = results3[0];
2283 for (
auto f1 = 0, f123 = 0; f1 != n1; ++f1)
2284 for (
auto f2 = 0; f2 != n2; ++f2)
2285 for (
auto f3 = 0; f3 != n3; ++f3, ++f123)
2286 assert(std::abs(buf4[f123] - buf3[f123]) < 1e-12 &&
2287 "3-center ints test failed");
2293 #if LIBINT2_DERIV_ERI_ORDER 2295 Engine eri4_engine(Operator::coulomb, obs.max_nprim(), obs.max_l(), 1);
2296 Engine eri2_engine = eri4_engine;
2297 eri2_engine.set(BraKet::xs_xs);
2298 auto shell2bf = obs.shell2bf();
2299 const auto& results4 = eri4_engine.results();
2300 const auto& results2 = eri2_engine.results();
2301 for (
auto s1 = 0; s1 != obs.size(); ++s1) {
2302 for (
auto s2 = 0; s2 != obs.size(); ++s2) {
2304 eri2_engine.compute(obs[s1], obs[s2]);
2306 auto bf1 = shell2bf[s1];
2307 auto n1 = obs[s1].size();
2308 auto bf2 = shell2bf[s2];
2309 auto n2 = obs[s2].size();
2312 for(
auto d=0; d!=6; ++d) {
2313 const auto* buf4 = results4[d<3 ? d : d+3];
2314 const auto* buf2 = results2[d];
2318 for (
auto f1 = 0, f12 = 0; f1 != n1; ++f1)
2319 for (
auto f2 = 0; f2 != n2; ++f2, ++f12)
2320 assert(std::abs(buf4[f12] - buf2[f12]) < 1e-12 &&
2321 "deriv 2-center ints test failed");
2329 #if LIBINT2_DERIV_ERI_ORDER > 1 2331 Engine eri4_engine(Operator::coulomb, obs.max_nprim(), obs.max_l(), 2);
2332 Engine eri2_engine = eri4_engine;
2333 eri2_engine.set(BraKet::xs_xs);
2334 auto shell2bf = obs.shell2bf();
2335 const auto& results4 = eri4_engine.results();
2336 const auto& results2 = eri2_engine.results();
2337 for (
auto s1 = 0; s1 != obs.size(); ++s1) {
2338 for (
auto s2 = 0; s2 != obs.size(); ++s2) {
2340 eri2_engine.compute(obs[s1], obs[s2]);
2342 auto bf1 = shell2bf[s1];
2343 auto n1 = obs[s1].size();
2344 auto bf2 = shell2bf[s2];
2345 auto n2 = obs[s2].size();
2348 for (
auto d1 = 0, d12 = 0; d1 != 6; ++d1) {
2349 const auto dd1 = d1 < 3 ? d1 : d1 + 3;
2350 for (
auto d2 = d1; d2 != 6; ++d2, ++d12) {
2351 const auto dd2 = d2 < 3 ? d2 : d2 + 3;
2352 const auto dd12 = dd1 * (24 - dd1 - 1) / 2 + dd2;
2353 const auto* buf4 = results4[dd12];
2354 const auto* buf2 = results2[d12];
2359 for (
auto f1 = 0, f12 = 0; f1 != n1; ++f1)
2360 for (
auto f2 = 0; f2 != n2; ++f2, ++f12)
2361 assert(std::abs(buf4[f12] - buf2[f12]) < 1e-12 &&
2362 "2nd deriv 2-center ints test failed");
std::vector< std::pair< double, std::array< double, 3 > > > make_point_charges(const std::vector< libint2::Atom > &atoms)
converts a vector of Atoms to a vector of point charges
Definition: atom.h:237
void clear()
resets timers to zero
Definition: timer.h:78
generally-contracted Solid-Harmonic/Cartesion Gaussian Shell
Definition: shell.h:81
static bool do_enforce_unit_normalization(defaultable_boolean flag=defaultable_boolean())
sets and/or reports whether the auto-renormalization to unity is set if called without arguments,...
Definition: shell.h:209
Defaults definitions for various parameters assumed by Libint.
Definition: algebra.cc:24
void parallel_do(Lambda &lambda)
fires off nthreads instances of lambda in parallel
Definition: 1body.h:168
DIIS (`‘direct inversion of iterative subspace’') extrapolation.
Definition: diis.h:132
dur_t stop(size_t t)
stops timer t
Definition: timer.h:67
void initialize(bool verbose=false)
initializes the libint library if not currently initialized, no-op otherwise
Definition: initialize.h:68
double read(size_t t) const
reads value (in seconds) of timer t , converted to double
Definition: timer.h:74
const std::vector< Real > & sto3g_ao_occupation_vector(size_t Z)
computes average orbital occupancies in the ground state of a neutral atoms
Definition: sto3g_atomic_density.h:84
std::vector< Atom > read_dotxyz(std::istream &is, const double bohr_to_angstrom=constants::codata_2010::bohr_to_angstrom)
reads the list of atoms from a file in the standard XYZ format supported by most chemistry software (...
Definition: atom.h:205
Array idential to C++0X arrays.
Definition: stdarray_bits.h:14
Iterates over all partitions of a non-negative integer into nonnegative integers in reverse lexicog...
Definition: intpart_iter.h:74
void start(size_t t)
starts timer t
Definition: timer.h:62
size_t sto3g_num_ao(size_t Z)
Definition: sto3g_atomic_density.h:56
void set_now_overhead(size_t ns)
use this to report the overhead of now() call; if set, the reported timings will be adjusted for this...
Definition: timer.h:57
static const Shell & unit()
Definition: shell.h:218
void finalize()
finalizes the libint library.
Definition: initialize.h:88
Timers aggregates N C++11 "timers"; used to high-resolution profile stages of integral computation.
Definition: timer.h:38