diff --git a/include/libint2/dfbs_generator.h b/include/libint2/dfbs_generator.h index 1a232307c..cdb13908c 100644 --- a/include/libint2/dfbs_generator.h +++ b/include/libint2/dfbs_generator.h @@ -149,6 +149,10 @@ namespace libint2 { return sorted_shells; } + /// @brief computes the reduced set of product functions via pivoted Cholesky decomposition + /// @param shells set of shells + /// @param cholesky_threshold threshold for choosing a product function via pivoted Cholesky decomposition + /// @return reduced set of product functions std::vector shell_pivoted_cholesky(const std::vector shells, const double cholesky_threshold) { auto n = shells.size(); // number of shells @@ -234,6 +238,16 @@ namespace libint2 { return candidate_shells_; } + /// @brief returns the candidate basis set (full set of product functions) + /// @warning generates huge and heavily linearly dependent basis sets + BasisSet product_basis(){ + std::vector product_shells; + for(auto shells: candidate_shells_){ + product_shells.insert(product_shells.end(), shells.begin(), shells.end()); + } + return BasisSet(std::move(product_shells)); + } + /// @brief returns the candidate shells sorted by angular momentum std::vector>> candidates_splitted_in_L() { std::vector>> sorted_shells; @@ -243,6 +257,7 @@ namespace libint2 { return sorted_shells; } + /// @brief returns the reduced shells (reduced set of product functions) computed via pivoted Cholesky decomposition std::vector> reduced_shells() { if (reduced_shells_.size() != 0) return reduced_shells_; @@ -261,6 +276,17 @@ namespace libint2 { return reduced_shells_; } + /// @brief returns the reduced basis set (reduced set of product functions) computed via pivoted Cholesky decomposition + BasisSet reduced_basis() { + auto reduced_cluster = reduced_shells(); + std::vector reduced_shells; + for (auto shells: reduced_cluster) { + reduced_shells.insert(reduced_shells.end(), shells.begin(), shells.end()); + } + return BasisSet(std::move(reduced_shells)); + } + + private: std::string obs_name_; //name of AO basis set std::vector atoms_; //vector of atoms diff --git a/tests/hartree-fock/hartree-fock++.cc b/tests/hartree-fock/hartree-fock++.cc index 44af2b2f9..edfd9f4d6 100644 --- a/tests/hartree-fock/hartree-fock++.cc +++ b/tests/hartree-fock/hartree-fock++.cc @@ -49,6 +49,8 @@ #include #include #include +#include + #if !LIBINT2_CONSTEXPR_STATICS # include #endif @@ -57,6 +59,7 @@ #include #endif + /// to use precomputed shell pair data must decide on max precision a priori const auto max_engine_precision = std::numeric_limits::epsilon() / 1e10; @@ -75,6 +78,8 @@ typedef Eigen::Matrix typedef Eigen::DiagonalMatrix DiagonalMatrix; +typedef btas::Tensor tensor; + using libint2::Shell; using libint2::Atom; using libint2::BasisSet; @@ -166,10 +171,8 @@ struct DFFockEngine { const BasisSet& dfbs; DFFockEngine(const BasisSet& _obs, const BasisSet& _dfbs) : obs(_obs), dfbs(_dfbs) {} - - typedef btas::RangeNd> Range3d; - typedef btas::Tensor Tensor3d; - Tensor3d xyK; + typedef btas::Tensor Tensor; + Tensor xyK; // a DF-based builder, using coefficients of occupied MOs Matrix compute_2body_fock_dfC(const Matrix& Cocc); @@ -216,21 +219,26 @@ int main(int argc, char* argv[]) { // filename (.xyz) from the command line const auto filename = (argc > 1) ? argv[1] : "h2o.xyz"; const auto basisname = (argc > 2) ? argv[2] : "aug-cc-pVDZ"; - bool do_density_fitting = false; + bool do_density_fitting = false; + double cholesky_threshold = 1e-4; #ifdef HAVE_DENSITY_FITTING - do_density_fitting = (argc > 3); - const auto dfbasisname = do_density_fitting ? argv[3] : ""; + do_density_fitting = (argc > 3); + const auto dfbasisname = do_density_fitting ? argv[3] : ""; + if ((strcmp(dfbasisname, "autodf") == 0) && argc > 4) { + cholesky_threshold = std::stod(argv[4]); + std::cout << "New Cholesky threshold for generating df basis set = " << cholesky_threshold << std::endl; + } #endif - std::vector atoms = read_geometry(filename); - - // set up thread pool - { - using libint2::nthreads; - auto nthreads_cstr = getenv("LIBINT_NUM_THREADS"); - nthreads = 1; - if (nthreads_cstr && strcmp(nthreads_cstr, "")) { - std::istringstream iss(nthreads_cstr); - iss >> nthreads; + std::vector atoms = read_geometry(filename); + + // set up thread pool + { + using libint2::nthreads; + auto nthreads_cstr = getenv("LIBINT_NUM_THREADS"); + nthreads = 1; + if (nthreads_cstr && strcmp(nthreads_cstr, "")) { + std::istringstream iss(nthreads_cstr); + iss >> nthreads; if (nthreads > 1 << 16 || nthreads <= 0) nthreads = 1; } #if defined(_OPENMP) @@ -254,33 +262,38 @@ int main(int argc, char* argv[]) { // compute the nuclear repulsion energy auto enuc = 0.0; for (auto i = 0; i < atoms.size(); i++) - for (auto j = i + 1; j < atoms.size(); j++) { - auto xij = atoms[i].x - atoms[j].x; - auto yij = atoms[i].y - atoms[j].y; - auto zij = atoms[i].z - atoms[j].z; - auto r2 = xij * xij + yij * yij + zij * zij; - auto r = sqrt(r2); - enuc += atoms[i].atomic_number * atoms[j].atomic_number / r; - } - cout << "Nuclear repulsion energy = " << std::setprecision(15) << enuc - << endl; + for (auto j = i + 1; j < atoms.size(); j++) { + auto xij = atoms[i].x - atoms[j].x; + auto yij = atoms[i].y - atoms[j].y; + auto zij = atoms[i].z - atoms[j].z; + auto r2 = xij * xij + yij * yij + zij * zij; + auto r = sqrt(r2); + enuc += atoms[i].atomic_number * atoms[j].atomic_number / r; + } + cout << "Nuclear repulsion energy = " << std::setprecision(15) << enuc + << endl; - libint2::Shell::do_enforce_unit_normalization(false); + libint2::Shell::do_enforce_unit_normalization(false); - cout << "Atomic Cartesian coordinates (a.u.):" << endl; - for (const auto& a : atoms) - std::cout << a.atomic_number << " " << a.x << " " << a.y << " " << a.z - << std::endl; + cout << "Atomic Cartesian coordinates (a.u.):" << endl; + for (const auto &a: atoms) + std::cout << a.atomic_number << " " << a.x << " " << a.y << " " << a.z + << std::endl; - BasisSet obs(basisname, atoms); - cout << "orbital basis set rank = " << obs.nbf() << endl; + BasisSet obs(basisname, atoms); + cout << "orbital basis set rank = " << obs.nbf() << endl; #ifdef HAVE_DENSITY_FITTING - BasisSet dfbs; - if (do_density_fitting) { - dfbs = BasisSet(dfbasisname, atoms); - cout << "density-fitting basis set rank = " << dfbs.nbf() << endl; - } + BasisSet dfbs; + if (do_density_fitting) { + if (strcmp(dfbasisname, "autodf") == 0) { + if (!libint2::initialized()) libint2::initialize(); + libint2::DFBasisSetGenerator dfbs_generator(basisname, atoms, cholesky_threshold); + dfbs = dfbs_generator.reduced_basis(); + } else + dfbs = BasisSet(dfbasisname, atoms); + cout << "density-fitting basis set rank = " << dfbs.nbf() << endl; + } #endif // HAVE_DENSITY_FITTING /*** =========================== ***/ @@ -288,34 +301,36 @@ int main(int argc, char* argv[]) { /*** =========================== ***/ // initializes the Libint integrals library ... now ready to compute - libint2::initialize(); - - // compute OBS non-negligible shell-pair list - { - const auto tstart = std::chrono::high_resolution_clock::now(); - std::tie(obs_shellpair_list, obs_shellpair_data) = compute_shellpairs(obs); - size_t nsp = 0; - for (auto& sp : obs_shellpair_list) { - nsp += sp.second.size(); + if (!libint2::initialized()) + libint2::initialize(); + + // compute OBS non-negligible shell-pair list + { + const auto tstart = std::chrono::high_resolution_clock::now(); + std::tie(obs_shellpair_list, obs_shellpair_data) = compute_shellpairs(obs); + size_t nsp = 0; + for (auto &sp: obs_shellpair_list) { + nsp += sp.second.size(); + } + const auto tstop = std::chrono::high_resolution_clock::now(); + const std::chrono::duration time_elapsed = tstop - tstart; + std::cout << "computed shell-pair data in " << time_elapsed.count() + << " seconds: # of {all,non-negligible} shell-pairs = {" + << obs.size() * (obs.size() + 1) / 2 << "," << nsp << "}" + << std::endl; } - const auto tstop = std::chrono::high_resolution_clock::now(); - const std::chrono::duration time_elapsed = tstop - tstart; - std::cout << "computed shell-pair data in " << time_elapsed.count() << " seconds: # of {all,non-negligible} shell-pairs = {" - << obs.size() * (obs.size() + 1) / 2 << "," << nsp << "}" - << std::endl; - } - // compute one-body integrals - auto S = compute_1body_ints(obs)[0]; - auto T = compute_1body_ints(obs)[0]; - auto V = compute_1body_ints(obs, libint2::make_point_charges(atoms))[0]; - Matrix H = T + V; - T.resize(0, 0); - V.resize(0, 0); - - // compute orthogonalizer X such that X.transpose() . S . X = I - Matrix X, Xinv; - double XtX_condition_number; // condition number of "re-conditioned" + // compute one-body integrals + auto S = compute_1body_ints(obs)[0]; + auto T = compute_1body_ints(obs)[0]; + auto V = compute_1body_ints(obs, libint2::make_point_charges(atoms))[0]; + Matrix H = T + V; + T.resize(0, 0); + V.resize(0, 0); + + // compute orthogonalizer X such that X.transpose() . S . X = I + Matrix X, Xinv; + double XtX_condition_number; // condition number of "re-conditioned" // overlap obtained as Xinv.transpose() . Xinv // one should think of columns of Xinv as the conditioned basis // Re: name ... cond # (Xinv.transpose() . Xinv) = cond # (X.transpose() . @@ -2162,11 +2177,6 @@ Matrix DFFockEngine::compute_2body_fock_dfC(const Matrix& Cocc) { std::vector> timers(nthreads); for(auto& timer: timers) timer.set_now_overhead(25); - typedef btas::RangeNd> Range1d; - typedef btas::RangeNd> Range2d; - typedef btas::Tensor Tensor1d; - typedef btas::Tensor Tensor2d; - // using first time? compute 3-center ints and transform to inv sqrt // representation if (xyK.size() == 0) { @@ -2192,7 +2202,7 @@ Matrix DFFockEngine::compute_2body_fock_dfC(const Matrix& Cocc) { auto shell2bf = obs.shell2bf(); auto shell2bf_df = dfbs.shell2bf(); - Tensor3d Zxy{ndf, n, n}; + Tensor Zxy{ndf, n, n}; auto lambda = [&](int thread_id) { @@ -2269,12 +2279,12 @@ Matrix DFFockEngine::compute_2body_fock_dfC(const Matrix& Cocc) { // std::cout << "||V^-1 - L^-1^t L^-1|| = " << (V.inverse() - Linv_t * // Linv_t.transpose()).norm() << std::endl; - Tensor2d K{ndf, ndf}; + Tensor K{ndf, ndf}; std::copy(Linv_t.data(), Linv_t.data() + ndf * ndf, K.begin()); - xyK = Tensor3d{n, n, ndf}; + xyK = Tensor{n, n, ndf}; btas::contract(1.0, Zxy, {1, 2, 3}, K, {1, 4}, 0.0, xyK, {2, 3, 4}); - Zxy = Tensor3d{0, 0, 0}; // release memory + Zxy = Tensor{0, 0, 0}; // release memory timers[0].stop(2); std::cout << "time for integrals metric tform = " << timers[0].read(2) @@ -2285,12 +2295,12 @@ Matrix DFFockEngine::compute_2body_fock_dfC(const Matrix& Cocc) { timers[0].start(3); const auto nocc = Cocc.cols(); - Tensor2d Co{n, nocc}; + Tensor Co{n, nocc}; std::copy(Cocc.data(), Cocc.data() + n * nocc, Co.begin()); - Tensor3d xiK{n, nocc, ndf}; + Tensor xiK{n, nocc, ndf}; btas::contract(1.0, xyK, {1, 2, 3}, Co, {2, 4}, 0.0, xiK, {1, 4, 3}); - Tensor2d G{n, n}; + Tensor G{n, n}; btas::contract(1.0, xiK, {1, 2, 3}, xiK, {4, 2, 3}, 0.0, G, {1, 4}); timers[0].stop(3); @@ -2299,9 +2309,9 @@ Matrix DFFockEngine::compute_2body_fock_dfC(const Matrix& Cocc) { // compute Coulomb timers[0].start(4); - Tensor1d Jtmp{ndf}; + Tensor Jtmp{ndf}; btas::contract(1.0, xiK, {1, 2, 3}, Co, {1, 2}, 0.0, Jtmp, {3}); - xiK = Tensor3d{0, 0, 0}; + xiK = Tensor{0, 0, 0}; btas::contract(2.0, xyK, {1, 2, 3}, Jtmp, {3}, -1.0, G, {1, 2}); timers[0].stop(4);