diff --git a/forte/CMakeLists.txt b/forte/CMakeLists.txt index 3c20c0f69..d08f33f5d 100644 --- a/forte/CMakeLists.txt +++ b/forte/CMakeLists.txt @@ -212,12 +212,15 @@ mrdsrg-spin-integrated/mrdsrg_pt.cc mrdsrg-spin-integrated/mrdsrg_smart_s.cc mrdsrg-spin-integrated/mrdsrg_srg.cc mrdsrg-spin-integrated/three_dsrg_mrpt2.cc +mrdsrg-spin-integrated/Laplace_dsrg.cc orbital-helpers/ao_helper.cc orbital-helpers/aosubspace.cc orbital-helpers/ci-no/ci-no.cc orbital-helpers/ci-no/mrci-no.cc orbital-helpers/fragment_projector.cc orbital-helpers/iao_builder.cc +orbital-helpers/Laplace.cc +orbital-helpers/LaplaceDenominator.cc orbital-helpers/localize.cc orbital-helpers/mp2_nos.cc orbital-helpers/mrpt2_nos.cc diff --git a/forte/mrdsrg-spin-integrated/Laplace_dsrg.cc b/forte/mrdsrg-spin-integrated/Laplace_dsrg.cc new file mode 100644 index 000000000..c17c129db --- /dev/null +++ b/forte/mrdsrg-spin-integrated/Laplace_dsrg.cc @@ -0,0 +1,1032 @@ +#include +#include +#include +#include +#include +#include + +#include "psi4/libpsi4util/PsiOutStream.h" +#include "psi4/libpsi4util/process.h" +#include "psi4/lib3index/dftensor.h" +#include "psi4/libmints/wavefunction.h" +#include "psi4/libmints/matrix.h" +#include "psi4/libmints/molecule.h" +#include "psi4/libmints/vector.h" +#include "psi4/libpsio/psio.h" +#include "psi4/libpsio/psio.hpp" +#include "psi4/libqt/qt.h" +#include "psi4/libmints/integral.h" +#include "psi4/psifiles.h" +#include "psi4/libfock/jk.h" + +#include "orbital-helpers/ao_helper.h" +#include "helpers/blockedtensorfactory.h" +#include "helpers/printing.h" +#include "helpers/timer.h" +#include "fci/fci_solver.h" +#include "fci/fci_vector.h" +#include "sci/aci.h" + +#include "orbital-helpers/Laplace.h" +#include "Laplace_dsrg.h" + +using namespace ambit; + +using namespace psi; + +namespace forte { +LaplaceDSRG::LaplaceDSRG(std::shared_ptr options, + std::shared_ptr ints, + std::shared_ptr mo_space_info, + psi::SharedVector epsilon_rdocc, psi::SharedVector epsilon_virtual, + psi::SharedVector epsilon_active, psi::SharedMatrix Gamma1, + psi::SharedMatrix Eta1) + : foptions_(options), ints_(ints), mo_space_info_(mo_space_info), eps_rdocc_(epsilon_rdocc), + eps_virtual_(epsilon_virtual), eps_active_(epsilon_active), Gamma1_mat_(Gamma1), + Eta1_mat_(Eta1) { + + /// Number of MO. + nmo_ = mo_space_info_->dimension("ALL").sum(); + nthree_ = ints_->nthree(); + Cwfn_ = ints_->Ca(); + nfrozen_ = mo_space_info_->dimension("FROZEN").sum(); + ncore_ = eps_rdocc_->dim(); + nactive_ = eps_active_->dim(); + nvirtual_ = eps_virtual_->dim(); + vir_tol_ = foptions_->get_double("VIR_TOL"); + /// CCVV + theta_NB_ = foptions_->get_double("THETA_NB"); + theta_NB_IAP_ = foptions_->get_double("THETA_NB_IAP"); + theta_ij_ = foptions_->get_double("THETA_IJ"); + Omega_ = foptions_->get_double("OMEGA"); + theta_schwarz_ = foptions_->get_double("THETA_SCHWARZ"); + laplace_threshold_ = foptions_->get_double("LAPLACE_THRESHOLD"); + theta_ij_sqrt_ = sqrt(theta_ij_); + /// CAVV + theta_NB_cavv_ = foptions_->get_double("THETA_NB_CAVV"); + theta_NB_IAP_cavv_ = foptions_->get_double("THETA_NB_IAP_CAVV"); + theta_ij_cavv_ = foptions_->get_double("THETA_IJ_CAVV"); + theta_schwarz_cavv_ = foptions_->get_double("THETA_SCHWARZ_CAVV"); + theta_ij_sqrt_cavv_ = sqrt(theta_ij_cavv_); + theta_XNB_cavv_ = foptions_->get_double("THETA_XNB_CAVV"); + theta_NB_XAP_cavv_ = foptions_->get_double("THETA_NB_XAP_CAVV"); + /// CCAV + theta_NB_ccav_ = foptions_->get_double("THETA_NB_CCAV"); + theta_NB_IAP_ccav_ = foptions_->get_double("THETA_NB_IAP_CCAV"); + theta_ij_ccav_ = foptions_->get_double("THETA_IJ_CCAV"); + theta_schwarz_ccav_ = foptions_->get_double("THETA_SCHWARZ_CCAV"); + theta_ij_sqrt_ccav_ = sqrt(theta_ij_ccav_); + + C_pq_ = erfc_metric(Omega_, ints_); + + outfile->Printf("\n Done with C_pq"); + S_ = ints_->wfn()->S(); + outfile->Printf("\n Done with S"); + + primary_ = ints_->wfn()->basisset(); + auxiliary_ = ints_->wfn()->get_basisset("DF_BASIS_MP2"); + + DiskDFJK disk_jk(primary_, auxiliary_); + disk_jk.erfc_three_disk(Omega_); + + // std::cout << ints_.use_count() << std::endl; + + // ints_->~ForteIntegrals(); + + // std::cout << ints_.use_count() << std::endl; + + print_header(); +} + +LaplaceDSRG::~LaplaceDSRG() { outfile->Printf("\n\n Done with Laplace_DSRG class"); } + +void LaplaceDSRG::print_header() { + print_method_banner({"LAPLACE-DSRG-MRPT2", "Shuhang Li"}); + std::vector> calculation_info_int; + std::vector> calculation_info_string; + std::vector> calculation_info_bool; + std::vector> calculation_info_double{ + {"laplace_threshold", laplace_threshold_}, + {"Omega", Omega_}, + {"theta_NB", theta_NB_}, + {"theta_NB_IAP", theta_NB_IAP_}, + {"theta_ij", theta_ij_}, + {"theta_schwarz", theta_schwarz_}, + {"theta_NB_cavv", theta_NB_cavv_}, + {"theta_NB_IAP_cavv", theta_NB_IAP_cavv_}, + {"theta_XNB_cavv", theta_XNB_cavv_}, + {"theta_NB_XAP_cavv", theta_NB_XAP_cavv_}, + {"theta_ij_cavv", theta_ij_cavv_}, + {"theta_schwarz_cavv", theta_schwarz_cavv_}, + {"theta_NB_ccav", theta_NB_ccav_}, + {"theta_NB_IAP_ccav", theta_NB_IAP_ccav_}, + {"theta_ij_ccav", theta_ij_ccav_}, + {"theta_schwarz_ccav", theta_schwarz_ccav_}, + {"vir_tol", vir_tol_}}; + print_selected_options("Calculation Information", calculation_info_string, + calculation_info_bool, calculation_info_double, calculation_info_int); +} + +// void LaplaceDSRG::prepare_cholesky_coeff(&psi::SharedMatrix C_wfn, &psi::SharedVector eps_rdocc, +// &psi::SharedVector eps_virtual, +// &double laplace_threshold, &size_t nactive, &size_t +// nfrozen, &double vir_tol) { + +// } + +void LaplaceDSRG::prepare_cholesky_coeff(std::shared_ptr ao_helper_ptr, + const std::string& algorithm) { + weights_ = ao_helper_ptr->Weights(); + vir_start_ = ao_helper_ptr->vir_start(); + ao_helper_ptr->Compute_Cholesky_Pseudo_Density(); + if (algorithm == "ccvv") { + ao_helper_ptr->Compute_Cholesky_Density(); + Cholesky_Occ_ = ao_helper_ptr->L_Occ_real(); + Cholesky_Occ_abs_ = + std::make_shared("LOcc_abs", nmo_, Cholesky_Occ_->coldim()); + for (size_t u = 0; u < nmo_; u++) { + for (size_t i = 0; i < Cholesky_Occ_->coldim(); i++) { + Cholesky_Occ_abs_->set(u, i, std::abs(Cholesky_Occ_->get(u, i))); + } + } + + } else if (algorithm == "cavv") { + // I assume that Cholesky_Occ_ has been initialized in the CCVV algorithm. + ao_helper_ptr->Compute_Cholesky_Active_Density(Gamma1_mat_); + Active_cholesky_ = ao_helper_ptr->LAct_list(); + Active_cholesky_abs_.resize(weights_); + } else if (algorithm == "ccav") { + ao_helper_ptr->Compute_Cholesky_Active_Density(Eta1_mat_); + Active_cholesky_ = ao_helper_ptr->LAct_list(); + Active_cholesky_abs_.resize(weights_); + } + + Occupied_cholesky_ = ao_helper_ptr->LOcc_list(); + Virtual_cholesky_ = ao_helper_ptr->LVir_list(); + Occupied_cholesky_abs_.resize(weights_); + Virtual_cholesky_abs_.resize(weights_); + + for (size_t nweight = 0; nweight < weights_; nweight++) { + Occupied_cholesky_abs_[nweight] = + std::make_shared("LOcc_abs", nmo_, Occupied_cholesky_[nweight]->coldim()); + Virtual_cholesky_abs_[nweight] = + std::make_shared("LVir_abs", nmo_, Virtual_cholesky_[nweight]->coldim()); + + if (algorithm != "ccvv") { + Active_cholesky_abs_[nweight] = std::make_shared( + "LAct_abs", nmo_, Active_cholesky_[nweight]->coldim()); + } + + for (size_t u = 0; u < nmo_; u++) { + for (size_t i = 0; i < Occupied_cholesky_[nweight]->coldim(); i++) { + Occupied_cholesky_abs_[nweight]->set( + u, i, std::abs(Occupied_cholesky_[nweight]->get(u, i))); + } + for (size_t a = 0; a < Virtual_cholesky_[nweight]->coldim(); a++) { + Virtual_cholesky_abs_[nweight]->set( + u, a, std::abs(Virtual_cholesky_[nweight]->get(u, a))); + } + if (algorithm != "ccvv") { + for (size_t v = 0; v < Active_cholesky_[nweight]->coldim(); v++) { + Active_cholesky_abs_[nweight]->set( + u, v, std::abs(Active_cholesky_[nweight]->get(u, v))); + } + } + } + } + + T_ibar_i_list_.clear(); + T_ibar_i_list_.resize(weights_); + for (size_t nweight = 0; nweight < weights_; nweight++) { + T_ibar_i_list_[nweight] = psi::linalg::triplet(Occupied_cholesky_[nweight], S_, + Cholesky_Occ_, true, false, false); + } +} + +void LaplaceDSRG::Load_int(double& theta_NB, const std::string& algorithm) { + local_timer loadingTime; + int file_unit = PSIF_DFSCF_BJ; + size_t n_func_pairs = (nmo_ + 1) * nmo_ / 2; + size_t max_rows = 1; + std::shared_ptr psio(new PSIO()); + psi::SharedMatrix loadAOtensor = + std::make_shared("DiskDF: Load (A|mn) from DF-SCF", max_rows, n_func_pairs); + psio->open(file_unit, PSIO_OPEN_OLD); + + std::vector amn; + std::vector amn_new; + + ao_list_per_q_.clear(); + ao_list_per_q_.resize(nthree_); + + N_pu_ = std::make_shared("N_pu_", nthree_, nmo_); + N_pu_->zero(); + double* N_pu_p = N_pu_->get_pointer(); + outfile->Printf("\n Start loading"); + + P_xbar_u_cavv_.assign(weights_, std::vector(nthree_)); + P_xbar_u_ccav_.assign(weights_, std::vector(nthree_)); + + for (int Q = 0; Q < nthree_; Q += max_rows) { + int naux = (nthree_ - Q <= max_rows ? nthree_ - Q : max_rows); + psio_address addr = + psio_get_address(PSIO_ZERO, (Q * (size_t)n_func_pairs) * sizeof(double)); + psio->read(file_unit, "ERFC Integrals", (char*)(loadAOtensor->pointer()[0]), + sizeof(double) * naux * n_func_pairs, addr, &addr); + + double* add = loadAOtensor->get_pointer(); + + amn.resize(naux); + amn_new.resize(naux); + + psi::SharedMatrix N_pu_batch = std::make_shared("N_pu_batch", naux, nmo_); + N_pu_batch->zero(); + double* N_pu_batch_p = N_pu_batch->get_pointer(); + + for (int q = 0; q < naux; q++) { + amn[q] = std::make_shared("(m * n)", nmo_, nmo_); + /// Fill amn + for (int u_idx = 0; u_idx < nmo_; u_idx++) { + for (int v_idx = 0; v_idx <= u_idx; v_idx++) { + amn[q]->set(u_idx, v_idx, *add); + amn[q]->set(v_idx, u_idx, *add); + add++; + } + } + /// Loop over amn + double* row_add = amn[q]->get_pointer(); + for (int u = 0; u < nmo_; u++) { + double max_per_uq = (double)0.0; + for (int v = 0; v < nmo_; v++) { + double row_add_abs = std::abs(*row_add); + max_per_uq = + std::abs(max_per_uq) > row_add_abs ? std::abs(max_per_uq) : row_add_abs; + //*M_uv_p = std::abs(*M_uv_p) > row_add_abs ? std::abs(*M_uv_p) : row_add_abs; + row_add++; + // M_uv_p++; + } + if (std::abs(max_per_uq) >= theta_NB) { + ao_list_per_q_[Q + q].emplace_back(u); + } + *N_pu_p = std::abs(*N_pu_p) > max_per_uq ? std::abs(*N_pu_p) : max_per_uq; + *N_pu_batch_p = + std::abs(*N_pu_batch_p) > max_per_uq ? std::abs(*N_pu_batch_p) : max_per_uq; + N_pu_p++; + N_pu_batch_p++; + } + /// Slice Amn -> Amn_new + amn_new[q] = + submatrix_rows_and_cols(*amn[q], ao_list_per_q_[Q + q], ao_list_per_q_[Q + q]); + } + if (algorithm == "ccvv") { + P_iu_.resize(nthree_); + SparseMap i_p_up(nthree_); + i_p_.resize(nthree_); + SparseMap i_p_for_i_up(nthree_); + psi::SharedMatrix N_pi_batch = + psi::linalg::doublet(N_pu_batch, Cholesky_Occ_abs_, false, false); + double* N_pi_batch_p = N_pi_batch->get_pointer(); + for (int q = 0; q < naux; q++) { + for (int i = 0; i < Cholesky_Occ_->coldim(); i++) { + if (std::abs(*N_pi_batch_p) >= theta_NB_) { + i_p_up[Q + q].emplace_back(i); + } + N_pi_batch_p++; + } + psi::SharedMatrix Cholesky_Occ_new = + submatrix_rows_and_cols(*Cholesky_Occ_, ao_list_per_q_[Q + q], i_p_up[Q + q]); + P_iu_[Q + q] = psi::linalg::doublet(Cholesky_Occ_new, amn_new[q], true, false); + for (int inew = 0; inew < i_p_up[Q + q].size(); inew++) { + for (int u = 0; u < ao_list_per_q_[Q + q].size(); u++) { + if (std::abs(P_iu_[Q + q]->get(inew, u)) >= theta_NB_) { + i_p_for_i_up[Q + q].emplace_back(inew); + i_p_[Q + q].emplace_back(i_p_up[Q + q][inew]); + break; + } + } + } + P_iu_[Q + q] = submatrix_rows(*P_iu_[Q + q], i_p_for_i_up[Q + q]); + } + } else if (algorithm == "cavv") { + std::vector xbar_p_up_cavv(weights_, SparseMap(nthree_)); + std::vector xbar_u_p_for_xbar_up(weights_, SparseMap(nthree_)); + xbar_u_p_.assign(weights_, SparseMap(nthree_)); + + for (int nweight = 0; nweight < weights_; nweight++) { + psi::SharedMatrix N_px_cavv_batch = + psi::linalg::doublet(N_pu_batch, Active_cholesky_abs_[nweight]); + double* N_px_batch_p = N_px_cavv_batch->get_pointer(); + for (int q = 0; q < naux; q++) { + for (int x = 0; x < Active_cholesky_abs_[nweight]->coldim(); x++) { + if (std::abs(*N_px_batch_p) >= theta_XNB_cavv_) { + xbar_p_up_cavv[nweight][Q + q].emplace_back(x); + } + N_px_batch_p++; + } + psi::SharedMatrix Active_cholesky_new = + submatrix_rows_and_cols(*Active_cholesky_[nweight], ao_list_per_q_[Q + q], + xbar_p_up_cavv[nweight][Q + q]); + P_xbar_u_cavv_[nweight][Q + q] = + psi::linalg::doublet(Active_cholesky_new, amn_new[q], true, false); + for (int inew = 0; inew < xbar_p_up_cavv[nweight][Q + q].size(); inew++) { + for (int u = 0; u < ao_list_per_q_[Q + q].size(); u++) { + if (std::abs(P_xbar_u_cavv_[nweight][Q + q]->get(inew, u)) >= + theta_XNB_cavv_) { + xbar_u_p_for_xbar_up[nweight][Q + q].emplace_back(inew); + xbar_u_p_[nweight][Q + q].emplace_back( + xbar_p_up_cavv[nweight][Q + q][inew]); + break; + } + } + } + P_xbar_u_cavv_[nweight][Q + q] = submatrix_rows( + *P_xbar_u_cavv_[nweight][Q + q], xbar_u_p_for_xbar_up[nweight][Q + q]); + } + } + } else if (algorithm == "ccav") { + std::vector xbar_p_up_ccav(weights_, SparseMap(nthree_)); + std::vector xbar_u_p_for_xbar_up(weights_, SparseMap(nthree_)); + std::vector xbar_u_p_(weights_, SparseMap(nthree_)); + + for (int nweight = 0; nweight < weights_; nweight++) { + psi::SharedMatrix N_px_ccav_batch = + psi::linalg::doublet(N_pu_batch, Active_cholesky_abs_[nweight]); + double* N_px_batch_p = N_px_ccav_batch->get_pointer(); + for (int q = 0; q < naux; q++) { + for (int x = 0; x < Active_cholesky_abs_[nweight]->coldim(); x++) { + if (std::abs(*N_px_batch_p) >= theta_NB) { + xbar_p_up_ccav[nweight][Q + q].emplace_back(x); + } + N_px_batch_p++; + } + psi::SharedMatrix Active_cholesky_new = + submatrix_rows_and_cols(*Active_cholesky_[nweight], ao_list_per_q_[Q + q], + xbar_p_up_ccav[nweight][Q + q]); + P_xbar_u_ccav_[nweight][Q + q] = + psi::linalg::doublet(Active_cholesky_new, amn_new[q], true, false); + for (int inew = 0; inew < xbar_p_up_ccav[nweight][Q + q].size(); inew++) { + for (int u = 0; u < ao_list_per_q_[Q + q].size(); u++) { + if (std::abs(P_xbar_u_ccav_[nweight][Q + q]->get(inew, u)) >= + theta_NB) { + xbar_u_p_for_xbar_up[nweight][Q + q].emplace_back(inew); + xbar_u_p_[nweight][Q + q].emplace_back( + xbar_p_up_ccav[nweight][Q + q][inew]); + break; + } + } + } + P_xbar_u_ccav_[nweight][Q + q] = submatrix_rows( + *P_xbar_u_ccav_[nweight][Q + q], xbar_u_p_for_xbar_up[nweight][Q + q]); + } + } + } + } + psio->close(file_unit, 1); + outfile->Printf("\n End loading"); + outfile->Printf("\n\n Loading takes %8.8f", loadingTime.get()); +} + +void LaplaceDSRG::clear_maps(const std::string& algorithm) { + i_bar_p_up_.clear(); + a_bar_p_up_.clear(); + ibar_p_.clear(); + xbar_p_.clear(); + ibar_x_p_.clear(); + P_ibar_.clear(); + P_xbar_.clear(); + P_ibar_x_.clear(); + abar_ibar_.clear(); + abar_xbar_.clear(); + xbar_ibar_.clear(); + P_ibar_u_.clear(); + P_ibar_abar_.clear(); + P_xbar_abar_.clear(); + P_ibar_xbar_.clear(); + i_bar_a_bar_P_.clear(); + x_bar_a_bar_P_.clear(); + i_bar_x_bar_P_.clear(); + i_bar_a_bar_P_sliced_.clear(); + x_bar_a_bar_P_sliced_.clear(); + i_bar_x_bar_P_sliced_.clear(); + B_ia_Q_.clear(); + B_xa_Q_.clear(); + B_ix_Q_.clear(); + vir_intersection_per_ij_.clear(); + aux_intersection_per_ij_.clear(); + aux_in_B_i_.clear(); + aux_in_B_j_.clear(); + aux_in_B_ix_.clear(); + vir_intersection_per_ix_.clear(); + aux_intersection_per_ix_.clear(); + aux_intersection_per_ijx_.clear(); + aux_intersection_per_ixj_.clear(); + + i_bar_p_up_.resize(nthree_); + a_bar_p_up_.resize(nthree_); + P_ibar_u_.resize(nthree_); + P_ibar_abar_.resize(nthree_); + ibar_p_.resize(nthree_); + Z_pq_ = std::make_shared("Z_pq", nthree_, nthree_); + + if (algorithm == "cavv") { + P_xbar_abar_.resize(nthree_); + xbar_p_.resize(nthree_); + ZA_pq_ = std::make_shared("ZA_pq", nthree_, nthree_); + } else if (algorithm == "ccav") { + P_ibar_xbar_.resize(nthree_); + ibar_x_p_.resize(nthree_); + ZA_pq_ = std::make_shared("ZA_pq", nthree_, nthree_); + } +} + +void LaplaceDSRG::fill_maps(const std::string& algorithm, const int& nweight, const int& nocc, + const int& nact, const int& nvir) { + double theta_NB; + double theta_NB_IAP; + + if (algorithm == "ccvv") { + theta_NB = theta_NB_; + theta_NB_IAP = theta_NB_IAP_; + } else if (algorithm == "cavv") { + theta_NB = theta_NB_cavv_; + theta_NB_IAP = theta_NB_IAP_cavv_; + } else if (algorithm == "ccav") { + theta_NB = theta_NB_ccav_; + theta_NB_IAP = theta_NB_IAP_ccav_; + } + + psi::SharedMatrix N_pi_bar = + psi::linalg::doublet(N_pu_, Occupied_cholesky_abs_[nweight], false, false); + double* N_pi_bar_p = N_pi_bar->get_pointer(); + psi::SharedMatrix N_pa_bar = + psi::linalg::doublet(N_pu_, Virtual_cholesky_abs_[nweight], false, false); + double* N_pa_bar_p = N_pa_bar->get_pointer(); + + for (int qa = 0; qa < nthree_; qa++) { + i_bar_p_up_[qa].clear(); + a_bar_p_up_[qa].clear(); + ibar_p_[qa].clear(); + for (int ibar = 0; ibar < nocc; ibar++) { + if (std::abs(*N_pi_bar_p) >= theta_NB) { + i_bar_p_up_[qa].emplace_back(ibar); + } + N_pi_bar_p++; + } + + psi::SharedMatrix T_ibar_i_new = + submatrix_rows_and_cols(*T_ibar_i_list_[nweight], i_bar_p_up_[qa], i_p_[qa]); + P_ibar_u_[qa] = psi::linalg::doublet(T_ibar_i_new, P_iu_[qa], false, false); + + for (int abar = 0; abar < nvir; abar++) { + if (std::abs(*N_pa_bar_p) >= theta_NB) { + a_bar_p_up_[qa].emplace_back(abar); + } + N_pa_bar_p++; + } + + psi::SharedMatrix Pseudo_Vir_Mo = submatrix_rows_and_cols( + *Virtual_cholesky_[nweight], ao_list_per_q_[qa], a_bar_p_up_[qa]); + P_ibar_abar_[qa] = psi::linalg::doublet(P_ibar_u_[qa], Pseudo_Vir_Mo, false, false); + + /// Construct {ibar}_p. + for (int i = 0; i < P_ibar_abar_[qa]->rowdim(); i++) { + for (int a = 0; a < P_ibar_abar_[qa]->coldim(); a++) { + if (std::abs(P_ibar_abar_[qa]->get(i, a)) >= theta_NB_IAP) { + ibar_p_[qa].emplace_back(i_bar_p_up_[qa][i]); + break; + } + } + } + if (algorithm == "cavv") { + xbar_p_[qa].clear(); + P_xbar_abar_[qa] = + psi::linalg::doublet(P_xbar_u_cavv_[nweight][qa], Pseudo_Vir_Mo, false, false); + /// Construct {xbar}_p. + for (int x = 0; x < P_xbar_abar_[qa]->rowdim(); x++) { + for (int a = 0; a < P_xbar_abar_[qa]->coldim(); a++) { + if (std::abs(P_xbar_abar_[qa]->get(x, a)) >= theta_NB_XAP_cavv_) { + xbar_p_[qa].emplace_back(xbar_u_p_[nweight][qa][x]); + break; + } + } + } + } else if (algorithm == "ccav") { + ibar_x_p_[qa].clear(); + psi::SharedMatrix Pseudo_Occ_Mo = submatrix_rows_and_cols( + *Occupied_cholesky_[nweight], ao_list_per_q_[qa], i_bar_p_up_[qa]); + P_ibar_xbar_[qa] = + psi::linalg::doublet(Pseudo_Occ_Mo, P_xbar_u_ccav_[nweight][qa], true, true); + /// Construct {ibar_x}_p. + for (int i = 0; i < P_ibar_xbar_[qa]->rowdim(); i++) { + for (int x = 0; x < P_ibar_xbar_[qa]->coldim(); x++) { + if (std::abs(P_ibar_xbar_[qa]->get(i, x)) >= theta_NB_IAP_ccav_) { + ibar_x_p_[qa].emplace_back(i_bar_p_up_[qa][i]); + break; + } + } + } + } + } + /// Reorder. Construct (ibar abar|P). Use all ibar and all abar. + i_bar_a_bar_P_.resize(nocc); + i_bar_a_bar_P_sliced_.resize(nocc); + if (algorithm == "ccav") { + i_bar_x_bar_P_.resize(nocc); + i_bar_x_bar_P_sliced_.resize(nocc); + xbar_ibar_.resize(nocc); + } + + for (int i = 0; i < nocc; i++) { + i_bar_a_bar_P_[i] = std::make_shared("(abar * P)", nvir, nthree_); + i_bar_a_bar_P_[i]->zero(); + if (algorithm == "ccav") { + i_bar_x_bar_P_[i] = std::make_shared("(xbar * P)", nact, nthree_); + i_bar_x_bar_P_[i]->zero(); + } + } + + for (int q = 0; q < nthree_; q++) { + for (int i = 0; i < P_ibar_abar_[q]->rowdim(); i++) { + for (int a = 0; a < P_ibar_abar_[q]->coldim(); a++) { + i_bar_a_bar_P_[i_bar_p_up_[q][i]]->set(a_bar_p_up_[q][a], q, + P_ibar_abar_[q]->get(i, a)); + } + if (algorithm == "ccav") { + for (int x = 0; x < P_ibar_xbar_[q]->coldim(); x++) { + i_bar_x_bar_P_[i_bar_p_up_[q][i]]->set(xbar_u_p_[nweight][q][x], q, + P_ibar_xbar_[q]->get(i, x)); + } + } + } + } + + /// Construct {abar}_ibar. + abar_ibar_.resize(nocc); + for (int i = 0; i < nocc; i++) { + abar_ibar_[i].clear(); + for (int a = 0; a < nvir; a++) { + for (int q = 0; q < nthree_; q++) { + if (std::abs(i_bar_a_bar_P_[i]->get(a, q)) >= theta_NB_IAP) { + abar_ibar_[i].emplace_back(a); + break; + } + } + } + /// Construct {xbar}_ibar. + if (algorithm == "ccav") { + xbar_ibar_[i].clear(); + for (int x = 0; x < nact; x++) { + for (int q = 0; q < nthree_; q++) { + if (std::abs(i_bar_x_bar_P_[i]->get(x, q)) >= theta_NB_IAP) { + xbar_ibar_[i].emplace_back(x); + break; + } + } + } + } + } + + if (algorithm == "cavv") { + x_bar_a_bar_P_.resize(nact); + x_bar_a_bar_P_sliced_.resize(nact); + + for (int x = 0; x < nact; x++) { + x_bar_a_bar_P_[x] = std::make_shared("(abar * P)", nvir, nthree_); + x_bar_a_bar_P_[x]->zero(); + } + for (int q = 0; q < nthree_; q++) { + for (int x = 0; x < P_xbar_abar_[q]->rowdim(); x++) { + for (int a = 0; a < P_xbar_abar_[q]->coldim(); a++) { + x_bar_a_bar_P_[xbar_u_p_[nweight][q][x]]->set(a_bar_p_up_[q][a], q, + P_xbar_abar_[q]->get(x, a)); + } + } + } + /// Construct {abar}_xbar. + abar_xbar_.resize(nact); + for (int x = 0; x < nact; x++) { + abar_xbar_[x].clear(); + for (int a = 0; a < nvir; a++) { + for (int q = 0; q < nthree_; q++) { + if (std::abs(x_bar_a_bar_P_[x]->get(a, q)) >= theta_NB_XAP_cavv_) { + abar_xbar_[x].emplace_back(a); + break; + } + } + } + } + } +} + +void LaplaceDSRG::compute_coulomb(const std::string& algorithm, const int& nocc, const int& nact, + const int& nvir) { + /// Construct Z_pq + P_ibar_ = invert_map(ibar_p_, nocc); + Z_pq_->zero(); + + if (algorithm == "ccav") { + P_ibar_x_ = invert_map(ibar_x_p_, nocc); + ZA_pq_->zero(); + } + + for (int i = 0; i < nocc; i++) { + i_bar_a_bar_P_sliced_[i] = + submatrix_rows_and_cols(*i_bar_a_bar_P_[i], abar_ibar_[i], P_ibar_[i]); + psi::SharedMatrix multi_per_i = + psi::linalg::doublet(i_bar_a_bar_P_sliced_[i], i_bar_a_bar_P_sliced_[i], true, false); + for (int p = 0; p < P_ibar_[i].size(); p++) { + for (int q = 0; q < P_ibar_[i].size(); q++) { + int row = P_ibar_[i][p]; + int col = P_ibar_[i][q]; + Z_pq_->add(row, col, multi_per_i->get(p, q)); + } + } + if (algorithm == "ccav") { + i_bar_x_bar_P_sliced_[i] = + submatrix_rows_and_cols(*i_bar_x_bar_P_[i], xbar_ibar_[i], P_ibar_x_[i]); + psi::SharedMatrix multi_per_ix = psi::linalg::doublet( + i_bar_x_bar_P_sliced_[i], i_bar_x_bar_P_sliced_[i], true, false); + for (int px = 0; px < P_ibar_x_[i].size(); px++) { + for (int qx = 0; qx < P_ibar_x_[i].size(); qx++) { + int row_x = P_ibar_x_[i][px]; + int col_x = P_ibar_x_[i][qx]; + ZA_pq_->add(row_x, col_x, multi_per_ix->get(px, qx)); + } + } + } + } + + /// Construct D_pq + psi::SharedMatrix D_pq = psi::linalg::doublet(Z_pq_, C_pq_, false, false); + + if (algorithm == "cavv") { + P_xbar_ = invert_map(xbar_p_, nact); + ZA_pq_->zero(); + for (int x = 0; x < nact; x++) { + x_bar_a_bar_P_sliced_[x] = + submatrix_rows_and_cols(*x_bar_a_bar_P_[x], abar_xbar_[x], P_xbar_[x]); + psi::SharedMatrix multi_per_x = psi::linalg::doublet( + x_bar_a_bar_P_sliced_[x], x_bar_a_bar_P_sliced_[x], true, false); + for (int p = 0; p < P_xbar_[x].size(); p++) { + for (int q = 0; q < P_xbar_[x].size(); q++) { + int row = P_xbar_[x][p]; + int col = P_xbar_[x][q]; + ZA_pq_->add(row, col, multi_per_x->get(p, q)); + } + } + } + } + + /// Compute Coulomb Energy + if (algorithm == "ccvv") { + for (int q = 0; q < nthree_; q++) { + E_J_ -= 2 * D_pq->get_row(0, q)->vector_dot(*D_pq->get_column(0, q)); + } + } else { + psi::SharedMatrix DA_pq = psi::linalg::doublet(ZA_pq_, C_pq_, false, false); + for (int q = 0; q < nthree_; q++) { + E_J_ -= 2 * D_pq->get_row(0, q)->vector_dot(*DA_pq->get_column(0, q)); + } + } +} + +void LaplaceDSRG::compute_exchange(const std::string& algorithm, const int& nocc, const int& nact, + const int& nvir) { + /// Exchange part + B_ia_Q_.resize(nocc); + Q_ia_ = std::make_shared("Q_ia", nocc, nvir); + Q_ia_->zero(); + + if (algorithm == "ccav") { + B_ix_Q_.resize(nocc); + Q_ix_ = std::make_shared("Q_ix", nocc, nact); + Q_ix_->zero(); + } else if (algorithm == "cavv") { + B_xa_Q_.resize(nact); + Q_xa_ = std::make_shared("Q_xa", nact, nvir); + Q_xa_->zero(); + } + + for (int i = 0; i < nocc; i++) { + psi::SharedMatrix sliced_C_pq = submatrix_rows_and_cols(*C_pq_, P_ibar_[i], P_ibar_[i]); + B_ia_Q_[i] = + psi::linalg::doublet(i_bar_a_bar_P_sliced_[i], sliced_C_pq, false, false); /// a*p + + if (algorithm == "ccvv") { + /// Compute (i_bar a_bar|i_bar b_bar) + psi::SharedMatrix iaib_per_i = + psi::linalg::doublet(B_ia_Q_[i], i_bar_a_bar_P_sliced_[i], false, true); + for (int abar = 0; abar < iaib_per_i->rowdim(); abar++) { + E_K_ += iaib_per_i->get_row(0, abar)->vector_dot(*iaib_per_i->get_column(0, abar)); + double Q_value_2 = iaib_per_i->get(abar, abar); + Q_ia_->set(i, abar_ibar_[i][abar], sqrt(Q_value_2)); + } + } else if (algorithm == "ccav") { + /// Cannot do (ia|ib) before screening. The reason is that P_ibar and P_ibar_x are + /// different. + for (int abar = 0; abar < B_ia_Q_[i]->rowdim(); abar++) { + double Q_value_2 = B_ia_Q_[i]->get_row(0, abar)->vector_dot( + *i_bar_a_bar_P_sliced_[i]->get_row(0, abar)); + Q_ia_->set(i, abar_ibar_[i][abar], sqrt(Q_value_2)); + } + psi::SharedMatrix sliced_C_pq_A = + submatrix_rows_and_cols(*C_pq_, P_ibar_x_[i], P_ibar_x_[i]); + B_ix_Q_[i] = + psi::linalg::doublet(i_bar_x_bar_P_sliced_[i], sliced_C_pq_A, false, false); /// x*p + for (int xbar = 0; xbar < B_ix_Q_[i]->rowdim(); xbar++) { + double Q_value_A_2 = B_ix_Q_[i]->get_row(0, xbar)->vector_dot( + *i_bar_x_bar_P_sliced_[i]->get_row(0, xbar)); + Q_ix_->set(i, xbar_ibar_[i][xbar], sqrt(Q_value_A_2)); + } + } else if (algorithm == "cavv") { + for (int abar = 0; abar < B_ia_Q_[i]->rowdim(); abar++) { + double Q_value_2 = B_ia_Q_[i]->get_row(0, abar)->vector_dot( + *i_bar_a_bar_P_sliced_[i]->get_row(0, abar)); + Q_ia_->set(i, abar_ibar_[i][abar], sqrt(Q_value_2)); + } + } + } + + if (algorithm == "cavv") { + for (int x = 0; x < nact; x++) { + psi::SharedMatrix sliced_C_pq_A = + submatrix_rows_and_cols(*C_pq_, P_xbar_[x], P_xbar_[x]); + B_xa_Q_[x] = + psi::linalg::doublet(x_bar_a_bar_P_sliced_[x], sliced_C_pq_A, false, false); + for (int abar = 0; abar < B_xa_Q_[x]->rowdim(); abar++) { + double Q_value_A_2 = B_xa_Q_[x]->get_row(0, abar)->vector_dot( + *x_bar_a_bar_P_sliced_[x]->get_row(0, abar)); + Q_xa_->set(x, abar_xbar_[x][abar], sqrt(Q_value_A_2)); + } + } + } + + /// Start ij-prescreening. + if (algorithm == "ccvv") { + psi::SharedMatrix A_ij = psi::linalg::doublet(Q_ia_, Q_ia_, false, true); + for (int i = 0; i < nocc; i++) { + for (int j = 0; j < i; j++) { + if (std::abs(A_ij->get(i, j)) >= theta_ij_sqrt_) { + vir_intersection_per_ij_.clear(); + aux_intersection_per_ij_.clear(); + std::set_intersection(P_ibar_[i].begin(), P_ibar_[i].end(), P_ibar_[j].begin(), + P_ibar_[j].end(), + std::back_inserter(aux_intersection_per_ij_)); + aux_in_B_i_.clear(); + std::set_intersection(abar_ibar_[i].begin(), abar_ibar_[i].end(), + abar_ibar_[j].begin(), abar_ibar_[j].end(), + std::back_inserter(vir_intersection_per_ij_)); + for (auto aux : aux_intersection_per_ij_) { + int idx_aux_i = + binary_search_recursive(P_ibar_[i], aux, 0, P_ibar_[i].size() - 1); + aux_in_B_i_.emplace_back(idx_aux_i); + } + for (int a_idx = 0; a_idx < vir_intersection_per_ij_.size(); + a_idx++) { // b <= a and j < i + for (int b_idx = 0; b_idx <= a_idx; b_idx++) { + int a = vir_intersection_per_ij_[a_idx]; + int b = vir_intersection_per_ij_[b_idx]; + double Schwarz = Q_ia_->get(i, a) * Q_ia_->get(j, b) * + Q_ia_->get(i, b) * Q_ia_->get(j, a); + if (std::abs(Schwarz) >= theta_schwarz_) { + int a_idx_i = binary_search_recursive(abar_ibar_[i], a, 0, + abar_ibar_[i].size() - 1); + int b_idx_i = binary_search_recursive(abar_ibar_[i], b, 0, + abar_ibar_[i].size() - 1); + + std::vector vec_a_i{a_idx_i}; + std::vector vec_b_i{b_idx_i}; + std::vector vec_a_j{a}; + std::vector vec_b_j{b}; + + psi::SharedMatrix ia = + submatrix_rows_and_cols(*B_ia_Q_[i], vec_a_i, aux_in_B_i_); + psi::SharedMatrix jb = submatrix_rows_and_cols( + *i_bar_a_bar_P_[j], vec_b_j, aux_intersection_per_ij_); + psi::SharedMatrix ib = + submatrix_rows_and_cols(*B_ia_Q_[i], vec_b_i, aux_in_B_i_); + psi::SharedMatrix ja = submatrix_rows_and_cols( + *i_bar_a_bar_P_[j], vec_a_j, aux_intersection_per_ij_); + + psi::SharedMatrix iajb_mat = + psi::linalg::doublet(ia, jb, false, true); + double* iajb = iajb_mat->get_pointer(); + psi::SharedMatrix ibja_mat = + psi::linalg::doublet(ib, ja, false, true); + double* ibja = ibja_mat->get_pointer(); + + if (a == b) { + E_K_ += 2 * (*iajb) * (*ibja); + } else { + E_K_ += 4 * (*iajb) * (*ibja); + } + } + } + } + } + } + } + } else if (algorithm == "ccav") { + psi::SharedMatrix Aa_ij = psi::linalg::doublet(Q_ix_, Q_ix_, false, true); + psi::SharedMatrix A_ij = psi::linalg::doublet(Q_ia_, Q_ia_, false, true); + for (int i = 0; i < nocc; i++) { + for (int j = 0; j <= i; j++) { + if (std::abs(Aa_ij->get(i, j)) * std::abs(A_ij->get(i, j)) >= theta_ij_ccav_) { + aux_intersection_per_ijx_.clear(); /// iajx + aux_intersection_per_ixj_.clear(); /// ixja + std::set_intersection(P_ibar_[i].begin(), P_ibar_[i].end(), + P_ibar_x_[j].begin(), P_ibar_x_[j].end(), + std::back_inserter(aux_intersection_per_ijx_)); + std::set_intersection(P_ibar_[j].begin(), P_ibar_[j].end(), + P_ibar_x_[i].begin(), P_ibar_x_[i].end(), + std::back_inserter(aux_intersection_per_ixj_)); + aux_in_B_i_.clear(); + for (auto aux : aux_intersection_per_ijx_) { + int idx_aux_i = + binary_search_recursive(P_ibar_[i], aux, 0, P_ibar_[i].size() - 1); + aux_in_B_i_.emplace_back(idx_aux_i); + } + aux_in_B_ix_.clear(); + for (auto aux : aux_intersection_per_ixj_) { + int idx_aux_i = + binary_search_recursive(P_ibar_x_[i], aux, 0, P_ibar_x_[i].size() - 1); + aux_in_B_ix_.emplace_back(idx_aux_i); + } + for (int a_idx = 0; a_idx < abar_ibar_[i].size(); a_idx++) { + for (int b_idx = 0; b_idx < xbar_ibar_[i].size(); b_idx++) { + int a = abar_ibar_[i][a_idx]; + int b = xbar_ibar_[i][b_idx]; + double Schwarz = Q_ia_->get(i, a) * Q_ix_->get(j, b) * + Q_ix_->get(i, b) * Q_ia_->get(j, a); + if (std::abs(Schwarz) >= theta_schwarz_ccav_) { + int a_idx_i = binary_search_recursive(abar_ibar_[i], a, 0, + abar_ibar_[i].size() - 1); + int b_idx_i = binary_search_recursive(xbar_ibar_[i], b, 0, + xbar_ibar_[i].size() - 1); + std::vector vec_a_i{a_idx_i}; + std::vector vec_b_i{b_idx_i}; + std::vector vec_a_j{a}; + std::vector vec_b_j{b}; + + psi::SharedMatrix ia = + submatrix_rows_and_cols(*B_ia_Q_[i], vec_a_i, aux_in_B_i_); + psi::SharedMatrix jb = submatrix_rows_and_cols( + *i_bar_x_bar_P_[j], vec_b_j, aux_intersection_per_ijx_); + psi::SharedMatrix ib = + submatrix_rows_and_cols(*B_ix_Q_[i], vec_b_i, aux_in_B_ix_); + psi::SharedMatrix ja = submatrix_rows_and_cols( + *i_bar_a_bar_P_[j], vec_a_j, aux_intersection_per_ixj_); + + psi::SharedMatrix iajb_mat = + psi::linalg::doublet(ia, jb, false, true); + double* iajb = iajb_mat->get_pointer(); + psi::SharedMatrix ibja_mat = + psi::linalg::doublet(ib, ja, false, true); + double* ibja = ibja_mat->get_pointer(); + + if (i == j) { + E_K_ += (*iajb) * (*ibja); + } else { + E_K_ += 2 * (*iajb) * (*ibja); + } + } + } + } + } + } + } + } else if (algorithm == "cavv") { + /// Start ix-prescreening. + psi::SharedMatrix Aa_ix = psi::linalg::doublet(Q_ia_, Q_xa_, false, true); // (occ * act) + for (int i = 0; i < nocc; i++) { + for (int x = 0; x < nact; x++) { + if (std::abs(Aa_ix->get(i, x)) >= theta_ij_sqrt_cavv_) { + vir_intersection_per_ix_.clear(); + aux_intersection_per_ix_.clear(); + std::set_intersection(P_ibar_[i].begin(), P_ibar_[i].end(), P_xbar_[x].begin(), + P_xbar_[x].end(), + std::back_inserter(aux_intersection_per_ix_)); + aux_in_B_i_.clear(); + std::set_intersection(abar_ibar_[i].begin(), abar_ibar_[i].end(), + abar_xbar_[x].begin(), abar_xbar_[x].end(), + std::back_inserter(vir_intersection_per_ix_)); + for (auto aux : aux_intersection_per_ix_) { + int idx_aux_i = + binary_search_recursive(P_ibar_[i], aux, 0, P_ibar_[i].size() - 1); + aux_in_B_i_.emplace_back(idx_aux_i); + } + for (int a_idx = 0; a_idx < vir_intersection_per_ix_.size(); a_idx++) { + for (int b_idx = 0; b_idx <= a_idx; b_idx++) { + int a = vir_intersection_per_ix_[a_idx]; + int b = vir_intersection_per_ix_[b_idx]; + double Schwarz = Q_ia_->get(i, a) * Q_xa_->get(x, b) * + Q_ia_->get(i, b) * Q_xa_->get(x, a); + if (std::abs(Schwarz) >= theta_schwarz_cavv_) { + int a_idx_i = binary_search_recursive(abar_ibar_[i], a, 0, + abar_ibar_[i].size() - 1); + int b_idx_i = binary_search_recursive(abar_ibar_[i], b, 0, + abar_ibar_[i].size() - 1); + std::vector vec_a_i{a_idx_i}; + std::vector vec_b_i{b_idx_i}; + std::vector vec_a_x{a}; + std::vector vec_b_x{b}; + + psi::SharedMatrix ia = + submatrix_rows_and_cols(*B_ia_Q_[i], vec_a_i, aux_in_B_i_); + psi::SharedMatrix xb = submatrix_rows_and_cols( + *x_bar_a_bar_P_[x], vec_b_x, aux_intersection_per_ix_); + psi::SharedMatrix ib = + submatrix_rows_and_cols(*B_ia_Q_[i], vec_b_i, aux_in_B_i_); + psi::SharedMatrix xa = submatrix_rows_and_cols( + *x_bar_a_bar_P_[x], vec_a_x, aux_intersection_per_ix_); + + psi::SharedMatrix iaxb_mat = + psi::linalg::doublet(ia, xb, false, true); + double* iaxb = iaxb_mat->get_pointer(); + psi::SharedMatrix ibxa_mat = + psi::linalg::doublet(ib, xa, false, true); + double* ibxa = ibxa_mat->get_pointer(); + + if (a == b) { + E_K_ += (*iaxb) * (*ibxa); + } else { + E_K_ += 2 * (*iaxb) * (*ibxa); + } + } + } + } + } + } + } + } +} + +double LaplaceDSRG::compute_ccvv() { + local_timer loadingTime; + E_J_ = 0.0; + E_K_ = 0.0; + std::shared_ptr ao_helper = std::make_shared( + Cwfn_, eps_rdocc_, eps_virtual_, laplace_threshold_, nactive_, nfrozen_, vir_tol_); + + prepare_cholesky_coeff(ao_helper, "ccvv"); + Load_int(theta_NB_, "ccvv"); + clear_maps("ccvv"); + + local_timer looptime; + + for (int nweight = 0; nweight < weights_; nweight++) { + int nocc = Occupied_cholesky_[nweight]->coldim(); + int nvir = Virtual_cholesky_[nweight]->coldim(); + fill_maps("ccvv", nweight, nocc, 0, nvir); + compute_coulomb("ccvv", nocc, 0, nvir); + compute_exchange("ccvv", nocc, 0, nvir); + } + return (E_J_ + E_K_); +} + +double LaplaceDSRG::compute_cavv() { + E_J_ = 0.0; + E_K_ = 0.0; + + std::shared_ptr ao_helper = std::make_shared( + Cwfn_, eps_rdocc_, eps_active_, eps_virtual_, laplace_threshold_, nactive_, nfrozen_, true, + vir_tol_); + + prepare_cholesky_coeff(ao_helper, "cavv"); + Load_int(theta_NB_cavv_, "cavv"); + clear_maps("cavv"); + + for (int nweight = 0; nweight < weights_; nweight++) { + int nocc = Occupied_cholesky_[nweight]->coldim(); + int nvir = Virtual_cholesky_[nweight]->coldim(); + int nact = Active_cholesky_[nweight]->coldim(); + fill_maps("cavv", nweight, nocc, nact, nvir); + compute_coulomb("cavv", nocc, nact, nvir); + compute_exchange("cavv", nocc, nact, nvir); + } + return (E_J_ + E_K_); +} + +double LaplaceDSRG::compute_ccav() { + E_J_ = 0.0; + E_K_ = 0.0; + std::shared_ptr ao_helper = std::make_shared( + Cwfn_, eps_rdocc_, eps_active_, eps_virtual_, laplace_threshold_, nactive_, nfrozen_, false, + vir_tol_); + + prepare_cholesky_coeff(ao_helper, "ccav"); + Load_int(theta_NB_ccav_, "ccav"); + clear_maps("ccav"); + + for (int nweight = 0; nweight < weights_; nweight++) { + int nocc = Occupied_cholesky_[nweight]->coldim(); + int nvir = Virtual_cholesky_[nweight]->coldim(); + int nact = Active_cholesky_[nweight]->coldim(); + fill_maps("ccav", nweight, nocc, nact, nvir); + compute_coulomb("ccav", nocc, nact, nvir); + compute_exchange("ccav", nocc, nact, nvir); + } + return (E_J_ + E_K_); +} + +} // namespace forte \ No newline at end of file diff --git a/forte/mrdsrg-spin-integrated/Laplace_dsrg.h b/forte/mrdsrg-spin-integrated/Laplace_dsrg.h new file mode 100644 index 000000000..8428900a4 --- /dev/null +++ b/forte/mrdsrg-spin-integrated/Laplace_dsrg.h @@ -0,0 +1,158 @@ +#ifndef _laplace_dsrg_h_ +#define _laplace_dsrg_h_ + +#include "master_mrdsrg.h" + +using namespace ambit; + +namespace forte { + +class LaplaceDSRG { + public: + LaplaceDSRG(std::shared_ptr options, std::shared_ptr ints, + std::shared_ptr mo_space_info, psi::SharedVector epsilon_rdocc, + psi::SharedVector epsilon_virtual, psi::SharedVector epsilon_active, + psi::SharedMatrix Gamma1, psi::SharedMatrix Eta1); + + ~LaplaceDSRG(); + + void prepare_cholesky_coeff(std::shared_ptr ao_helper_ptr, + const std::string& algorithm); + void Load_int(double& theta_NB, const std::string& algorithm); + void clear_maps(const std::string& algorithm); + void fill_maps(const std::string& algorithm, const int& nweight, const int& nocc, + const int& nact, const int& nvir); + void compute_coulomb(const std::string& algorithm, const int& nocc, const int& nact, + const int& nvir); + void compute_exchange(const std::string& algorithm, const int& nocc, const int& nact, + const int& nvir); + double compute_ccvv(); + double compute_cavv(); + double compute_ccav(); + void print_header(); + size_t vir_start() { return vir_start_; } + + protected: + std::shared_ptr foptions_; + std::shared_ptr ints_; + std::shared_ptr mo_space_info_; + psi::SharedVector eps_rdocc_; + psi::SharedVector eps_virtual_; + psi::SharedVector eps_active_; + psi::SharedMatrix Cwfn_; + psi::SharedMatrix C_pq_; + psi::SharedMatrix S_; + SparseMap i_p_; + + std::shared_ptr primary_; + std::shared_ptr auxiliary_; + + size_t nmo_; + size_t nthree_; + size_t nfrozen_; + size_t ncore_; + size_t nactive_; + size_t nvirtual_; + + /// Thresholds + double Omega_; + double laplace_threshold_; + double vir_tol_; + size_t vir_start_; + size_t weights_; + + /// CCVV + double theta_NB_; + double theta_NB_IAP_; + double theta_ij_; + double theta_schwarz_; + double theta_ij_sqrt_; + + /// CAVV + double theta_NB_cavv_; + double theta_XNB_cavv_; + double theta_NB_IAP_cavv_; + double theta_NB_XAP_cavv_; + double theta_ij_cavv_; + double theta_schwarz_cavv_; + double theta_ij_sqrt_cavv_; + + /// CCAV + double theta_NB_ccav_; + double theta_NB_IAP_ccav_; + double theta_ij_ccav_; + double theta_schwarz_ccav_; + double theta_ij_sqrt_ccav_; + + /// Energy + double E_J_; + double E_K_; + + /// SparseMaps related + SparseMap ao_list_per_q_; + SparseMap i_bar_p_up_; /// Construct [ibar]_p + SparseMap a_bar_p_up_; /// Construct [abar]_p + SparseMap ibar_p_; /// Construct {ibar}_p + SparseMap xbar_p_; /// Construct {xbar}_p + SparseMap ibar_x_p_; /// Construct {ibar_x}_p + SparseMap P_ibar_; /// Construct {P}_ibar. Obtain this by "inversion" of {ibar}_p + SparseMap P_xbar_; /// Construct {P}_xbar + SparseMap P_ibar_x_; /// Construct {p}_ibar_x + SparseMap abar_ibar_; /// Construct {a_bar}_ibar + SparseMap abar_xbar_; /// Construct {a_bar}_xbar + SparseMap xbar_ibar_; /// Construct {x_bar}_ibar + + psi::SharedMatrix N_pu_; + std::vector P_iu_; + std::vector T_ibar_i_list_; + std::vector xbar_u_p_; + std::vector> P_xbar_u_cavv_; + std::vector> P_xbar_u_ccav_; + std::vector P_ibar_u_; /// Construct (P|ibar u) + std::vector P_ibar_abar_; /// Construct (P|ibar abar) + std::vector P_xbar_abar_; /// Construct (P|xbar abar) + std::vector P_ibar_xbar_; /// Construct (P|ibar xbar) + std::vector i_bar_a_bar_P_; /// Construct (ibar abar|P) + std::vector x_bar_a_bar_P_; /// Construct (xbar abar|P) + std::vector i_bar_x_bar_P_; /// Construct (ibar xbar|P) + std::vector i_bar_a_bar_P_sliced_; /// Constrcut Sliced (ibar abar|P); + std::vector x_bar_a_bar_P_sliced_; /// Construct Sliced (xbar abar|P) + std::vector i_bar_x_bar_P_sliced_; /// Construct Sliced (ibar xbar|P) + psi::SharedMatrix Z_pq_; /// Construct Z_pq + psi::SharedMatrix ZA_pq_; /// Construct ZA_pq + std::vector B_ia_Q_; /// Construct B_ia_Q + std::vector B_xa_Q_; /// Construct B_xa_Q + std::vector B_ix_Q_; /// Construct B_ix_Q + psi::SharedMatrix Q_ia_; /// Construc Q_ia square; + psi::SharedMatrix Q_xa_; + psi::SharedMatrix Q_ix_; + std::vector vir_intersection_per_ij_; + std::vector aux_intersection_per_ij_; + std::vector aux_in_B_i_; + std::vector aux_in_B_j_; + std::vector aux_in_B_ix_; + std::vector vir_intersection_per_ix_; + std::vector aux_intersection_per_ix_; + std::vector aux_intersection_per_ijx_; + std::vector aux_intersection_per_ixj_; + + /// CCVV Cholesky + std::vector Occupied_cholesky_; + std::vector Virtual_cholesky_; + psi::SharedMatrix Cholesky_Occ_; + + std::vector Occupied_cholesky_abs_; + std::vector Virtual_cholesky_abs_; + psi::SharedMatrix Cholesky_Occ_abs_; + + /// CAVV Cholesky + std::vector Active_cholesky_; + std::vector Active_cholesky_abs_; + psi::SharedMatrix Gamma1_mat_; + + /// CCAV Cholesky + psi::SharedMatrix Eta1_mat_; +}; +} // namespace forte + +#endif diff --git a/forte/mrdsrg-spin-integrated/three_dsrg_mrpt2.cc b/forte/mrdsrg-spin-integrated/three_dsrg_mrpt2.cc index 3cb26e22c..4a7395fa0 100644 --- a/forte/mrdsrg-spin-integrated/three_dsrg_mrpt2.cc +++ b/forte/mrdsrg-spin-integrated/three_dsrg_mrpt2.cc @@ -54,16 +54,23 @@ #include "psi4/libpsio/psio.h" #include "psi4/libpsio/psio.hpp" #include "psi4/libqt/qt.h" +#include "psi4/libmints/integral.h" +#include "psi4/psifiles.h" +#include "psi4/libfock/jk.h" #include "orbital-helpers/ao_helper.h" #include "helpers/blockedtensorfactory.h" #include "helpers/printing.h" #include "helpers/timer.h" +#include "helpers/helpers.h" #include "fci/fci_solver.h" #include "fci/fci_vector.h" #include "sci/aci.h" #include "three_dsrg_mrpt2.h" +#include "orbital-helpers/Laplace.h" +#include "Laplace_dsrg.h" + using namespace ambit; using namespace psi; @@ -86,6 +93,7 @@ THREE_DSRG_MRPT2::THREE_DSRG_MRPT2(std::shared_ptr rdms, std::shared_ptrget_bool("AO_DSRG_MRPT2")) { - double Eccvv_ao = E_VT2_2_AO_Slow(); - Eccvv = Eccvv_ao; - outfile->Printf("\n Eccvv_ao: %8.10f", Eccvv_ao); - } + // if (foptions_->get_bool("AO_DSRG_MRPT2")) { + // double Eccvv_ao = E_ccvv_lt_ao(); + // Eccvv = Eccvv_ao; + // outfile->Printf("\n Eccvv_ao: %8.10f", Eccvv_ao); + // } if (my_proc == 0) { outfile->Printf("... Done. Timing %15.6f s", ccvv_timer.get()); - outfile->Printf("\n Eccvv: %8.10f", Eccvv); + if (ccvv_algorithm != "LT-DSRG") { + outfile->Printf("\n Eccvv: %8.10f", Eccvv); + } } // double all_e = 0.0; @@ -1813,7 +1826,7 @@ double THREE_DSRG_MRPT2::E_VT2_2_ambit() { BefJKVec.push_back(ambit::Tensor::build(tensor_type_, "BefJK", {nvirtual_, nvirtual_})); RDVec.push_back(ambit::Tensor::build(tensor_type_, "RDVec", {nvirtual_, nvirtual_})); } - bool ao_dsrg_check = foptions_->get_bool("AO_DSRG_MRPT2"); + bool ao_dsrg_check = false; #pragma omp parallel for num_threads(num_threads_) reduction(+ : Ealpha, Ebeta, Emixed) for (size_t m = 0; m < ncore_; ++m) { @@ -1925,6 +1938,8 @@ double THREE_DSRG_MRPT2::E_VT2_2_ambit() { RDVec.push_back(ambit::Tensor::build(tensor_type_, "RD", {nvirtual_, nvirtual_})); } bool ao_dsrg_check = foptions_->get_bool("AO_DSRG_MRPT2"); + // bool ao_dsrg_check = true; + #pragma omp parallel for num_threads(num_threads_) schedule(dynamic) \ reduction(+ : Ealpha, Ebeta, Emixed) shared(Ba, Bb) @@ -2243,6 +2258,290 @@ double THREE_DSRG_MRPT2::E_VT2_2_batch_core() { return (Ealpha + Ebeta + Emixed); } +double THREE_DSRG_MRPT2::E_ccvv_lt_ao() { + bool laplace_cavv = foptions_->get_bool("LAPLACE_CAVV"); + bool laplace_ccav = foptions_->get_bool("LAPLACE_CCAV"); + + double E_cavv = 0.0; + double E_ccav = 0.0; + double E_ccvv = 0.0; + + if (mo_space_info_->nirrep() != 1) { + throw psi::PSIEXCEPTION("LT-DSRG-MRPT2 does not work with symmetry."); + } + if (integral_type_ != DiskDF) { + throw psi::PSIEXCEPTION("LT-DSRG-MRPT2 requires DiskDF."); + } + + ncore_ = core_mos_.size(); + nactive_ = actv_mos_.size(); + nvirtual_ = virt_mos_.size(); + nthree_ = ints_->nthree(); + int nfrozen = mo_space_info_->dimension("FROZEN").sum(); + + psi::SharedVector epsilon_rdocc(new psi::Vector("EPS_RDOCC", ncore_)); + psi::SharedVector epsilon_virtual(new psi::Vector("EPS_VIRTUAL", nvirtual_)); + psi::SharedVector epsilon_active(new psi::Vector("EPS_ACTIVE", nactive_)); + + int core_count = 0; + for (auto m : core_mos_) { + epsilon_rdocc->set(core_count, Fa_[m]); + core_count++; + } + + int virtual_count = 0; + for (auto e : virt_mos_) { + epsilon_virtual->set(virtual_count, Fa_[e]); + virtual_count++; + } + + int actv_count = 0; + for (auto u : actv_mos_) { + epsilon_active->set(actv_count, Fa_[u]); + actv_count++; + } + + epsilon_active->print(); + + ambit::Tensor Gamma1a = Gamma1_.block("aa").clone(); + Gamma1a("pq") = rdms_->SF_L1()("pq"); + + ambit::Tensor Eta1a = Eta1_.block("aa").clone(); + Eta1a.iterate( + [&](const std::vector& i, double& value) { value = i[0] == i[1] ? 2.0 : 0.0; }); + Eta1a("pq") -= Gamma1a("pq"); + + psi::SharedMatrix Gamma1_mat = ambit_to_matrix(Gamma1a); + psi::SharedMatrix Eta1_mat = tensor_to_matrix(Eta1a); + + auto LaplaceDSRG = + std::make_shared(foptions_, ints_, mo_space_info_, epsilon_rdocc, + epsilon_virtual, epsilon_active, Gamma1_mat, Eta1_mat); + + size_t vir_start = 0; + size_t vir_start_ccav = 0; + size_t vir_start_cavv = 0; + + local_timer timer1; + E_ccvv = LaplaceDSRG->compute_ccvv(); + outfile->Printf("\n\n LAPLACE: ccvv takes %8.8f", timer1.get()); + outfile->Printf("\n LAPLACE: E_ccvv %8.10f", E_ccvv); + + vir_start = LaplaceDSRG->vir_start(); + + if (laplace_cavv) { + local_timer timer2; + E_cavv = LaplaceDSRG->compute_cavv(); + outfile->Printf("\n\n LAPLACE: cavv takes %8.8f", timer2.get()); + outfile->Printf("\n LAPLACE: E_cavv %8.10f", E_cavv); + vir_start_cavv = LaplaceDSRG->vir_start(); + } + + if (laplace_ccav) { + local_timer timer3; + E_ccav = LaplaceDSRG->compute_ccav(); + outfile->Printf("\n\n LAPLACE: ccav takes %8.8f", timer3.get()); + outfile->Printf("\n LAPLACE: E_ccav %8.10f", E_ccav); + vir_start_ccav = LaplaceDSRG->vir_start(); + } + + LaplaceDSRG.reset(); + + if (vir_start != 0) { + size_t vir_laplace = nvirtual_ - vir_start; + + std::vector virt_mos(virt_mos_.begin(), virt_mos_.begin() + vir_start); + std::vector virt_mos_laplace(virt_mos_.begin() + vir_start, virt_mos_.end()); + + local_timer ccvvTimer; + + ambit::Tensor B_QME = ints_->three_integral_block(aux_mos_, core_mos_, virt_mos_laplace); + ambit::Tensor B_QMe = ints_->three_integral_block(aux_mos_, core_mos_, virt_mos); + + ambit::Tensor V_MeNf = + ambit::Tensor::build(tensor_type_, "V", {ncore_, vir_start, ncore_, vir_start}); + V_MeNf("m,e,n,f") = B_QMe("Q, m, e") * B_QMe("Q, n, f"); + + ambit::Tensor V_MeNF = + ambit::Tensor::build(tensor_type_, "V", {ncore_, vir_start, ncore_, vir_laplace}); + V_MeNF("m,e,n,f") = B_QMe("Q, m, e") * B_QME("Q, n, f"); + + ambit::Tensor T_MeNf = + ambit::Tensor::build(tensor_type_, "T2", {ncore_, vir_start, ncore_, vir_start}); + T_MeNf.data() = V_MeNf.data(); + + ambit::Tensor T_MeNF = + ambit::Tensor::build(tensor_type_, "T2", {ncore_, vir_start, ncore_, vir_laplace}); + T_MeNF.data() = V_MeNF.data(); + + V_MeNf.iterate([&](const std::vector& i, double& value) { + double Exp = Fa_[virt_mos[i[3]]] + Fa_[virt_mos[i[1]]] - Fa_[core_mos_[i[2]]] - + Fa_[core_mos_[i[0]]]; + double D = -1.0 * (Exp); + value = value + value * dsrg_source_->compute_renormalized(Exp); + T_MeNf.data()[i[0] * vir_start * ncore_ * vir_start + i[1] * ncore_ * vir_start + + i[2] * vir_start + i[3]] *= + dsrg_source_->compute_renormalized_denominator(D); + }); + + V_MeNF.iterate([&](const std::vector& i, double& value) { + double Exp = Fa_[virt_mos_laplace[i[3]]] + Fa_[virt_mos[i[1]]] - Fa_[core_mos_[i[2]]] - + Fa_[core_mos_[i[0]]]; + double D = -1.0 * (Exp); + value = value + value * dsrg_source_->compute_renormalized(Exp); + T_MeNF.data()[i[0] * vir_start * ncore_ * vir_laplace + i[1] * ncore_ * vir_laplace + + i[2] * vir_laplace + i[3]] *= + dsrg_source_->compute_renormalized_denominator(D); + }); + + E_ccvv += 2 * V_MeNf("m,e,n,f") * T_MeNf("m,e,n,f") - V_MeNf("m,e,n,f") * T_MeNf("m,f,n,e"); + E_ccvv += + 4 * V_MeNF("m,e,n,f") * T_MeNF("m,e,n,f") - 2 * V_MeNF("m,e,n,f") * T_MeNF("n,e,m,f"); + + if (print_ > 0) { + outfile->Printf("\n\n Leftover CCVV computation takes %8.8f", ccvvTimer.get()); + outfile->Printf("\n Corrected LAPLACE: E_ccvv %8.10f", E_ccvv); + } + } + + if (vir_start_cavv != 0) { + size_t vir_laplace_cavv = nvirtual_ - vir_start_cavv; + + std::vector virt_mos_cavv(virt_mos_.begin(), virt_mos_.begin() + vir_start_cavv); + std::vector virt_mos_laplace_cavv(virt_mos_.begin() + vir_start_cavv, + virt_mos_.end()); + + local_timer cavvTimer; + + ambit::Tensor temp_aa = ambit::Tensor::build(tensor_type_, "TEMPaa", {nactive_, nactive_}); + + ambit::Tensor B_QUF = + ints_->three_integral_block(aux_mos_, actv_mos_, virt_mos_laplace_cavv); + ambit::Tensor B_QUf = ints_->three_integral_block(aux_mos_, actv_mos_, virt_mos_cavv); + ambit::Tensor B_QME = + ints_->three_integral_block(aux_mos_, core_mos_, virt_mos_laplace_cavv); + ambit::Tensor B_QMe = ints_->three_integral_block(aux_mos_, core_mos_, virt_mos_cavv); + + ambit::Tensor V_MeUf = ambit::Tensor::build( + tensor_type_, "V", {ncore_, vir_start_cavv, nactive_, vir_start_cavv}); + V_MeUf("m,e,u,f") += B_QMe("Q, m, e") * B_QUf("Q, u, f"); + + ambit::Tensor V_MeUF_laplace = ambit::Tensor::build( + tensor_type_, "V", {ncore_, vir_start_cavv, nactive_, vir_laplace_cavv}); + ambit::Tensor temp = ambit::Tensor::build( + tensor_type_, "temp", {ncore_, vir_start_cavv, nactive_, vir_laplace_cavv}); + V_MeUF_laplace("m,e,u,f") += B_QMe("Q, m, e") * B_QUF("Q, u, f"); + + ambit::Tensor V_MEUf_laplace = ambit::Tensor::build( + tensor_type_, "V", {ncore_, vir_laplace_cavv, nactive_, vir_start_cavv}); + ambit::Tensor temp2 = ambit::Tensor::build( + tensor_type_, "temp2", {ncore_, vir_laplace_cavv, nactive_, vir_start_cavv}); + V_MEUf_laplace("m,e,u,f") += B_QME("Q, m, e") * B_QUf("Q, u, f"); + + ambit::Tensor T_MeUf = ambit::Tensor::build( + tensor_type_, "T2", {ncore_, vir_start_cavv, nactive_, vir_start_cavv}); + T_MeUf.data() = V_MeUf.data(); + + ambit::Tensor T_MeUF_laplace = ambit::Tensor::build( + tensor_type_, "T2", {ncore_, vir_start_cavv, nactive_, vir_laplace_cavv}); + T_MeUF_laplace.data() = V_MeUF_laplace.data(); + + ambit::Tensor T_MEUf_laplace = ambit::Tensor::build( + tensor_type_, "T2", {ncore_, vir_laplace_cavv, nactive_, vir_start_cavv}); + T_MEUf_laplace.data() = V_MEUf_laplace.data(); + + V_MeUf.iterate([&](const std::vector& i, double& value) { + double Exp = Fa_[virt_mos_cavv[i[3]]] + Fa_[virt_mos_cavv[i[1]]] - + Fa_[actv_mos_[i[2]]] - Fa_[core_mos_[i[0]]]; + double D = -1.0 * (Exp); + value = value + value * dsrg_source_->compute_renormalized(Exp); + T_MeUf.data()[i[0] * vir_start_cavv * nactive_ * vir_start_cavv + + i[1] * nactive_ * vir_start_cavv + i[2] * vir_start_cavv + i[3]] *= + dsrg_source_->compute_renormalized_denominator(D); + }); + + V_MeUF_laplace.iterate([&](const std::vector& i, double& value) { + double Exp = Fa_[virt_mos_laplace_cavv[i[3]]] + Fa_[virt_mos_cavv[i[1]]] - + Fa_[actv_mos_[i[2]]] - Fa_[core_mos_[i[0]]]; + double D = -1.0 * (Exp); + value = value + value * dsrg_source_->compute_renormalized(Exp); + T_MeUF_laplace + .data()[i[0] * vir_start_cavv * nactive_ * vir_laplace_cavv + + i[1] * nactive_ * vir_laplace_cavv + i[2] * vir_laplace_cavv + i[3]] *= + dsrg_source_->compute_renormalized_denominator(D); + }); + + V_MEUf_laplace.iterate([&](const std::vector& i, double& value) { + double Exp = Fa_[virt_mos_cavv[i[3]]] + Fa_[virt_mos_laplace_cavv[i[1]]] - + Fa_[actv_mos_[i[2]]] - Fa_[core_mos_[i[0]]]; + double D = -1.0 * (Exp); + value = value + value * dsrg_source_->compute_renormalized(Exp); + T_MEUf_laplace + .data()[i[0] * vir_laplace_cavv * nactive_ * vir_start_cavv + + i[1] * nactive_ * vir_start_cavv + i[2] * vir_start_cavv + i[3]] *= + dsrg_source_->compute_renormalized_denominator(D); + }); + + temp_aa("u,v") += 2 * V_MeUf("m,e,u,f") * T_MeUf("m,e,v,f"); + temp_aa("u,v") -= V_MeUf("m,e,u,f") * T_MeUf("m,f,v,e"); + + temp("m,e,u,f") = 2 * V_MeUF_laplace("m,e,u,f") - V_MEUf_laplace("m,f,u,e"); + + temp2("m,e,u,f") = 2 * V_MEUf_laplace("m,e,u,f") - V_MeUF_laplace("m,f,u,e"); + + temp_aa("u,v") += temp("m,e,u,f") * T_MeUF_laplace("m,e,v,f"); + temp_aa("u,v") += temp2("m,e,u,f") * T_MEUf_laplace("m,e,v,f"); + + E_cavv += temp_aa("u,v") * Gamma1a("u,v"); + + if (print_ > 0) { + outfile->Printf("\n\n Leftover CAVV computation takes %8.8f", cavvTimer.get()); + outfile->Printf("\n Corrected LAPLACE: E_cavv %8.10f", E_cavv); + } + } + + if (vir_start_ccav != 0) { + std::vector virt_mos_ccav(virt_mos_.begin(), virt_mos_.begin() + vir_start_ccav); + + local_timer ccavTimer; + + ambit::Tensor temp_aa = ambit::Tensor::build(tensor_type_, "TEMPaa", {nactive_, nactive_}); + + ambit::Tensor B_QMU = ints_->three_integral_block(aux_mos_, core_mos_, actv_mos_); + ambit::Tensor B_QMe = ints_->three_integral_block(aux_mos_, core_mos_, virt_mos_ccav); + + ambit::Tensor V_MUMe = + ambit::Tensor::build(tensor_type_, "V", {ncore_, nactive_, ncore_, vir_start_ccav}); + V_MUMe("m,u,n,f") += B_QMU("Q, m, u") * B_QMe("Q, n, f"); + + ambit::Tensor T_MUMe = + ambit::Tensor::build(tensor_type_, "T2", {ncore_, nactive_, ncore_, vir_start_ccav}); + T_MUMe.data() = V_MUMe.data(); + V_MUMe.iterate([&](const std::vector& i, double& value) { + double Exp = Fa_[virt_mos_ccav[i[3]]] + Fa_[actv_mos_[i[1]]] - Fa_[core_mos_[i[2]]] - + Fa_[core_mos_[i[0]]]; + double D = -1.0 * (Exp); + value = value + value * dsrg_source_->compute_renormalized(Exp); + T_MUMe.data()[i[0] * nactive_ * ncore_ * vir_start_ccav + + i[1] * ncore_ * vir_start_ccav + i[2] * vir_start_ccav + i[3]] *= + dsrg_source_->compute_renormalized_denominator(D); + }); + + temp_aa("u,v") += 2 * V_MUMe("m,u,n,f") * T_MUMe("m,v,n,f"); + temp_aa("u,v") -= V_MUMe("m,u,n,f") * T_MUMe("n,v,m,f"); + + E_ccav += temp_aa("u,v") * Eta1a("u,v"); + + if (print_ > 0) { + outfile->Printf("\n\n Leftover CCAV computation takes %8.8f", ccavTimer.get()); + outfile->Printf("\n Corrected LAPLACE: E_ccav %8.10f", E_ccav); + } + } + + return (E_ccvv + E_cavv + E_ccav); +} + +/* double THREE_DSRG_MRPT2::E_VT2_2_AO_Slow() { /// E_{DSRG} -> Conventional basis @@ -2360,6 +2659,8 @@ double THREE_DSRG_MRPT2::E_VT2_2_AO_Slow() { return (0.25 * Ealpha + 0.25 * Ebeta + Emixed); } +*/ + double THREE_DSRG_MRPT2::E_VT2_2_batch_virtual() { bool debug_print = foptions_->get_bool("DSRG_MRPT2_DEBUG"); double Ealpha = 0.0; @@ -2666,293 +2967,600 @@ double THREE_DSRG_MRPT2::E_VT2_2_core() { double THREE_DSRG_MRPT2::E_VT2_2_one_active() { double Eccva = 0; double Eacvv = 0; + bool laplace_cavv = foptions_->get_bool("LAPLACE_CAVV"); + bool laplace_ccav = foptions_->get_bool("LAPLACE_CCAV"); + if (laplace_cavv && laplace_ccav) { + outfile->Printf( + "\n CAVV and CCAV have not been calculated. Laplace algorithm will be used."); + return (Eacvv + Eccva); + } + int nthread = 1; #ifdef _OPENMP nthread = omp_get_max_threads(); // thread = omp_get_thread_num(); #endif - /// This block of code assumes that ThreeIntegral are not stored as a member variable. - /// Requires the reading from aptei_block which makes code - ambit::Tensor Gamma1_aa = Gamma1_.block("aa"); - ambit::Tensor Gamma1_AA = Gamma1_.block("AA"); - - std::vector Bm_Qe; - std::vector Bm_Qf; + // /// This block of code assumes that ThreeIntegral are not stored as a member variable. + // /// Requires the reading from aptei_block which makes code + // ambit::Tensor Gamma1_aa = Gamma1_.block("aa"); + // ambit::Tensor Gamma1_AA = Gamma1_.block("AA"); - std::vector Vefu; - std::vector Tefv; - std::vector tempTaa; - std::vector tempTAA; + // std::vector Bm_Qe; + // std::vector Bm_Qf; - local_timer ccvaTimer; - for (int thread = 0; thread < nthread; thread++) { - Bm_Qe.push_back(ambit::Tensor::build(tensor_type_, "BemQ", {nthree_, nvirtual_})); - Bm_Qf.push_back(ambit::Tensor::build(tensor_type_, "Bmq", {nthree_, nvirtual_})); + // std::vector Vefu; + // std::vector Tefv; + // std::vector tempTaa; + // std::vector tempTAA; - Vefu.push_back( - ambit::Tensor::build(tensor_type_, "muJK", {nvirtual_, nvirtual_, nactive_})); - Tefv.push_back(ambit::Tensor::build(tensor_type_, "T2", {nvirtual_, nvirtual_, nactive_})); - - tempTaa.push_back(ambit::Tensor::build(tensor_type_, "TEMPaa", {nactive_, nactive_})); - tempTAA.push_back(ambit::Tensor::build(tensor_type_, "TEMPAA", {nactive_, nactive_})); - } - // ambit::Tensor BemQ = ints_->three_integral_block(naux, acore_mos_, - // avirt_mos_); - // ambit::Tensor BeuQ = ints_->three_integral_block(naux, aactv_mos_, - // avirt_mos_); + if (laplace_cavv) { + outfile->Printf("\n CAVV has not been calculated. Laplace algorithm will be used."); + } else { + /// This block of code assumes that ThreeIntegral are not stored as a member variable. + /// Requires the reading from aptei_block which makes code + ambit::Tensor Gamma1_aa = Gamma1_.block("aa"); + ambit::Tensor Gamma1_AA = Gamma1_.block("AA"); + + std::vector Bm_Qe; + std::vector Bm_Qf; + + std::vector Vefu; + std::vector Tefv; + std::vector tempTaa; + std::vector tempTAA; + local_timer cavvTimer; + for (int thread = 0; thread < nthread; thread++) { + Bm_Qe.push_back(ambit::Tensor::build(tensor_type_, "BemQ", {nthree_, nvirtual_})); + Bm_Qf.push_back(ambit::Tensor::build(tensor_type_, "Bmq", {nthree_, nvirtual_})); + + Vefu.push_back( + ambit::Tensor::build(tensor_type_, "muJK", {nvirtual_, nvirtual_, nactive_})); + Tefv.push_back( + ambit::Tensor::build(tensor_type_, "T2", {nvirtual_, nvirtual_, nactive_})); + + tempTaa.push_back(ambit::Tensor::build(tensor_type_, "TEMPaa", {nactive_, nactive_})); + tempTAA.push_back(ambit::Tensor::build(tensor_type_, "TEMPAA", {nactive_, nactive_})); + } + // ambit::Tensor BemQ = ints_->three_integral_block(naux, acore_mos_, + // avirt_mos_); + // ambit::Tensor BeuQ = ints_->three_integral_block(naux, aactv_mos_, + // avirt_mos_); - // Loop over e and f to compute V + // Loop over e and f to compute V - ambit::Tensor BeuQ = ints_->three_integral_block(aux_mos_, virt_mos_, actv_mos_); + ambit::Tensor BeuQ = ints_->three_integral_block(aux_mos_, virt_mos_, actv_mos_); -// std::vector& BemQ_data = BemQ.data(); + // std::vector& BemQ_data = BemQ.data(); -// I think this loop is typically too small to allow efficient use of -// OpenMP. Should probably test this assumption. + // I think this loop is typically too small to allow efficient use of + // OpenMP. Should probably test this assumption. #pragma omp parallel for num_threads(num_threads_) - for (size_t m = 0; m < ncore_; m++) { - int thread = 0; + for (size_t m = 0; m < ncore_; m++) { + int thread = 0; #ifdef _OPENMP - thread = omp_get_thread_num(); + thread = omp_get_thread_num(); #endif - size_t ma = core_mos_[m]; + size_t ma = core_mos_[m]; -// V[efu]_m = B_{em}^Q * B_{fu}^Q - B_{eu}^Q B_{fm}^Q -// V[efu]_m = V[efmu] + V[efmu] * exp[efmu] -// T2["mvef"] = V["mvef"] * D["mvef"] -// temp["uv"] = V * T2 + // V[efu]_m = B_{em}^Q * B_{fu}^Q - B_{eu}^Q B_{fm}^Q + // V[efu]_m = V[efmu] + V[efmu] * exp[efmu] + // T2["mvef"] = V["mvef"] * D["mvef"] + // temp["uv"] = V * T2 #pragma omp critical - { Bm_Qe[thread] = ints_->three_integral_block_two_index(aux_mos_, ma, virt_mos_); } + { Bm_Qe[thread] = ints_->three_integral_block_two_index(aux_mos_, ma, virt_mos_); } - Vefu[thread]("e, f, u") = Bm_Qe[thread]("Q, e") * BeuQ("Q, f, u"); - Vefu[thread]("e, f, u") -= BeuQ("Q, e, u") * Bm_Qe[thread]("Q, f"); + Vefu[thread]("e, f, u") = Bm_Qe[thread]("Q, e") * BeuQ("Q, f, u"); + Vefu[thread]("e, f, u") -= BeuQ("Q, e, u") * Bm_Qe[thread]("Q, f"); - // E = V["efmu"] (1 + Exp(-s * D^{ef}_{mu}) * V^{mv}_{ef} * - // Denom^{mv}_{ef} - Tefv[thread].data() = Vefu[thread].data(); + // E = V["efmu"] (1 + Exp(-s * D^{ef}_{mu}) * V^{mv}_{ef} * + // Denom^{mv}_{ef} + Tefv[thread].data() = Vefu[thread].data(); + + std::vector& T_mv_data = Tefv[thread].data(); + Vefu[thread].iterate([&](const std::vector& i, double& value) { + double Exp = + Fa_[virt_mos_[i[0]]] + Fa_[virt_mos_[i[1]]] - Fa_[actv_mos_[i[2]]] - Fa_[ma]; + double D = -1.0 * (Fa_[virt_mos_[i[0]]] + Fa_[virt_mos_[i[1]]] - + Fa_[actv_mos_[i[2]]] - Fa_[ma]); + value = value + value * dsrg_source_->compute_renormalized(Exp); + T_mv_data[i[0] * nvirtual_ * nactive_ + i[1] * nactive_ + i[2]] *= + dsrg_source_->compute_renormalized_denominator(D); + }); - std::vector& T_mv_data = Tefv[thread].data(); - Vefu[thread].iterate([&](const std::vector& i, double& value) { - double Exp = - Fa_[virt_mos_[i[0]]] + Fa_[virt_mos_[i[1]]] - Fa_[actv_mos_[i[2]]] - Fa_[ma]; - double D = -1.0 * (Fa_[virt_mos_[i[0]]] + Fa_[virt_mos_[i[1]]] - Fa_[actv_mos_[i[2]]] - - Fa_[ma]); - value = value + value * dsrg_source_->compute_renormalized(Exp); - T_mv_data[i[0] * nvirtual_ * nactive_ + i[1] * nactive_ + i[2]] *= - dsrg_source_->compute_renormalized_denominator(D); - }); + // T_mv[thread].iterate([&](const std::vector& i,double& + // value){ + // double D = Fa_[aactv_mos_[i[1]]] + Fa_[acore_mos_[i[0]]] - + // Fa_[ea] - Fa_[fa]; + // value = value * + // dsrg_source_->compute_renormalized_denominator(D);}); - // T_mv[thread].iterate([&](const std::vector& i,double& - // value){ - // double D = Fa_[aactv_mos_[i[1]]] + Fa_[acore_mos_[i[0]]] - - // Fa_[ea] - Fa_[fa]; - // value = value * - // dsrg_source_->compute_renormalized_denominator(D);}); - - tempTaa[thread]("u,v") += 0.5 * Vefu[thread]("e, f, u") * Tefv[thread]("e, f, v"); - Vefu[thread].zero(); - Tefv[thread].zero(); - - Vefu[thread].zero(); - Vefu[thread]("e, f, u") = Bm_Qe[thread]("Q, e") * BeuQ("Q, f, u"); - - // E = V["efmu"] (1 + Exp(-s * D^{ef}_{mu}) * V^{mv}_{ef} * - // Denom^{mv}_{ef} - Tefv[thread].data() = Vefu[thread].data(); - T_mv_data = Tefv[thread].data(); - T_mv_data = Tefv[thread].data(); - Vefu[thread].iterate([&](const std::vector& i, double& value) { - double Exp = - Fa_[virt_mos_[i[0]]] + Fb_[virt_mos_[i[1]]] - Fa_[actv_mos_[i[2]]] - Fb_[ma]; - double D = -1.0 * (Fa_[virt_mos_[i[0]]] + Fb_[virt_mos_[i[1]]] - Fa_[actv_mos_[i[2]]] - - Fb_[ma]); - value = value + value * dsrg_source_->compute_renormalized(Exp); - T_mv_data[i[0] * nvirtual_ * nactive_ + i[1] * nactive_ + i[2]] *= - dsrg_source_->compute_renormalized_denominator(D); - }); + tempTaa[thread]("u,v") += 0.5 * Vefu[thread]("e, f, u") * Tefv[thread]("e, f, v"); + Vefu[thread].zero(); + Tefv[thread].zero(); - // T_mv[thread].iterate([&](const std::vector& i,double& - // value){ - // double D = Fa_[aactv_mos_[i[1]]] + Fa_[acore_mos_[i[0]]] - - // Fa_[ea] - Fa_[fa]; - // value = value * - // dsrg_source_->compute_renormalized_denominator(D);}); - - tempTAA[thread]("vu") += Vefu[thread]("e, f, u") * Tefv[thread]("e,f, v"); - tempTaa[thread]("vu") += Vefu[thread]("e,f, u") * Tefv[thread]("e,f, v"); - Vefu[thread].zero(); - Tefv[thread].zero(); - - Vefu[thread]("e, f, u") = Bm_Qe[thread]("Q, e") * BeuQ("Q, f, u"); - Vefu[thread]("e, f, u") -= BeuQ("Q, e, u") * Bm_Qe[thread]("Q, f"); - - // E = V["efmu"] (1 + Exp(-s * D^{ef}_{mu}) * V^{mv}_{ef} * - // Denom^{mv}_{ef} - - Tefv[thread].data() = Vefu[thread].data(); - T_mv_data = Tefv[thread].data(); - - T_mv_data = Tefv[thread].data(); - Vefu[thread].iterate([&](const std::vector& i, double& value) { - double Exp = - Fa_[virt_mos_[i[0]]] + Fb_[virt_mos_[i[1]]] - Fb_[actv_mos_[i[2]]] - Fb_[ma]; - double D = -1.0 * (Fa_[virt_mos_[i[0]]] + Fa_[virt_mos_[i[1]]] - Fb_[actv_mos_[i[2]]] - - Fb_[ma]); - value = value + value * dsrg_source_->compute_renormalized(Exp); - T_mv_data[i[0] * nvirtual_ * nactive_ + i[1] * nactive_ + i[2]] *= - dsrg_source_->compute_renormalized_denominator(D); - }); + Vefu[thread].zero(); + Vefu[thread]("e, f, u") = Bm_Qe[thread]("Q, e") * BeuQ("Q, f, u"); - tempTaa[thread]("u,v") += 0.5 * Vefu[thread]("e, f, u") * Tefv[thread]("e, f, v"); - } + // E = V["efmu"] (1 + Exp(-s * D^{ef}_{mu}) * V^{mv}_{ef} * + // Denom^{mv}_{ef} + Tefv[thread].data() = Vefu[thread].data(); + T_mv_data = Tefv[thread].data(); + T_mv_data = Tefv[thread].data(); + Vefu[thread].iterate([&](const std::vector& i, double& value) { + double Exp = + Fa_[virt_mos_[i[0]]] + Fb_[virt_mos_[i[1]]] - Fa_[actv_mos_[i[2]]] - Fb_[ma]; + double D = -1.0 * (Fa_[virt_mos_[i[0]]] + Fb_[virt_mos_[i[1]]] - + Fa_[actv_mos_[i[2]]] - Fb_[ma]); + value = value + value * dsrg_source_->compute_renormalized(Exp); + T_mv_data[i[0] * nvirtual_ * nactive_ + i[1] * nactive_ + i[2]] *= + dsrg_source_->compute_renormalized_denominator(D); + }); - ambit::Tensor tempTAA_all = - ambit::Tensor::build(tensor_type_, "tempTAA_all", {nactive_, nactive_}); - ambit::Tensor tempTaa_all = - ambit::Tensor::build(tensor_type_, "tempTaa_all", {nactive_, nactive_}); - for (int thread = 0; thread < nthread; thread++) { - tempTAA_all("v, u") += tempTAA[thread]("v, u"); - tempTaa_all("v, u") += tempTaa[thread]("v, u"); - } + // T_mv[thread].iterate([&](const std::vector& i,double& + // value){ + // double D = Fa_[aactv_mos_[i[1]]] + Fa_[acore_mos_[i[0]]] - + // Fa_[ea] - Fa_[fa]; + // value = value * + // dsrg_source_->compute_renormalized_denominator(D);}); - Eacvv += tempTAA_all("v,u") * Gamma1_AA("v,u"); - Eacvv += tempTaa_all("v,u") * Gamma1_aa("v,u"); + tempTAA[thread]("vu") += Vefu[thread]("e, f, u") * Tefv[thread]("e,f, v"); + tempTaa[thread]("vu") += Vefu[thread]("e,f, u") * Tefv[thread]("e,f, v"); + Vefu[thread].zero(); + Tefv[thread].zero(); - if (print_ > 1) { - outfile->Printf("\n\n CAVV computation takes %8.8f", ccvaTimer.get()); - } + Vefu[thread]("e, f, u") = Bm_Qe[thread]("Q, e") * BeuQ("Q, f, u"); + Vefu[thread]("e, f, u") -= BeuQ("Q, e, u") * Bm_Qe[thread]("Q, f"); + + // E = V["efmu"] (1 + Exp(-s * D^{ef}_{mu}) * V^{mv}_{ef} * + // Denom^{mv}_{ef} + + Tefv[thread].data() = Vefu[thread].data(); + T_mv_data = Tefv[thread].data(); - std::vector Bm_vQ; - std::vector Bn_eQ; - std::vector Bm_eQ; - std::vector Bn_vQ; + T_mv_data = Tefv[thread].data(); + Vefu[thread].iterate([&](const std::vector& i, double& value) { + double Exp = + Fa_[virt_mos_[i[0]]] + Fb_[virt_mos_[i[1]]] - Fb_[actv_mos_[i[2]]] - Fb_[ma]; + double D = -1.0 * (Fa_[virt_mos_[i[0]]] + Fa_[virt_mos_[i[1]]] - + Fb_[actv_mos_[i[2]]] - Fb_[ma]); + value = value + value * dsrg_source_->compute_renormalized(Exp); + T_mv_data[i[0] * nvirtual_ * nactive_ + i[1] * nactive_ + i[2]] *= + dsrg_source_->compute_renormalized_denominator(D); + }); - std::vector V_eu; - std::vector T_ev; - std::vector tempTaa_e; - std::vector tempTAA_e; + tempTaa[thread]("u,v") += 0.5 * Vefu[thread]("e, f, u") * Tefv[thread]("e, f, v"); + } - ambit::Tensor BmvQ = ints_->three_integral_block(aux_mos_, core_mos_, actv_mos_); - ambit::Tensor BmvQ_swapped = - ambit::Tensor::build(tensor_type_, "Bm_vQ", {ncore_, nthree_, nactive_}); - BmvQ_swapped("m, Q, u") = BmvQ("Q, m, u"); - local_timer cavvTimer; - for (int thread = 0; thread < nthread; thread++) { - Bm_vQ.push_back(ambit::Tensor::build(tensor_type_, "BemQ", {nthree_, nactive_})); - Bn_eQ.push_back(ambit::Tensor::build(tensor_type_, "Bf_uQ", {nthree_, nvirtual_})); - Bm_eQ.push_back(ambit::Tensor::build(tensor_type_, "Bmq", {nthree_, nvirtual_})); - Bn_vQ.push_back(ambit::Tensor::build(tensor_type_, "Bmq", {nthree_, nactive_})); + ambit::Tensor tempTAA_all = + ambit::Tensor::build(tensor_type_, "tempTAA_all", {nactive_, nactive_}); + ambit::Tensor tempTaa_all = + ambit::Tensor::build(tensor_type_, "tempTaa_all", {nactive_, nactive_}); + for (int thread = 0; thread < nthread; thread++) { + tempTAA_all("v, u") += tempTAA[thread]("v, u"); + tempTaa_all("v, u") += tempTaa[thread]("v, u"); + } - V_eu.push_back(ambit::Tensor::build(tensor_type_, "muJK", {nvirtual_, nactive_})); - T_ev.push_back(ambit::Tensor::build(tensor_type_, "T2", {nvirtual_, nactive_})); + Eacvv += tempTAA_all("v,u") * Gamma1_AA("v,u"); + Eacvv += tempTaa_all("v,u") * Gamma1_aa("v,u"); - tempTaa_e.push_back(ambit::Tensor::build(tensor_type_, "TEMPaa", {nactive_, nactive_})); - tempTAA_e.push_back(ambit::Tensor::build(tensor_type_, "TEMPAA", {nactive_, nactive_})); + if (print_ > 0) { + outfile->Printf("\n\n CAVV computation takes %8.8f", cavvTimer.get()); + } } - ambit::Tensor Eta1_aa = Eta1_.block("aa"); - ambit::Tensor Eta1_AA = Eta1_.block("AA"); + // local_timer cavvTimer; + // for (int thread = 0; thread < nthread; thread++) { + // Bm_Qe.push_back(ambit::Tensor::build(tensor_type_, "BemQ", {nthree_, nvirtual_})); + // Bm_Qf.push_back(ambit::Tensor::build(tensor_type_, "Bmq", {nthree_, nvirtual_})); + + // Vefu.push_back( + // ambit::Tensor::build(tensor_type_, "muJK", {nvirtual_, nvirtual_, nactive_})); + // Tefv.push_back(ambit::Tensor::build(tensor_type_, "T2", {nvirtual_, nvirtual_, + // nactive_})); + + // tempTaa.push_back(ambit::Tensor::build(tensor_type_, "TEMPaa", {nactive_, + // nactive_})); tempTAA.push_back(ambit::Tensor::build(tensor_type_, "TEMPAA", + // {nactive_, nactive_})); + // } + // // ambit::Tensor BemQ = ints_->three_integral_block(naux, acore_mos_, + // // avirt_mos_); + // // ambit::Tensor BeuQ = ints_->three_integral_block(naux, aactv_mos_, + // // avirt_mos_); + + // // Loop over e and f to compute V + + // ambit::Tensor BeuQ = ints_->three_integral_block(aux_mos_, virt_mos_, actv_mos_); + + // // std::vector& BemQ_data = BemQ.data(); + + // // I think this loop is typically too small to allow efficient use of + // // OpenMP. Should probably test this assumption. + // #pragma omp parallel for num_threads(num_threads_) + // for (size_t m = 0; m < ncore_; m++) { + // int thread = 0; + // #ifdef _OPENMP + // thread = omp_get_thread_num(); + // #endif + // size_t ma = core_mos_[m]; + + // // V[efu]_m = B_{em}^Q * B_{fu}^Q - B_{eu}^Q B_{fm}^Q + // // V[efu]_m = V[efmu] + V[efmu] * exp[efmu] + // // T2["mvef"] = V["mvef"] * D["mvef"] + // // temp["uv"] = V * T2 + // #pragma omp critical + // { Bm_Qe[thread] = ints_->three_integral_block_two_index(aux_mos_, ma, virt_mos_); } + + // Vefu[thread]("e, f, u") = Bm_Qe[thread]("Q, e") * BeuQ("Q, f, u"); + // Vefu[thread]("e, f, u") -= BeuQ("Q, e, u") * Bm_Qe[thread]("Q, f"); + + // // E = V["efmu"] (1 + Exp(-s * D^{ef}_{mu}) * V^{mv}_{ef} * + // // Denom^{mv}_{ef} + // Tefv[thread].data() = Vefu[thread].data(); + + // std::vector& T_mv_data = Tefv[thread].data(); + // Vefu[thread].iterate([&](const std::vector& i, double& value) { + // double Exp = + // Fa_[virt_mos_[i[0]]] + Fa_[virt_mos_[i[1]]] - Fa_[actv_mos_[i[2]]] - Fa_[ma]; + // double D = -1.0 * (Fa_[virt_mos_[i[0]]] + Fa_[virt_mos_[i[1]]] - + // Fa_[actv_mos_[i[2]]] - + // Fa_[ma]); + // value = value + value * dsrg_source_->compute_renormalized(Exp); + // T_mv_data[i[0] * nvirtual_ * nactive_ + i[1] * nactive_ + i[2]] *= + // dsrg_source_->compute_renormalized_denominator(D); + // }); + + // // T_mv[thread].iterate([&](const std::vector& i,double& + // // value){ + // // double D = Fa_[aactv_mos_[i[1]]] + Fa_[acore_mos_[i[0]]] - + // // Fa_[ea] - Fa_[fa]; + // // value = value * + // // dsrg_source_->compute_renormalized_denominator(D);}); + + // tempTaa[thread]("u,v") += 0.5 * Vefu[thread]("e, f, u") * Tefv[thread]("e, f, v"); + // Vefu[thread].zero(); + // Tefv[thread].zero(); + + // Vefu[thread].zero(); + // Vefu[thread]("e, f, u") = Bm_Qe[thread]("Q, e") * BeuQ("Q, f, u"); + + // // E = V["efmu"] (1 + Exp(-s * D^{ef}_{mu}) * V^{mv}_{ef} * + // // Denom^{mv}_{ef} + // Tefv[thread].data() = Vefu[thread].data(); + // T_mv_data = Tefv[thread].data(); + // T_mv_data = Tefv[thread].data(); + // Vefu[thread].iterate([&](const std::vector& i, double& value) { + // double Exp = + // Fa_[virt_mos_[i[0]]] + Fb_[virt_mos_[i[1]]] - Fa_[actv_mos_[i[2]]] - Fb_[ma]; + // double D = -1.0 * (Fa_[virt_mos_[i[0]]] + Fb_[virt_mos_[i[1]]] - + // Fa_[actv_mos_[i[2]]] - + // Fb_[ma]); + // value = value + value * dsrg_source_->compute_renormalized(Exp); + // T_mv_data[i[0] * nvirtual_ * nactive_ + i[1] * nactive_ + i[2]] *= + // dsrg_source_->compute_renormalized_denominator(D); + // }); + + // // T_mv[thread].iterate([&](const std::vector& i,double& + // // value){ + // // double D = Fa_[aactv_mos_[i[1]]] + Fa_[acore_mos_[i[0]]] - + // // Fa_[ea] - Fa_[fa]; + // // value = value * + // // dsrg_source_->compute_renormalized_denominator(D);}); + + // tempTAA[thread]("vu") += Vefu[thread]("e, f, u") * Tefv[thread]("e,f, v"); + // tempTaa[thread]("vu") += Vefu[thread]("e,f, u") * Tefv[thread]("e,f, v"); + // Vefu[thread].zero(); + // Tefv[thread].zero(); + + // Vefu[thread]("e, f, u") = Bm_Qe[thread]("Q, e") * BeuQ("Q, f, u"); + // Vefu[thread]("e, f, u") -= BeuQ("Q, e, u") * Bm_Qe[thread]("Q, f"); + + // // E = V["efmu"] (1 + Exp(-s * D^{ef}_{mu}) * V^{mv}_{ef} * + // // Denom^{mv}_{ef} + + // Tefv[thread].data() = Vefu[thread].data(); + // T_mv_data = Tefv[thread].data(); + + // T_mv_data = Tefv[thread].data(); + // Vefu[thread].iterate([&](const std::vector& i, double& value) { + // double Exp = + // Fa_[virt_mos_[i[0]]] + Fb_[virt_mos_[i[1]]] - Fb_[actv_mos_[i[2]]] - Fb_[ma]; + // double D = -1.0 * (Fa_[virt_mos_[i[0]]] + Fa_[virt_mos_[i[1]]] - + // Fb_[actv_mos_[i[2]]] - + // Fb_[ma]); + // value = value + value * dsrg_source_->compute_renormalized(Exp); + // T_mv_data[i[0] * nvirtual_ * nactive_ + i[1] * nactive_ + i[2]] *= + // dsrg_source_->compute_renormalized_denominator(D); + // }); + + // tempTaa[thread]("u,v") += 0.5 * Vefu[thread]("e, f, u") * Tefv[thread]("e, f, v"); + // } + + // ambit::Tensor tempTAA_all = + // ambit::Tensor::build(tensor_type_, "tempTAA_all", {nactive_, nactive_}); + // ambit::Tensor tempTaa_all = + // ambit::Tensor::build(tensor_type_, "tempTaa_all", {nactive_, nactive_}); + // for (int thread = 0; thread < nthread; thread++) { + // tempTAA_all("v, u") += tempTAA[thread]("v, u"); + // tempTaa_all("v, u") += tempTaa[thread]("v, u"); + // } + + // Eacvv += tempTAA_all("v,u") * Gamma1_AA("v,u"); + // Eacvv += tempTaa_all("v,u") * Gamma1_aa("v,u"); + + // if (print_ > 0) { + // outfile->Printf("\n\n CAVV computation takes %8.8f", cavvTimer.get()); + // } + + if (laplace_ccav) { + outfile->Printf("\n CCAV has not been calculated. Laplace algorithm will be used."); + } else { + std::vector Bm_vQ; + std::vector Bn_eQ; + std::vector Bm_eQ; + std::vector Bn_vQ; + + std::vector V_eu; + std::vector T_ev; + std::vector tempTaa_e; + std::vector tempTAA_e; + + ambit::Tensor BmvQ = ints_->three_integral_block(aux_mos_, core_mos_, actv_mos_); + ambit::Tensor BmvQ_swapped = + ambit::Tensor::build(tensor_type_, "Bm_vQ", {ncore_, nthree_, nactive_}); + BmvQ_swapped("m, Q, u") = BmvQ("Q, m, u"); + local_timer ccvaTimer; + for (int thread = 0; thread < nthread; thread++) { + Bm_vQ.push_back(ambit::Tensor::build(tensor_type_, "BemQ", {nthree_, nactive_})); + Bn_eQ.push_back(ambit::Tensor::build(tensor_type_, "Bf_uQ", {nthree_, nvirtual_})); + Bm_eQ.push_back(ambit::Tensor::build(tensor_type_, "Bmq", {nthree_, nvirtual_})); + Bn_vQ.push_back(ambit::Tensor::build(tensor_type_, "Bmq", {nthree_, nactive_})); + + V_eu.push_back(ambit::Tensor::build(tensor_type_, "muJK", {nvirtual_, nactive_})); + T_ev.push_back(ambit::Tensor::build(tensor_type_, "T2", {nvirtual_, nactive_})); + + tempTaa_e.push_back(ambit::Tensor::build(tensor_type_, "TEMPaa", {nactive_, nactive_})); + tempTAA_e.push_back(ambit::Tensor::build(tensor_type_, "TEMPAA", {nactive_, nactive_})); + } + ambit::Tensor Eta1_aa = Eta1_.block("aa"); + ambit::Tensor Eta1_AA = Eta1_.block("AA"); #pragma omp parallel for num_threads(num_threads_) - for (size_t m = 0; m < ncore_; ++m) { - size_t ma = core_mos_[m]; - size_t mb = core_mos_[m]; - int thread = 0; + for (size_t m = 0; m < ncore_; ++m) { + size_t ma = core_mos_[m]; + size_t mb = core_mos_[m]; + int thread = 0; #ifdef _OPENMP - thread = omp_get_thread_num(); + thread = omp_get_thread_num(); #endif #pragma omp critical - { Bm_eQ[thread] = ints_->three_integral_block_two_index(aux_mos_, ma, virt_mos_); } - std::copy(&BmvQ_swapped.data()[m * nthree_ * nactive_], - &BmvQ_swapped.data()[m * nthree_ * nactive_ + nthree_ * nactive_], - Bm_vQ[thread].data().begin()); - - for (size_t n = 0; n < ncore_; ++n) { - // alpha-aplha - size_t na = core_mos_[n]; - size_t nb = core_mos_[n]; - - std::copy(&BmvQ_swapped.data()[n * nthree_ * nactive_], - &BmvQ_swapped.data()[n * nthree_ * nactive_ + nthree_ * nactive_], - Bn_vQ[thread].data().begin()); -// Bn_vQ[thread].iterate([&](const std::vector& i,double& value){ -// value = BmvQ_data[i[0] * core_ * active_ + n * active_ + i[1] ]; -// }); -#pragma omp critical - { Bn_eQ[thread] = ints_->three_integral_block_two_index(aux_mos_, na, virt_mos_); } + { Bm_eQ[thread] = ints_->three_integral_block_two_index(aux_mos_, ma, virt_mos_); } + std::copy(&BmvQ_swapped.data()[m * nthree_ * nactive_], + &BmvQ_swapped.data()[m * nthree_ * nactive_ + nthree_ * nactive_], + Bm_vQ[thread].data().begin()); - // B_{mv}^{Q} * B_{ne}^{Q} - B_{me}^Q * B_{nv} - V_eu[thread]("e, u") = Bm_vQ[thread]("Q, u") * Bn_eQ[thread]("Q, e"); - V_eu[thread]("e, u") -= Bm_eQ[thread]("Q, e") * Bn_vQ[thread]("Q, u"); - // E = V["efmu"] (1 + Exp(-s * D^{ef}_{mu}) * V^{mv}_{ef} * - // Denom^{mv}_{ef} - T_ev[thread].data() = V_eu[thread].data(); + for (size_t n = 0; n < ncore_; ++n) { + // alpha-aplha + size_t na = core_mos_[n]; + size_t nb = core_mos_[n]; - V_eu[thread].iterate([&](const std::vector& i, double& value) { - double Exp = Fa_[actv_mos_[i[1]]] + Fa_[virt_mos_[i[0]]] - Fa_[ma] - Fa_[na]; - value = value + value * dsrg_source_->compute_renormalized(Exp); - double D = Fa_[ma] + Fa_[na] - Fa_[actv_mos_[i[1]]] - Fa_[virt_mos_[i[0]]]; - T_ev[thread].data()[i[0] * nactive_ + i[1]] *= - dsrg_source_->compute_renormalized_denominator(D); - ; - }); + std::copy(&BmvQ_swapped.data()[n * nthree_ * nactive_], + &BmvQ_swapped.data()[n * nthree_ * nactive_ + nthree_ * nactive_], + Bn_vQ[thread].data().begin()); + // Bn_vQ[thread].iterate([&](const std::vector& i,double& value){ + // value = BmvQ_data[i[0] * core_ * active_ + n * active_ + i[1] ]; + // }); +#pragma omp critical + { Bn_eQ[thread] = ints_->three_integral_block_two_index(aux_mos_, na, virt_mos_); } + + // B_{mv}^{Q} * B_{ne}^{Q} - B_{me}^Q * B_{nv} + V_eu[thread]("e, u") = Bm_vQ[thread]("Q, u") * Bn_eQ[thread]("Q, e"); + V_eu[thread]("e, u") -= Bm_eQ[thread]("Q, e") * Bn_vQ[thread]("Q, u"); + // E = V["efmu"] (1 + Exp(-s * D^{ef}_{mu}) * V^{mv}_{ef} * + // Denom^{mv}_{ef} + T_ev[thread].data() = V_eu[thread].data(); + + V_eu[thread].iterate([&](const std::vector& i, double& value) { + double Exp = Fa_[actv_mos_[i[1]]] + Fa_[virt_mos_[i[0]]] - Fa_[ma] - Fa_[na]; + value = value + value * dsrg_source_->compute_renormalized(Exp); + double D = Fa_[ma] + Fa_[na] - Fa_[actv_mos_[i[1]]] - Fa_[virt_mos_[i[0]]]; + T_ev[thread].data()[i[0] * nactive_ + i[1]] *= + dsrg_source_->compute_renormalized_denominator(D); + ; + }); - tempTaa_e[thread]("u,v") += 0.5 * V_eu[thread]("e,u") * T_ev[thread]("e,v"); - V_eu[thread].zero(); - T_ev[thread].zero(); - - // alpha-beta - // temp["vu"] += V_["vEmN"] * T2_["mNuE"]; - // - V_eu[thread]("E,u") = Bm_vQ[thread]("Q, u") * Bn_eQ[thread]("Q, E"); - T_ev[thread].data() = V_eu[thread].data(); - V_eu[thread].iterate([&](const std::vector& i, double& value) { - double Exp = Fa_[actv_mos_[i[1]]] + Fb_[virt_mos_[i[0]]] - Fa_[ma] - Fb_[nb]; - value = value + value * dsrg_source_->compute_renormalized(Exp); - double D = Fa_[ma] + Fb_[nb] - Fa_[actv_mos_[i[1]]] - Fb_[virt_mos_[i[0]]]; - T_ev[thread].data()[i[0] * nactive_ + i[1]] *= - dsrg_source_->compute_renormalized_denominator(D); - ; - }); + tempTaa_e[thread]("u,v") += 0.5 * V_eu[thread]("e,u") * T_ev[thread]("e,v"); + V_eu[thread].zero(); + T_ev[thread].zero(); - tempTAA_e[thread]("vu") += V_eu[thread]("M,v") * T_ev[thread]("M,u"); - tempTaa_e[thread]("vu") += V_eu[thread]("M,v") * T_ev[thread]("M, u"); + // alpha-beta + // temp["vu"] += V_["vEmN"] * T2_["mNuE"]; + // + V_eu[thread]("E,u") = Bm_vQ[thread]("Q, u") * Bn_eQ[thread]("Q, E"); + T_ev[thread].data() = V_eu[thread].data(); + V_eu[thread].iterate([&](const std::vector& i, double& value) { + double Exp = Fa_[actv_mos_[i[1]]] + Fb_[virt_mos_[i[0]]] - Fa_[ma] - Fb_[nb]; + value = value + value * dsrg_source_->compute_renormalized(Exp); + double D = Fa_[ma] + Fb_[nb] - Fa_[actv_mos_[i[1]]] - Fb_[virt_mos_[i[0]]]; + T_ev[thread].data()[i[0] * nactive_ + i[1]] *= + dsrg_source_->compute_renormalized_denominator(D); + ; + }); - // beta-beta - V_eu[thread].zero(); - T_ev[thread].zero(); - V_eu[thread]("E,U") = Bm_vQ[thread]("Q, U") * Bn_eQ[thread]("Q,E"); - V_eu[thread]("E,U") -= Bm_eQ[thread]("Q, E") * Bn_vQ[thread]("Q, U"); - T_ev[thread].data() = V_eu[thread].data(); + tempTAA_e[thread]("vu") += V_eu[thread]("M,v") * T_ev[thread]("M,u"); + tempTaa_e[thread]("vu") += V_eu[thread]("M,v") * T_ev[thread]("M, u"); - V_eu[thread].iterate([&](const std::vector& i, double& value) { - double Exp = Fb_[mb] + Fb_[nb] - Fb_[actv_mos_[i[1]]] - Fb_[virt_mos_[i[0]]]; - value = value + value * dsrg_source_->compute_renormalized(Exp); - double D = Fb_[mb] + Fb_[nb] - Fb_[actv_mos_[i[1]]] - Fb_[virt_mos_[i[0]]]; - T_ev[thread].data()[i[0] * nactive_ + i[1]] *= - dsrg_source_->compute_renormalized_denominator(D); - ; - }); + // beta-beta + V_eu[thread].zero(); + T_ev[thread].zero(); + V_eu[thread]("E,U") = Bm_vQ[thread]("Q, U") * Bn_eQ[thread]("Q,E"); + V_eu[thread]("E,U") -= Bm_eQ[thread]("Q, E") * Bn_vQ[thread]("Q, U"); + T_ev[thread].data() = V_eu[thread].data(); + + V_eu[thread].iterate([&](const std::vector& i, double& value) { + double Exp = Fb_[mb] + Fb_[nb] - Fb_[actv_mos_[i[1]]] - Fb_[virt_mos_[i[0]]]; + value = value + value * dsrg_source_->compute_renormalized(Exp); + double D = Fb_[mb] + Fb_[nb] - Fb_[actv_mos_[i[1]]] - Fb_[virt_mos_[i[0]]]; + T_ev[thread].data()[i[0] * nactive_ + i[1]] *= + dsrg_source_->compute_renormalized_denominator(D); + ; + }); - tempTAA_e[thread]("v,u") += 0.5 * V_eu[thread]("M,v") * T_ev[thread]("M,u"); - V_eu[thread].zero(); - T_ev[thread].zero(); + tempTAA_e[thread]("v,u") += 0.5 * V_eu[thread]("M,v") * T_ev[thread]("M,u"); + V_eu[thread].zero(); + T_ev[thread].zero(); + } } - } - tempTAA_all = ambit::Tensor::build(tensor_type_, "tempTAA_all", {nactive_, nactive_}); - tempTaa_all = ambit::Tensor::build(tensor_type_, "tempTaa_all", {nactive_, nactive_}); - for (int thread = 0; thread < nthread; thread++) { - tempTAA_all("u, v") += tempTAA_e[thread]("u,v"); - tempTaa_all("u, v") += tempTaa_e[thread]("u,v"); - } - Eccva += tempTaa_all("vu") * Eta1_aa("uv"); - Eccva += tempTAA_all("VU") * Eta1_AA("UV"); - if (print_ > 1) { - outfile->Printf("\n\n CCVA takes %8.8f", cavvTimer.get()); + ambit::Tensor tempTAA_all = + ambit::Tensor::build(tensor_type_, "tempTAA_all", {nactive_, nactive_}); + ambit::Tensor tempTaa_all = + ambit::Tensor::build(tensor_type_, "tempTaa_all", {nactive_, nactive_}); + for (int thread = 0; thread < nthread; thread++) { + tempTAA_all("u, v") += tempTAA_e[thread]("u,v"); + tempTaa_all("u, v") += tempTaa_e[thread]("u,v"); + } + Eccva += tempTaa_all("vu") * Eta1_aa("uv"); + Eccva += tempTAA_all("VU") * Eta1_AA("UV"); + if (print_ > 0) { + outfile->Printf("\n\n CCVA takes %8.8f", ccvaTimer.get()); + } } + // std::vector Bm_vQ; + // std::vector Bn_eQ; + // std::vector Bm_eQ; + // std::vector Bn_vQ; + + // std::vector V_eu; + // std::vector T_ev; + // std::vector tempTaa_e; + // std::vector tempTAA_e; + + // ambit::Tensor BmvQ = ints_->three_integral_block(aux_mos_, core_mos_, actv_mos_); + // ambit::Tensor BmvQ_swapped = + // ambit::Tensor::build(tensor_type_, "Bm_vQ", {ncore_, nthree_, nactive_}); + // BmvQ_swapped("m, Q, u") = BmvQ("Q, m, u"); + // local_timer ccvaTimer; + // for (int thread = 0; thread < nthread; thread++) { + // Bm_vQ.push_back(ambit::Tensor::build(tensor_type_, "BemQ", {nthree_, nactive_})); + // Bn_eQ.push_back(ambit::Tensor::build(tensor_type_, "Bf_uQ", {nthree_, nvirtual_})); + // Bm_eQ.push_back(ambit::Tensor::build(tensor_type_, "Bmq", {nthree_, nvirtual_})); + // Bn_vQ.push_back(ambit::Tensor::build(tensor_type_, "Bmq", {nthree_, nactive_})); + + // V_eu.push_back(ambit::Tensor::build(tensor_type_, "muJK", {nvirtual_, nactive_})); + // T_ev.push_back(ambit::Tensor::build(tensor_type_, "T2", {nvirtual_, nactive_})); + + // tempTaa_e.push_back(ambit::Tensor::build(tensor_type_, "TEMPaa", {nactive_, + // nactive_})); tempTAA_e.push_back(ambit::Tensor::build(tensor_type_, "TEMPAA", + // {nactive_, nactive_})); + // } + // ambit::Tensor Eta1_aa = Eta1_.block("aa"); + // ambit::Tensor Eta1_AA = Eta1_.block("AA"); + + // #pragma omp parallel for num_threads(num_threads_) + // for (size_t m = 0; m < ncore_; ++m) { + // size_t ma = core_mos_[m]; + // size_t mb = core_mos_[m]; + // int thread = 0; + // #ifdef _OPENMP + // thread = omp_get_thread_num(); + // #endif + + // #pragma omp critical + // { Bm_eQ[thread] = ints_->three_integral_block_two_index(aux_mos_, ma, virt_mos_); } + // std::copy(&BmvQ_swapped.data()[m * nthree_ * nactive_], + // &BmvQ_swapped.data()[m * nthree_ * nactive_ + nthree_ * nactive_], + // Bm_vQ[thread].data().begin()); + + // for (size_t n = 0; n < ncore_; ++n) { + // // alpha-aplha + // size_t na = core_mos_[n]; + // size_t nb = core_mos_[n]; + + // std::copy(&BmvQ_swapped.data()[n * nthree_ * nactive_], + // &BmvQ_swapped.data()[n * nthree_ * nactive_ + nthree_ * nactive_], + // Bn_vQ[thread].data().begin()); + // // Bn_vQ[thread].iterate([&](const std::vector& i,double& value){ + // // value = BmvQ_data[i[0] * core_ * active_ + n * active_ + i[1] ]; + // // }); + // #pragma omp critical + // { Bn_eQ[thread] = ints_->three_integral_block_two_index(aux_mos_, na, virt_mos_); + // } + + // // B_{mv}^{Q} * B_{ne}^{Q} - B_{me}^Q * B_{nv} + // V_eu[thread]("e, u") = Bm_vQ[thread]("Q, u") * Bn_eQ[thread]("Q, e"); + // V_eu[thread]("e, u") -= Bm_eQ[thread]("Q, e") * Bn_vQ[thread]("Q, u"); + // // E = V["efmu"] (1 + Exp(-s * D^{ef}_{mu}) * V^{mv}_{ef} * + // // Denom^{mv}_{ef} + // T_ev[thread].data() = V_eu[thread].data(); + + // V_eu[thread].iterate([&](const std::vector& i, double& value) { + // double Exp = Fa_[actv_mos_[i[1]]] + Fa_[virt_mos_[i[0]]] - Fa_[ma] - Fa_[na]; + // value = value + value * dsrg_source_->compute_renormalized(Exp); + // double D = Fa_[ma] + Fa_[na] - Fa_[actv_mos_[i[1]]] - Fa_[virt_mos_[i[0]]]; + // T_ev[thread].data()[i[0] * nactive_ + i[1]] *= + // dsrg_source_->compute_renormalized_denominator(D); + // ; + // }); + + // tempTaa_e[thread]("u,v") += 0.5 * V_eu[thread]("e,u") * T_ev[thread]("e,v"); + // V_eu[thread].zero(); + // T_ev[thread].zero(); + + // // alpha-beta + // // temp["vu"] += V_["vEmN"] * T2_["mNuE"]; + // // + // V_eu[thread]("E,u") = Bm_vQ[thread]("Q, u") * Bn_eQ[thread]("Q, E"); + // T_ev[thread].data() = V_eu[thread].data(); + // V_eu[thread].iterate([&](const std::vector& i, double& value) { + // double Exp = Fa_[actv_mos_[i[1]]] + Fb_[virt_mos_[i[0]]] - Fa_[ma] - Fb_[nb]; + // value = value + value * dsrg_source_->compute_renormalized(Exp); + // double D = Fa_[ma] + Fb_[nb] - Fa_[actv_mos_[i[1]]] - Fb_[virt_mos_[i[0]]]; + // T_ev[thread].data()[i[0] * nactive_ + i[1]] *= + // dsrg_source_->compute_renormalized_denominator(D); + // ; + // }); + + // tempTAA_e[thread]("vu") += V_eu[thread]("M,v") * T_ev[thread]("M,u"); + // tempTaa_e[thread]("vu") += V_eu[thread]("M,v") * T_ev[thread]("M, u"); + + // // beta-beta + // V_eu[thread].zero(); + // T_ev[thread].zero(); + // V_eu[thread]("E,U") = Bm_vQ[thread]("Q, U") * Bn_eQ[thread]("Q,E"); + // V_eu[thread]("E,U") -= Bm_eQ[thread]("Q, E") * Bn_vQ[thread]("Q, U"); + // T_ev[thread].data() = V_eu[thread].data(); + + // V_eu[thread].iterate([&](const std::vector& i, double& value) { + // double Exp = Fb_[mb] + Fb_[nb] - Fb_[actv_mos_[i[1]]] - Fb_[virt_mos_[i[0]]]; + // value = value + value * dsrg_source_->compute_renormalized(Exp); + // double D = Fb_[mb] + Fb_[nb] - Fb_[actv_mos_[i[1]]] - Fb_[virt_mos_[i[0]]]; + // T_ev[thread].data()[i[0] * nactive_ + i[1]] *= + // dsrg_source_->compute_renormalized_denominator(D); + // ; + // }); + + // tempTAA_e[thread]("v,u") += 0.5 * V_eu[thread]("M,v") * T_ev[thread]("M,u"); + // V_eu[thread].zero(); + // T_ev[thread].zero(); + // } + // } + + // tempTAA_all = ambit::Tensor::build(tensor_type_, "tempTAA_all", {nactive_, nactive_}); + // tempTaa_all = ambit::Tensor::build(tensor_type_, "tempTaa_all", {nactive_, nactive_}); + // for (int thread = 0; thread < nthread; thread++) { + // tempTAA_all("u, v") += tempTAA_e[thread]("u,v"); + // tempTaa_all("u, v") += tempTaa_e[thread]("u,v"); + // } + // Eccva += tempTaa_all("vu") * Eta1_aa("uv"); + // Eccva += tempTAA_all("VU") * Eta1_AA("UV"); + // if (print_ > 0) { + // outfile->Printf("\n\n CCVA takes %8.8f", ccvaTimer.get()); + // } + + outfile->Printf("\n E_cavv: %8.6f", Eacvv); + outfile->Printf("\n E_ccav: %8.6f", Eccva); return (Eacvv + Eccva); } diff --git a/forte/mrdsrg-spin-integrated/three_dsrg_mrpt2.h b/forte/mrdsrg-spin-integrated/three_dsrg_mrpt2.h index ff17084ac..371a79a92 100644 --- a/forte/mrdsrg-spin-integrated/three_dsrg_mrpt2.h +++ b/forte/mrdsrg-spin-integrated/three_dsrg_mrpt2.h @@ -209,12 +209,19 @@ class THREE_DSRG_MRPT2 : public MASTER_DSRG { double E_VT2_2_batch_virtual_mpi(); double E_VT2_2_batch_virtual_ga(); double E_VT2_2_batch_virtual_rep(); - double E_VT2_2_AO_Slow(); + //double E_VT2_2_AO_Slow(); double E_VT2_4PP(); double E_VT2_4HH(); double E_VT2_4PH(); double E_VT2_6(); + // double E_cavv_lt(); + // double E_ccav_lt(); + + double E_ccvv_df_ao(); + double E_ccvv_lt_ao(); + double E_ccvv_diskdf_ao(); + void de_normal_order(); /// Form Hbar for reference relaxation diff --git a/forte/orbital-helpers/Laplace.cc b/forte/orbital-helpers/Laplace.cc new file mode 100644 index 000000000..666cfed59 --- /dev/null +++ b/forte/orbital-helpers/Laplace.cc @@ -0,0 +1,294 @@ +#include "psi4/libpsi4util/PsiOutStream.h" +#include "psi4/libpsio/psio.h" +#include "psi4/libpsio/psio.hpp" +#include "psi4/libpsi4util/process.h" +#include "psi4/libmints/basisset.h" +#include "psi4/libmints/wavefunction.h" +#include "psi4/libmints/integral.h" +#include "psi4/libmints/matrix.h" +#include "psi4/libmints/sieve.h" +#include "psi4/lib3index/cholesky.h" +#include "psi4/libqt/qt.h" +#include "psi4/psifiles.h" +#include "psi4/lib3index/dftensor.h" + +#include "helpers/timer.h" +#include "helpers/printing.h" +#include "helpers/memory.h" + +#include "base_classes/forte_options.h" + +#include "Laplace.h" + +#include +#include + + +using namespace psi; + +namespace forte { + +std::vector merge_lists(const std::vector &l1, const std::vector &l2) { + + std::vector l12; + + int i1 = 0, i2 = 0; + while(i1 < l1.size() || i2 < l2.size()) { + if(i1 == l1.size()) { + l12.push_back(l2[i2]); + ++i2; + } else if(i2 == l2.size()) { + l12.push_back(l1[i1]); + ++i1; + } else if(l1[i1] == l2[i2]) { + l12.push_back(l1[i1]); + ++i1; + ++i2; + } else if(l1[i1] < l2[i2]) { + l12.push_back(l1[i1]); + ++i1; + } else { + l12.push_back(l2[i2]); + ++i2; + } + } + + return l12; + +} + +std::vector contract_lists(const std::vector &y, const std::vector> &A_to_y) { + + // TODO: runtime is proportional to A_to_y size (system size, O(N)) + // could maybe reduce to &y size (domain size, O(1)), probably doesn't matter + std::vector yA; + + for(int a = 0, y_ind = 0; a < A_to_y.size(); ++a) { + + bool is_a = false; + for(auto y_val : A_to_y[a]) { + if (y_ind < y.size() && y[y_ind] == y_val) { + y_ind++; + is_a = true; + } + } + + if(is_a) { + for(auto y_val : A_to_y[a]) { + yA.push_back(y_val); + } + } + + } + + return yA; + +} + +std::vector block_list(const std::vector &x_list, const std::vector &x_to_y_map) { + + std::vector y_list; + + for(int x_val : x_list) { + int y_val = x_to_y_map[x_val]; + if(y_list.size() == 0) { + y_list.push_back(y_val); + } else if(y_list[y_list.size() - 1] != y_val) { + y_list.push_back(y_val); + } + } + + return y_list; + +} + +std::vector> invert_map(const std::vector> &x_to_y, int ny) { + + int nx = x_to_y.size(); + std::vector> y_to_x(ny); + + for(int x = 0; x < nx; x++) { + for(auto y : x_to_y[x]) { + y_to_x[y].push_back(x); + } + } + + return y_to_x; + +} + +std::vector> chain_maps(const std::vector> &x_to_y, const std::vector> &y_to_z) { + + int nx = x_to_y.size(); + std::vector> x_to_z(nx); + + for(int x = 0; x < nx; x++) { + for(auto y : x_to_y[x]) { + for(auto z : y_to_z[y]) { + x_to_z[x].push_back(z); + } + } + std::sort(x_to_z[x].begin(), x_to_z[x].end()); + x_to_z[x].erase(std::unique(x_to_z[x].begin(), x_to_z[x].end()), x_to_z[x].end()); + } + + return x_to_z; + +} + +std::vector> extend_maps(const std::vector> &x_to_y, const std::vector> &xpairs) { + + int nx = x_to_y.size(); + std::vector> xext_to_y(nx); + + for(auto xpair : xpairs) { + size_t x1, x2; + std::tie(x1,x2) = xpair; + xext_to_y[x1] = merge_lists(xext_to_y[x1], x_to_y[x2]); + } + + return xext_to_y; + +} + +psi::SharedMatrix submatrix_rows(const psi::Matrix &mat, const std::vector &row_inds) { + + psi::SharedMatrix mat_new = std::make_shared(mat.name(), row_inds.size(), mat.colspi(0)); + double** mat_newp = mat_new->pointer(); + double** matp = mat.pointer(); + for(int r_new = 0; r_new < row_inds.size(); r_new++) { + int r_old = row_inds[r_new]; + ::memcpy(&mat_newp[r_new][0], &matp[r_old][0], sizeof(double) * mat.colspi(0)); + } + return mat_new; +} + +psi::SharedMatrix submatrix_cols(const psi::Matrix &mat, const std::vector &col_inds) { + + psi::SharedMatrix mat_new = std::make_shared(mat.name(), mat.rowspi(0), col_inds.size()); + double** mat_newp = mat_new->pointer(); + double** matp = mat.pointer(); + for(int c_new = 0; c_new < col_inds.size(); c_new++) { + int c_old = col_inds[c_new]; + C_DCOPY(mat.rowspi(0), &matp[0][c_old], mat.colspi(0), &mat_newp[0][c_new], col_inds.size()); + } + return mat_new; +} + +psi::SharedMatrix submatrix_rows_and_cols(const psi::Matrix &mat, const std::vector &row_inds, const std::vector &col_inds) { + + psi::SharedMatrix mat_new = std::make_shared(mat.name(), row_inds.size(), col_inds.size()); + for(int r_new = 0; r_new < row_inds.size(); r_new++) { + int r_old = row_inds[r_new]; + for(int c_new = 0; c_new < col_inds.size(); c_new++) { + int c_old = col_inds[c_new]; + mat_new->set(r_new, c_new, mat.get(r_old, c_old)); + } + } + return mat_new; +} + +psi::SharedMatrix load_Amn(const size_t A, const size_t mn) { + std::shared_ptr psio(new PSIO()); + auto Amn = std::make_shared("Load (A|mn) from DF-SCF", A, mn); + outfile->Printf("\t Will attempt to load (A|mn) from file %d.\n\n", PSIF_DFSCF_BJ); + double** Amnp = Amn->pointer(); + int file_unit = PSIF_DFSCF_BJ; + psio->open(file_unit, PSIO_OPEN_OLD); + psio->read_entry(file_unit, "ERFC Integrals", (char*)Amnp[0], sizeof(double) * A * mn); + psio->close(file_unit, 1); + return Amn; +} + +psi::SharedMatrix load_Jinv_full(const size_t P, const size_t Q) { + std::shared_ptr psio(new PSIO()); + auto Jinv_full = std::make_shared("Load full inverse metric from DF-SCF", P, Q); + outfile->Printf("\t Will attempt to load full inverse metric from file %d.\n\n", PSIF_DFSCF_BJ); + double** Jinv_fullp = Jinv_full->pointer(); + int file_unit = PSIF_DFSCF_BJ; + psio->open(file_unit, PSIO_OPEN_OLD); + psio->read_entry(file_unit, "Jinv_full", (char*)Jinv_fullp[0], sizeof(double) * P * Q); + psio->close(file_unit, 1); + return Jinv_full; +} + +psi::SharedMatrix initialize_erfc_integral(double Omega, int n_func_pairs, std::shared_ptr ints_forte) { + std::shared_ptr zero = psi::BasisSet::zero_ao_basis_set(); + std::shared_ptr primary = ints_forte->wfn()->basisset(); + std::shared_ptr auxiliary = ints_forte->wfn()->get_basisset("DF_BASIS_MP2"); + + std::shared_ptr integral = std::make_shared(auxiliary,zero,primary,primary); + std::shared_ptr ints(integral->erf_complement_eri(Omega)); + ///std::shared_ptr ints(integral->eri()); + + int nthree = auxiliary->nbf(); + int nbf = primary->nbf(); + + auto I = std::make_shared("erfc integral", nthree, n_func_pairs); + double **Ip = I->pointer(); + + int numP, Pshell, MU, NU, P, PHI, mu, nu, nummu, numnu, omu, onu; + + for (MU = 0; MU < primary->nshell(); ++MU) { + nummu = primary->shell(MU).nfunction(); + for (NU = 0; NU <= MU; ++NU) { + numnu = primary->shell(NU).nfunction(); + for (Pshell = 0; Pshell < auxiliary->nshell(); ++Pshell) { + numP = auxiliary->shell(Pshell).nfunction(); + ints->compute_shell(Pshell, 0, MU, NU); + const double *buffer = ints->buffer(); + for (mu = 0; mu < nummu; ++mu) { + omu = primary->shell(MU).function_index() + mu; + for (nu = 0; nu < numnu; ++nu) { + onu = primary->shell(NU).function_index() + nu; + size_t addr = omu > onu ? omu * (omu + 1) / 2 + onu : + onu * (onu + 1) / 2 + omu; + for (P = 0; P < numP; ++P) { + PHI = auxiliary->shell(Pshell).function_index() + P; + Ip[PHI][addr] = buffer[P * nummu * numnu + mu * numnu + nu]; + } + } + } + } + } + } + return I; +} + +psi::SharedMatrix erfc_metric (double Omega, std::shared_ptr ints_forte) { + std::shared_ptr auxiliary = ints_forte->wfn()->get_basisset("DF_BASIS_MP2"); + auto Jinv = std::make_shared(auxiliary, true); + auto Jinv_no_erfc = std::make_shared(auxiliary, true); + Jinv->form_full_eig_inverse_erfc(Omega, 1E-12); + Jinv_no_erfc->form_fitting_metric(); + ///Jinv->form_full_eig_inverse(1E-12); + psi::SharedMatrix Jinv_metric = Jinv->get_metric(); + psi::SharedMatrix Jinv_metric_no_erfc = Jinv_no_erfc->get_metric(); + psi::SharedMatrix Jinv_combine = psi::linalg::triplet(Jinv_metric, Jinv_metric_no_erfc, Jinv_metric, false, false, false); + return Jinv_combine; +} + +int binary_search_recursive(std::vector A, int key, int low, int high) { + if (low > high) { + return -1; + } + + int mid = low + ((high - low) / 2); + if (A[mid] == key) { + return mid; + } else if (key < A[mid]) { + return binary_search_recursive(A, key, low, mid - 1); + } + + return binary_search_recursive(A, key, mid + 1, high); +} + +psi::SharedMatrix ambit_to_matrix(ambit::Tensor t) { + size_t size1 = t.dim(0); + size_t size2 = t.dim(1); + auto M = std::make_shared("M", size1, size2); + t.iterate([&](const std::vector& i, double& value) { M->set(i[0], i[1], value); }); + return M; +} + +} // namespace forte \ No newline at end of file diff --git a/forte/orbital-helpers/Laplace.h b/forte/orbital-helpers/Laplace.h new file mode 100644 index 000000000..9fc9c4b2f --- /dev/null +++ b/forte/orbital-helpers/Laplace.h @@ -0,0 +1,70 @@ +#ifndef __LAPLACE_H__ +#define __LAPLACE_H__ + +#include "psi4/libmints/matrix.h" +#include "mrdsrg-spin-integrated/master_mrdsrg.h" + + +#include + +typedef std::vector> SparseMap; + +namespace psi { +class PSIO; +} // namespace psi + +namespace forte { +/// Written by Shuhang Li + +std::vector merge_lists(const std::vector &l1, const std::vector &l2); + +/* Args: sorted list of y, sparse map from A to another list of y (assume sorted, each possible value of y appears exactly once in entire map) + * Returns: the union of lists in A_to_y where at least one element is in y + */ +std::vector contract_lists(const std::vector &y, const std::vector> &A_to_y); + +/* Args: x is a list of values (sorted), y is a map from values of x to values of y + * Returns: a list of y values + * + * Multiple values in x may map to the same value in y (i.e. x is a list of bf, y is atoms) + */ +std::vector block_list(const std::vector &x_list, const std::vector &x_to_y_map); + +/* Args: SparseMap from x to y, maximum possible y value +* Returns: SparseMap from y to x +*/ +SparseMap invert_map(const SparseMap &x_to_y, int ny); + +/* Args: SparseMap from x to y, SparseMap from y to z +* Returns: SparseMap from x to z +*/ +SparseMap chain_maps(const SparseMap &x_to_y, const SparseMap &y_to_z); + +/* Args: SparseMap from x to y, list of pairs of type x +* Returns: extended SparseMap from x to y +*/ +SparseMap extend_maps(const SparseMap &i_to_y, const std::vector> &ij_to_i_j); + +/* Args: Matrix, list of row indices */ +psi::SharedMatrix submatrix_rows(const psi::Matrix &mat, const std::vector &row_inds); + +/* Args: Matrix, list of column indices */ +psi::SharedMatrix submatrix_cols(const psi::Matrix &mat, const std::vector &col_inds); + +/* Args: Matrix, list of row and column indices */ +psi::SharedMatrix submatrix_rows_and_cols(const psi::Matrix &mat, const std::vector &row_inds, const std::vector &col_inds); + +psi::SharedMatrix load_Amn(const size_t A, const size_t mn); + +psi::SharedMatrix load_Jinv_full(const size_t P, const size_t Q); + +psi::SharedMatrix initialize_erfc_integral(double Omega, int n_func_pairs, std::shared_ptr ints_forte); + +psi::SharedMatrix erfc_metric (double Omega, std::shared_ptr ints_forte); + +int binary_search_recursive(std::vector A, int key, int low, int high); + +psi::SharedMatrix ambit_to_matrix(ambit::Tensor t); +} // namespace forte + +#endif \ No newline at end of file diff --git a/forte/orbital-helpers/LaplaceDenominator.cc b/forte/orbital-helpers/LaplaceDenominator.cc new file mode 100644 index 000000000..1ec8fec47 --- /dev/null +++ b/forte/orbital-helpers/LaplaceDenominator.cc @@ -0,0 +1,689 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "psi4/lib3index/dftensor.h" +#include "psi4/psifiles.h" +#include "psi4/libpsio/psio.h" +#include "psi4/libqt/qt.h" +#include "psi4/libciomr/libciomr.h" +#include "psi4/libmints/matrix.h" +#include "psi4/libmints/vector.h" +#include "psi4/libpsi4util/PsiOutStream.h" +#include "psi4/libpsi4util/process.h" + +#include "LaplaceDenominator.h" + +// MKL Header +#ifdef USING_LAPACK_MKL +#include +#endif + +// OpenMP Header +//_OPENMP is defined by the compiler if it exists +#ifdef _OPENMP +#include +#endif + +using namespace psi; + +namespace forte { + +// the third parameter of from_string() should be +// one of std::hex, std::dec or std::oct +template +bool from_string(T &t, const std::string &s, std::ios_base &(*f)(std::ios_base &)) { + std::istringstream iss(s); + return !(iss >> f >> t).fail(); +} + + +LaplaceDenominator::LaplaceDenominator(std::shared_ptr eps_occ, std::shared_ptr eps_vir, double delta, double vir_tol) + : eps_occ_(eps_occ), eps_vir_(eps_vir), delta_(delta), vir_tol_(vir_tol){ + decompose_ccvv(); +} + +LaplaceDenominator::LaplaceDenominator(std::shared_ptr eps_occ, std::shared_ptr eps_act, std::shared_ptr eps_vir, double delta, bool cavv, double vir_tol) + : eps_occ_(eps_occ), eps_vir_(eps_vir), eps_act_(eps_act), delta_(delta), cavv_(cavv), vir_tol_(vir_tol) { + if (cavv_) { + decompose_cavv(); + } else { + decompose_ccav(); + } +} + +LaplaceDenominator::~LaplaceDenominator() {} + +void LaplaceDenominator::decompose_ccvv() { + outfile->Printf("\n\n ==> FORTE Laplace Denominator <==\n\n"); + int nocc = eps_occ_->dimpi()[0]; + int nvir = eps_vir_->dimpi()[0]; + + vir_start_ = 0; + + double E_LOMO = eps_occ_->get(0, 0); + double E_HOMO = eps_occ_->get(0, nocc - 1); + //double E_LUMO = eps_vir_->get(0, 0); + double E_HUMO = eps_vir_->get(0, nvir - 1); + + for (int i = 0; i < nvir; i++) { + double E_LUMO = eps_vir_->get(0, i); + double E_test = 2.0 * (E_LUMO - E_HOMO); + if (E_test < vir_tol_) { + vir_start_ += 1; + } else { + break; + } + } + + double E_LUMO = eps_vir_->get(0, vir_start_); + + double A = 2.0 * (E_LUMO - E_HOMO); + double B = 2.0 * (E_HUMO - E_LOMO); + double R = B / A; + //double R = B / 1.044679; + + nvir -= vir_start_; + + outfile->Printf(" The smallest denominator: %f. \n", A); + outfile->Printf(" The largest denominator: %f. \n", B); + outfile->Printf(" Number of unselected virtual orbitals: %d. \n", vir_start_); + + // Pick appropriate quadrature file and read contents + std::string PSIDATADIR = Process::environment.get_datadir(); + std::string err_table_filename = PSIDATADIR + "/quadratures/1_x/error.bin"; + std::string R_filename = PSIDATADIR + "/quadratures/1_x/R_avail.bin"; + + std::ifstream err_table_file(err_table_filename.c_str(), std::ios::in | std::ios::binary); + std::ifstream R_avail_file(R_filename.c_str(), std::ios::in | std::ios::binary); + + if (!err_table_file) + throw PSIEXCEPTION( + "LaplaceQuadrature: Cannot locate error property file for quadrature rules (should be " + "PSIDATADIR/quadratures/1_x/error.bin)"); + if (!R_avail_file) + throw PSIEXCEPTION( + "LaplaceQuadrature: Cannot locate R property file for quadrature rules (should be " + "PSIDATADIR/quadratures/1_x/R_avail.bin)"); + + int nk = 53; + int nR = 99; + + // Read in the R available + auto *R_availp = new double[nR]; + R_avail_file.read((char *)R_availp, nR * sizeof(double)); + + auto err_table = std::make_shared("Error Table (nR x nk)", nR, nk); + double **err_tablep = err_table->pointer(); + err_table_file.read((char *)err_tablep[0], static_cast (nR) * nk * sizeof(double)); + + R_avail_file.close(); + err_table_file.close(); + + // for (int r2 = 0; r2 < nR; r2++) + // outfile->Printf( " R[%4d] = %20.14E\n", r2+1, R_availp[r2]); + // err_table->print(); + + int indR; + for (indR = 0; indR < nR; indR++) { + if (R < R_availp[indR]) break; + } + if (indR == nR) { + // TODO: Relax this + throw PSIEXCEPTION( + "Laplace Quadrature requested for (E_HUMO - E_LOMO)/(E_LUMO-E_HOMO) > 7.0 * 10^12, quadratures are not " + "designed for this range."); + } + + double accuracy; + int k, r; + bool found = false; + for (k = 0; k < nk; k++) { + for (r = indR; r < nR; r++) { + double err = err_tablep[r][k]; + if (err != 0.0 && err < delta_) { + accuracy = err; + found = true; + break; + } + } + if (found) break; + } + + if (!found) { + throw PSIEXCEPTION("Laplace Quadrature rule could not be found with specified accuracy for this system"); + } + + nvector_ = k + 1; + + // A bit hacky, but OK + int exponent = (int)floor(log(R_availp[r]) / log(10.0)); + int mantissa = (int)round(R_availp[r] / pow(10.0, exponent)); + if (mantissa == 10) { + exponent++; + mantissa = 1; + } + + std::stringstream st; + st << std::setfill('0'); + st << "1_xk" << std::setw(2) << nvector_; + st << "_" << mantissa; + st << "E" << exponent; + + std::string quadfile = PSIDATADIR + "/quadratures/1_x/" + st.str().c_str(); + + outfile->Printf(" This system has an intrinsic R = (E_HUMO - E_LOMO)/(E_LUMO - E_HOMO) of %7.4E.\n", R); + outfile->Printf(" A %d point minimax quadrature with R of %1.0E will be used for the denominator.\n", nvector_, + R_availp[r]); + outfile->Printf(" The worst-case Chebyshev norm for this quadrature rule is %7.4E.\n", accuracy); + outfile->Printf(" Quadrature rule read from file %s.\n\n", quadfile.c_str()); + + // The quadrature is defined as \omega_v exp(-\alpha_v x) = 1/x + auto *alpha = new double[nvector_]; + auto *omega = new double[nvector_]; + + std::vector lines; + std::string text; + std::ifstream infile(quadfile.c_str()); + if (!infile) throw PSIEXCEPTION("LaplaceDenominator: Unable to open quadrature rule file: " + quadfile); + while (infile.good()) { + std::getline(infile, text); + lines.push_back(text); + } + +#define NUMBER "((?:[-+]?\\d*\\.\\d+(?:[DdEe][-+]?\\d+)?)|(?:[-+]?\\d+\\.\\d*(?:[DdEe][-+]?\\d+)?))" + std::regex numberline("^\\s*(" NUMBER ").*"); + std::smatch what; + + // We'll be rigorous, the files are extremely well defined + int lineno = 0; + for (int index = 0; index < nvector_; index++) { + std::string line = lines[lineno++]; + if (!std::regex_match(line, what, numberline)) + throw PSIEXCEPTION("LaplaceDenominator: Unable to read grid file line: \n" + line); + if (!from_string(omega[index], what[1], std::dec)) + throw PSIEXCEPTION("LaplaceDenominator: Unable to convert grid file line: \n" + line); + } + for (int index = 0; index < nvector_; index++) { + std::string line = lines[lineno++]; + if (!std::regex_match(line, what, numberline)) + throw PSIEXCEPTION("LaplaceDenominator: Unable to read grid file line: \n" + line); + if (!from_string(alpha[index], what[1], std::dec)) + throw PSIEXCEPTION("LaplaceDenominator: Unable to convert grid file line: \n" + line); + } + + // for (int k = 0; k < nvector_; k++) + // printf(" %24.16E, %24.16E\n", omega[k], alpha[k]); + + //Cast weights back to problem size + for (int k = 0; k < nvector_; k++) { + alpha[k] /= A; + omega[k] /= A; + } + + // Fermi-level + double Fermi = (E_LUMO + E_HOMO)/2; + // double Fermi = 0.0; + + denominator_occ_ = std::make_shared("Occupied Laplace Delta Tensor", nvector_, nocc); + denominator_vir_ = std::make_shared("Virtual Laplace Delta Tensor", nvector_, nvir); + // denominator_ = std::make_shared("OV Laplace Delta Tensor", nvector_, nocc * nvir); + + double **dop = denominator_occ_->pointer(); + double **dvp = denominator_vir_->pointer(); + // double **dovp = denominator_->pointer(); + + double *e_o = eps_occ_->pointer(); + double *e_v = eps_vir_->pointer(); + + for (int k = 0; k < nvector_; k++) { + for (int i = 0; i < nocc; i++) { + dop[k][i] = pow(omega[k], 0.25) * exp(alpha[k] * (e_o[i] - Fermi)); + } + for (int a = 0; a < nvir; a++) { + dvp[k][a] = pow(omega[k], 0.25) * exp(-alpha[k] * (e_v[a + vir_start_] - Fermi)); + } + // for (int i = 0; i < nocc; i++) { + // for (int a = 0; a < nvir; a++) { + // dovp[k][i * nvir + a] = pow(omega[k], 0.5) * exp(-alpha[k] * (e_v[a] - e_o[i])); + // } + // } + } + + // denominator_occ_->print(); + // denominator_vir_->print(); + + delete[] alpha; + delete[] omega; + delete[] R_availp; + outfile->Printf("\n ==>Finish: FORTE Laplace Denominator <==\n\n"); +} + +void LaplaceDenominator::decompose_cavv() { + outfile->Printf("\n\n ==> FORTE Laplace Denominator (CAVV) <==\n\n"); + int nocc = eps_occ_->dimpi()[0]; + int nvir = eps_vir_->dimpi()[0]; + int nact = eps_act_->dimpi()[0]; + + // double E_LOMO = eps_occ_->get(0, 0); + // double E_HOMO = eps_occ_->get(0, nocc - 1); + // double E_LUMO = eps_vir_->get(0, 0); + // double E_HUMO = eps_vir_->get(0, nvir - 1); + + // double A = 2.0 * (E_LUMO - E_HOMO); // min + // double B = 2.0 * (E_HUMO - E_LOMO); // max + vir_start_ = 0; + + double E_vir_max = eps_vir_->get(0, nvir - 1); + //double E_vir_min = eps_vir_->get(0, 0); + double E_occ_max = eps_occ_->get(0, nocc - 1); + double E_occ_min = eps_occ_->get(0, 0); + double E_act_max = eps_act_->get(0, nact - 1); + double E_act_min = eps_act_->get(0, 0); + + for (int i = 0; i < nvir; i++) { + double E_vir_min_test = eps_vir_->get(0, i); + double E_test = 2 * E_vir_min_test - E_occ_max - E_act_max; + if (E_test < vir_tol_) { + vir_start_ += 1; + } else { + break; + } + } + + double E_vir_min = eps_vir_->get(0, vir_start_); + double A = 2 * E_vir_min - E_occ_max - E_act_max; + double B = 2 * E_vir_max - E_occ_min - E_act_min; + + nvir -= vir_start_; + + double R = B / A; + // double R = B; // no scale + + outfile->Printf(" The smallest denominator: %f. \n", A); + outfile->Printf(" The largest denominator: %f. \n", B); + outfile->Printf(" Number of unselected virtual orbitals: %d. \n", vir_start_); + + // Pick appropriate quadrature file and read contents + std::string PSIDATADIR = Process::environment.get_datadir(); + std::string err_table_filename = PSIDATADIR + "/quadratures/1_x/error.bin"; + std::string R_filename = PSIDATADIR + "/quadratures/1_x/R_avail.bin"; + + std::ifstream err_table_file(err_table_filename.c_str(), std::ios::in | std::ios::binary); + std::ifstream R_avail_file(R_filename.c_str(), std::ios::in | std::ios::binary); + + if (!err_table_file) + throw PSIEXCEPTION( + "LaplaceQuadrature: Cannot locate error property file for quadrature rules (should be " + "PSIDATADIR/quadratures/1_x/error.bin)"); + if (!R_avail_file) + throw PSIEXCEPTION( + "LaplaceQuadrature: Cannot locate R property file for quadrature rules (should be " + "PSIDATADIR/quadratures/1_x/R_avail.bin)"); + + int nk = 53; + int nR = 99; + + // Read in the R available + auto *R_availp = new double[nR]; + R_avail_file.read((char *)R_availp, nR * sizeof(double)); + + auto err_table = std::make_shared("Error Table (nR x nk)", nR, nk); + double **err_tablep = err_table->pointer(); + err_table_file.read((char *)err_tablep[0], static_cast (nR) * nk * sizeof(double)); + + R_avail_file.close(); + err_table_file.close(); + + // for (int r2 = 0; r2 < nR; r2++) + // outfile->Printf( " R[%4d] = %20.14E\n", r2+1, R_availp[r2]); + // err_table->print(); + + int indR; + for (indR = 0; indR < nR; indR++) { + if (R < R_availp[indR]) break; + } + if (indR == nR) { + // TODO: Relax this + throw PSIEXCEPTION( + "Laplace Quadrature requested for (E_HUMO - E_LOMO)/(E_LUMO-E_HOMO) > 7.0 * 10^12, quadratures are not " + "designed for this range."); + } + + double accuracy; + int k, r; + bool found = false; + for (k = 0; k < nk; k++) { + for (r = indR; r < nR; r++) { + double err = err_tablep[r][k]; + if (err != 0.0 && err < delta_) { + accuracy = err; + found = true; + break; + } + } + if (found) break; + } + + if (!found) { + throw PSIEXCEPTION("Laplace Quadrature rule could not be found with specified accuracy for this system"); + } + + nvector_ = k + 1; + + // A bit hacky, but OK + int exponent = (int)floor(log(R_availp[r]) / log(10.0)); + int mantissa = (int)round(R_availp[r] / pow(10.0, exponent)); + if (mantissa == 10) { + exponent++; + mantissa = 1; + } + + std::stringstream st; + st << std::setfill('0'); + st << "1_xk" << std::setw(2) << nvector_; + st << "_" << mantissa; + st << "E" << exponent; + + std::string quadfile = PSIDATADIR + "/quadratures/1_x/" + st.str().c_str(); + + outfile->Printf(" This system has an intrinsic R = (E_HUMO - E_LOMO)/(E_LUMO - E_HOMO) of %7.4E.\n", R); + outfile->Printf(" A %d point minimax quadrature with R of %1.0E will be used for the denominator.\n", nvector_, + R_availp[r]); + outfile->Printf(" The worst-case Chebyshev norm for this quadrature rule is %7.4E.\n", accuracy); + outfile->Printf(" Quadrature rule read from file %s.\n\n", quadfile.c_str()); + + // The quadrature is defined as \omega_v exp(-\alpha_v x) = 1/x + auto *alpha = new double[nvector_]; + auto *omega = new double[nvector_]; + + std::vector lines; + std::string text; + std::ifstream infile(quadfile.c_str()); + if (!infile) throw PSIEXCEPTION("LaplaceDenominator: Unable to open quadrature rule file: " + quadfile); + while (infile.good()) { + std::getline(infile, text); + lines.push_back(text); + } + +#define NUMBER "((?:[-+]?\\d*\\.\\d+(?:[DdEe][-+]?\\d+)?)|(?:[-+]?\\d+\\.\\d*(?:[DdEe][-+]?\\d+)?))" + std::regex numberline("^\\s*(" NUMBER ").*"); + std::smatch what; + + // We'll be rigorous, the files are extremely well defined + int lineno = 0; + for (int index = 0; index < nvector_; index++) { + std::string line = lines[lineno++]; + if (!std::regex_match(line, what, numberline)) + throw PSIEXCEPTION("LaplaceDenominator: Unable to read grid file line: \n" + line); + if (!from_string(omega[index], what[1], std::dec)) + throw PSIEXCEPTION("LaplaceDenominator: Unable to convert grid file line: \n" + line); + } + for (int index = 0; index < nvector_; index++) { + std::string line = lines[lineno++]; + if (!std::regex_match(line, what, numberline)) + throw PSIEXCEPTION("LaplaceDenominator: Unable to read grid file line: \n" + line); + if (!from_string(alpha[index], what[1], std::dec)) + throw PSIEXCEPTION("LaplaceDenominator: Unable to convert grid file line: \n" + line); + } + + // for (int k = 0; k < nvector_; k++) + // printf(" %24.16E, %24.16E\n", omega[k], alpha[k]); + + // Cast weights back to problem size + for (int k = 0; k < nvector_; k++) { + alpha[k] /= A; + omega[k] /= A; + } + + // // Fermi-level + double Fermi = (E_vir_min + E_act_max)/2; //???? How to calculate fermi-level in this case. + + denominator_occ_ = std::make_shared("Occupied Laplace Delta Tensor", nvector_, nocc); + denominator_vir_ = std::make_shared("Virtual Laplace Delta Tensor", nvector_, nvir); + denominator_act_ = std::make_shared("Active Laplace Delta Tensor", nvector_, nact); + + double **dop = denominator_occ_->pointer(); + double **dvp = denominator_vir_->pointer(); + double **dap = denominator_act_->pointer(); + + double *e_o = eps_occ_->pointer(); + double *e_v = eps_vir_->pointer(); + double *e_a = eps_act_->pointer(); + + for (int k = 0; k < nvector_; k++) { + for (int i = 0; i < nocc; i++) { + dop[k][i] = pow(omega[k], 0.25) * exp(alpha[k] * (e_o[i] - Fermi)); + } + for (int a = 0; a < nvir; a++) { + dvp[k][a] = pow(omega[k], 0.25) * exp(-alpha[k] * (e_v[a + vir_start_] - Fermi)); + } + for (int u = 0; u < nact; u++) { + dap[k][u] = pow(omega[k], 0.25) * exp(alpha[k] * (e_a[u] - Fermi)); + } + } + + delete[] alpha; + delete[] omega; + delete[] R_availp; + outfile->Printf("\n ==>Finish: FORTE Laplace Denominator (CAVV) <==\n\n"); +} + +void LaplaceDenominator::decompose_ccav() { + outfile->Printf("\n\n ==> FORTE Laplace Denominator (CCAV) <==\n\n"); + int nocc = eps_occ_->dimpi()[0]; + int nvir = eps_vir_->dimpi()[0]; + int nact = eps_act_->dimpi()[0]; + + // double E_LOMO = eps_occ_->get(0, 0); + // double E_HOMO = eps_occ_->get(0, nocc - 1); + // double E_LUMO = eps_vir_->get(0, 0); + // double E_HUMO = eps_vir_->get(0, nvir - 1); + + // double A = 2.0 * (E_LUMO - E_HOMO); // min + // double B = 2.0 * (E_HUMO - E_LOMO); // max + vir_start_ = 0; + + double E_vir_max = eps_vir_->get(0, nvir - 1); + //double E_vir_min = eps_vir_->get(0, 0); + double E_occ_max = eps_occ_->get(0, nocc - 1); + double E_occ_min = eps_occ_->get(0, 0); + double E_act_max = eps_act_->get(0, nact - 1); + double E_act_min = eps_act_->get(0, 0); + + for (int i = 0; i < nvir; i++) { + double E_vir_min_test = eps_vir_->get(0, i); + double E_test = E_vir_min_test + E_act_min - 2 * E_occ_max; + if (E_test < vir_tol_) { + vir_start_ += 1; + } else { + break; + } + } + + double E_vir_min = eps_vir_->get(0, vir_start_); + + double A = E_vir_min + E_act_min - 2 * E_occ_max; + double B = E_vir_max + E_act_max - 2 * E_occ_min; + + double R = B / A; + // double R = B; // no scale + + nvir -= vir_start_; + + outfile->Printf(" The smallest denominator: %f. \n", A); + outfile->Printf(" The largest denominator: %f. \n", B); + outfile->Printf(" Number of unselected virtual orbitals: %d. \n", vir_start_); + + // Pick appropriate quadrature file and read contents + std::string PSIDATADIR = Process::environment.get_datadir(); + std::string err_table_filename = PSIDATADIR + "/quadratures/1_x/error.bin"; + std::string R_filename = PSIDATADIR + "/quadratures/1_x/R_avail.bin"; + + std::ifstream err_table_file(err_table_filename.c_str(), std::ios::in | std::ios::binary); + std::ifstream R_avail_file(R_filename.c_str(), std::ios::in | std::ios::binary); + + if (!err_table_file) + throw PSIEXCEPTION( + "LaplaceQuadrature: Cannot locate error property file for quadrature rules (should be " + "PSIDATADIR/quadratures/1_x/error.bin)"); + if (!R_avail_file) + throw PSIEXCEPTION( + "LaplaceQuadrature: Cannot locate R property file for quadrature rules (should be " + "PSIDATADIR/quadratures/1_x/R_avail.bin)"); + + int nk = 53; + int nR = 99; + + // Read in the R available + auto *R_availp = new double[nR]; + R_avail_file.read((char *)R_availp, nR * sizeof(double)); + + auto err_table = std::make_shared("Error Table (nR x nk)", nR, nk); + double **err_tablep = err_table->pointer(); + err_table_file.read((char *)err_tablep[0], static_cast (nR) * nk * sizeof(double)); + + R_avail_file.close(); + err_table_file.close(); + + // for (int r2 = 0; r2 < nR; r2++) + // outfile->Printf( " R[%4d] = %20.14E\n", r2+1, R_availp[r2]); + // err_table->print(); + + int indR; + for (indR = 0; indR < nR; indR++) { + if (R < R_availp[indR]) break; + } + if (indR == nR) { + // TODO: Relax this + throw PSIEXCEPTION( + "Laplace Quadrature requested for (E_HUMO - E_LOMO)/(E_LUMO-E_HOMO) > 7.0 * 10^12, quadratures are not " + "designed for this range."); + } + + double accuracy; + int k, r; + bool found = false; + for (k = 0; k < nk; k++) { + for (r = indR; r < nR; r++) { + double err = err_tablep[r][k]; + if (err != 0.0 && err < delta_) { + accuracy = err; + found = true; + break; + } + } + if (found) break; + } + + if (!found) { + throw PSIEXCEPTION("Laplace Quadrature rule could not be found with specified accuracy for this system"); + } + + nvector_ = k + 1; + + // A bit hacky, but OK + int exponent = (int)floor(log(R_availp[r]) / log(10.0)); + int mantissa = (int)round(R_availp[r] / pow(10.0, exponent)); + if (mantissa == 10) { + exponent++; + mantissa = 1; + } + + std::stringstream st; + st << std::setfill('0'); + st << "1_xk" << std::setw(2) << nvector_; + st << "_" << mantissa; + st << "E" << exponent; + + std::string quadfile = PSIDATADIR + "/quadratures/1_x/" + st.str().c_str(); + + outfile->Printf(" This system has an intrinsic R = (E_HUMO - E_LOMO)/(E_LUMO - E_HOMO) of %7.4E.\n", R); + outfile->Printf(" A %d point minimax quadrature with R of %1.0E will be used for the denominator.\n", nvector_, + R_availp[r]); + outfile->Printf(" The worst-case Chebyshev norm for this quadrature rule is %7.4E.\n", accuracy); + outfile->Printf(" Quadrature rule read from file %s.\n\n", quadfile.c_str()); + + // The quadrature is defined as \omega_v exp(-\alpha_v x) = 1/x + auto *alpha = new double[nvector_]; + auto *omega = new double[nvector_]; + + std::vector lines; + std::string text; + std::ifstream infile(quadfile.c_str()); + if (!infile) throw PSIEXCEPTION("LaplaceDenominator: Unable to open quadrature rule file: " + quadfile); + while (infile.good()) { + std::getline(infile, text); + lines.push_back(text); + } + +#define NUMBER "((?:[-+]?\\d*\\.\\d+(?:[DdEe][-+]?\\d+)?)|(?:[-+]?\\d+\\.\\d*(?:[DdEe][-+]?\\d+)?))" + std::regex numberline("^\\s*(" NUMBER ").*"); + std::smatch what; + + // We'll be rigorous, the files are extremely well defined + int lineno = 0; + for (int index = 0; index < nvector_; index++) { + std::string line = lines[lineno++]; + if (!std::regex_match(line, what, numberline)) + throw PSIEXCEPTION("LaplaceDenominator: Unable to read grid file line: \n" + line); + if (!from_string(omega[index], what[1], std::dec)) + throw PSIEXCEPTION("LaplaceDenominator: Unable to convert grid file line: \n" + line); + } + for (int index = 0; index < nvector_; index++) { + std::string line = lines[lineno++]; + if (!std::regex_match(line, what, numberline)) + throw PSIEXCEPTION("LaplaceDenominator: Unable to read grid file line: \n" + line); + if (!from_string(alpha[index], what[1], std::dec)) + throw PSIEXCEPTION("LaplaceDenominator: Unable to convert grid file line: \n" + line); + } + + // for (int k = 0; k < nvector_; k++) + // printf(" %24.16E, %24.16E\n", omega[k], alpha[k]); + + // Cast weights back to problem size + for (int k = 0; k < nvector_; k++) { + alpha[k] /= A; + omega[k] /= A; + } + + denominator_occ_ = std::make_shared("Occupied Laplace Delta Tensor", nvector_, nocc); + denominator_vir_ = std::make_shared("Virtual Laplace Delta Tensor", nvector_, nvir); + denominator_act_ = std::make_shared("Active Laplace Delta Tensor", nvector_, nact); + + double **dop = denominator_occ_->pointer(); + double **dvp = denominator_vir_->pointer(); + double **dap = denominator_act_->pointer(); + + double *e_o = eps_occ_->pointer(); + double *e_v = eps_vir_->pointer(); + double *e_a = eps_act_->pointer(); + + for (int k = 0; k < nvector_; k++) { + for (int i = 0; i < nocc; i++) { + dop[k][i] = pow(omega[k], 0.25) * exp(alpha[k] * e_o[i]); + } + for (int a = 0; a < nvir; a++) { + dvp[k][a] = pow(omega[k], 0.25) * exp(-alpha[k] * e_v[a + vir_start_]); + } + for (int u = 0; u < nact; u++) { + dap[k][u] = pow(omega[k], 0.25) * exp(-alpha[k] * e_a[u]); + } + } + + delete[] alpha; + delete[] omega; + delete[] R_availp; + outfile->Printf("\n ==>Finish: FORTE Laplace Denominator (CCAV) <==\n\n"); +} + +} // namespace forte \ No newline at end of file diff --git a/forte/orbital-helpers/LaplaceDenominator.h b/forte/orbital-helpers/LaplaceDenominator.h new file mode 100644 index 000000000..52c4ab69f --- /dev/null +++ b/forte/orbital-helpers/LaplaceDenominator.h @@ -0,0 +1,51 @@ +#ifndef __LAPLACEDENOMINATOR_H__ +#define __LAPLACEDENOMINATOR_H__ + +#include "psi4/libmints/matrix.h" +#include "psi4/libmints/vector.h" + +namespace forte { +class LaplaceDenominator { + protected: + // Fully split denominator (w in rows, i in columns) + psi::SharedMatrix denominator_occ_; + // Fully split denominator (w in rows, a in columns) + psi::SharedMatrix denominator_vir_; + // Fully split denominator (w in rows, u in columns) + psi::SharedMatrix denominator_act_; + + // Pointer to active occupied orbital eigenvalues + std::shared_ptr eps_occ_; + // Pointer to active virtual orbital eigenvalues + std::shared_ptr eps_vir_; + // Pointer to active orbital eigenvalues + std::shared_ptr eps_act_; + // Number of vectors required to obtain given accuracy + int nvector_; + // Maximum error norm allowed in denominator + double delta_; + // The first selected virtual orbital + int vir_start_ = 0; + + double vir_tol_; + + bool cavv_; + + void decompose_ccvv(); + + void decompose_cavv(); + + void decompose_ccav(); + + public: + LaplaceDenominator(std::shared_ptr eps_occ_, std::shared_ptr eps_vir, double delta, double vir_tol); + LaplaceDenominator(std::shared_ptr eps_occ, std::shared_ptr eps_act, std::shared_ptr eps_vir, double delta, bool cavv, double vir_tol); + ~LaplaceDenominator(); + psi::SharedMatrix denominator_occ() const { return denominator_occ_; } + psi::SharedMatrix denominator_vir() const { return denominator_vir_; } + psi::SharedMatrix denominator_act() const { return denominator_act_; } + int vir_start() const { return vir_start_; } +}; +} + +#endif \ No newline at end of file diff --git a/forte/orbital-helpers/ao_helper.cc b/forte/orbital-helpers/ao_helper.cc index 80d6807f1..0ae46f59e 100644 --- a/forte/orbital-helpers/ao_helper.cc +++ b/forte/orbital-helpers/ao_helper.cc @@ -33,107 +33,151 @@ #include "psi4/libpsi4util/PsiOutStream.h" #include "psi4/libpsi4util/process.h" -#include "psi4/lib3index/denominator.h" +#include "LaplaceDenominator.h" #include "psi4/libfock/jk.h" +#include "Laplace.h" #include "ao_helper.h" +#include + using namespace psi; namespace forte { - -AtomicOrbitalHelper::AtomicOrbitalHelper(std::shared_ptr CMO, - std::shared_ptr eps_occ, - std::shared_ptr eps_vir, - double laplace_tolerance) - : CMO_(CMO), eps_rdocc_(eps_occ), eps_virtual_(eps_vir), laplace_tolerance_(laplace_tolerance) { - psi::LaplaceDenominator laplace(eps_rdocc_, eps_virtual_, laplace_tolerance_); + +AtomicOrbitalHelper::AtomicOrbitalHelper(psi::SharedMatrix CMO, psi::SharedVector eps_occ, + psi::SharedVector eps_vir, double laplace_tolerance, + int shift, int nfrozen, double vir_tol) + : CMO_(CMO), eps_rdocc_(eps_occ), eps_virtual_(eps_vir), laplace_tolerance_(laplace_tolerance), + shift_(shift), nfrozen_(nfrozen) { + LaplaceDenominator laplace(eps_rdocc_, eps_virtual_, laplace_tolerance_, vir_tol); Occupied_Laplace_ = laplace.denominator_occ(); Virtual_Laplace_ = laplace.denominator_vir(); weights_ = Occupied_Laplace_->rowspi()[0]; nrdocc_ = eps_rdocc_->dim(); - nvir_ = eps_virtual_->dim(); + nvir_ = Virtual_Laplace_->colspi()[0]; nbf_ = CMO_->rowspi()[0]; - shift_ = 0; + vir_start_ = laplace.vir_start(); + shift_ += vir_start_; } -AtomicOrbitalHelper::AtomicOrbitalHelper(std::shared_ptr CMO, - std::shared_ptr eps_occ, - std::shared_ptr eps_vir, - double laplace_tolerance, int shift) - : CMO_(CMO), eps_rdocc_(eps_occ), eps_virtual_(eps_vir), laplace_tolerance_(laplace_tolerance), - shift_(shift) { - psi::LaplaceDenominator laplace(eps_rdocc_, eps_virtual_, laplace_tolerance_); + +AtomicOrbitalHelper::AtomicOrbitalHelper(psi::SharedMatrix CMO, psi::SharedVector eps_occ, psi::SharedVector eps_act, + psi::SharedVector eps_vir, double laplace_tolerance, + int shift, int nfrozen, bool cavv, double vir_tol) + : CMO_(CMO), eps_rdocc_(eps_occ), eps_active_(eps_act), eps_virtual_(eps_vir), laplace_tolerance_(laplace_tolerance), + shift_(shift), nfrozen_(nfrozen) { + LaplaceDenominator laplace(eps_rdocc_, eps_active_, eps_virtual_, laplace_tolerance_, cavv, vir_tol); Occupied_Laplace_ = laplace.denominator_occ(); Virtual_Laplace_ = laplace.denominator_vir(); + Active_Laplace_ = laplace.denominator_act(); weights_ = Occupied_Laplace_->rowspi()[0]; nrdocc_ = eps_rdocc_->dim(); - nvir_ = eps_virtual_->dim(); + nact_ = eps_active_->dim(); + nvir_ = Virtual_Laplace_->colspi()[0]; + vir_start_ = laplace.vir_start(); + shift_ += vir_start_; nbf_ = CMO_->rowspi()[0]; } -AtomicOrbitalHelper::~AtomicOrbitalHelper() { outfile->Printf("\n Done with AO helper class"); } -void AtomicOrbitalHelper::Compute_Psuedo_Density() { - int nmo_ = nbf_; - auto Xocc = std::make_shared("DensityOccupied", weights_, nbf_ * nbf_); - auto Yvir = std::make_shared("DensityVirtual", weights_, nmo_ * nmo_); +AtomicOrbitalHelper::~AtomicOrbitalHelper() { outfile->Printf("\n Done with AO helper class"); } + +void AtomicOrbitalHelper::Compute_Cholesky_Pseudo_Density() { + psi::SharedMatrix POcc_single(new psi::Matrix("Single_POcc", nbf_, nbf_)); + psi::SharedMatrix PVir_single(new psi::Matrix("Single_PVir", nbf_, nbf_)); + + double value_occ, value_vir = 0.0; for (int w = 0; w < weights_; w++) { for (int mu = 0; mu < nbf_; mu++) { for (int nu = 0; nu < nbf_; nu++) { - double value_occ = 0.0; for (int i = 0; i < nrdocc_; i++) { - value_occ += CMO_->get(mu, i) * CMO_->get(nu, i) * Occupied_Laplace_->get(w, i); + value_occ += CMO_->get(mu, i + nfrozen_) * CMO_->get(nu, i + nfrozen_) * Occupied_Laplace_->get(w, i); } - Xocc->set(w, mu * nmo_ + nu, value_occ); - double value_vir = 0.0; + POcc_single->set(mu, nu, value_occ); for (int a = 0; a < nvir_; a++) { - value_vir += CMO_->get(mu, nrdocc_ + shift_ + a) * - CMO_->get(nu, nrdocc_ + shift_ + a) * Virtual_Laplace_->get(w, a); + value_vir += CMO_->get(mu, nfrozen_ + nrdocc_ + shift_ + a) * + CMO_->get(nu, nfrozen_ + nrdocc_ + shift_ + a) * Virtual_Laplace_->get(w, a); } - Yvir->set(w, mu * nmo_ + nu, value_vir); + PVir_single->set(mu, nu, value_vir); + value_occ = 0.0; + value_vir = 0.0; } } + psi::SharedMatrix LOcc = POcc_single->partial_cholesky_factorize(1e-10); + psi::SharedMatrix LVir = PVir_single->partial_cholesky_factorize(1e-10); + LOcc_list_.push_back(LOcc); + LVir_list_.push_back(LVir); + //POcc_list_.push_back(POcc_single); + //PVir_list_.push_back(PVir_single); + //n_pseudo_occ_list_.push_back(LOcc->coldim()); + //n_pseudo_vir_list_.push_back(LVir->coldim()); } - POcc_ = Xocc->clone(); - PVir_ = Yvir->clone(); } -void AtomicOrbitalHelper::Compute_AO_Screen(std::shared_ptr& primary) { - auto integral = std::make_shared(primary, primary, primary, primary); - auto twobodyaoints = std::shared_ptr(integral->eri()); - // ERISieve sieve(primary, 1e-10); - // auto my_function_pair_values = twobodyaoints->function_pair_values(); - auto AO_Screen = std::make_shared("Z", nbf_, nbf_); - for (int mu = 0; mu < nbf_; mu++) - for (int nu = 0; nu < nbf_; nu++) { - // AO_Screen->set(mu, nu, my_function_pair_values[mu * nbf_ + nu]); - double value = std::sqrt(twobodyaoints->function_ceiling2(mu, nu, mu, nu)); - AO_Screen->set(mu, nu, value); - } - AO_Screen_ = AO_Screen; - AO_Screen_->set_name("ScwartzAOInts"); -} -void AtomicOrbitalHelper::Estimate_TransAO_Screen(std::shared_ptr& primary, - std::shared_ptr& auxiliary) { - Compute_Psuedo_Density(); - MemDFJK jk(primary, auxiliary, psi::Process::environment.options); - jk.initialize(); - jk.compute(); - auto AO_Trans_Screen = std::make_shared("AOTrans", weights_, nbf_ * nbf_); +void AtomicOrbitalHelper::Compute_Cholesky_Active_Density(psi::SharedMatrix RDM) { + psi::SharedMatrix PAct_single(new psi::Matrix("Single_PAct", nbf_, nbf_)); + std::vector Act_idx(nact_); + std::iota(Act_idx.begin(), Act_idx.end(), nfrozen_+nrdocc_); + psi::SharedMatrix C_Act = submatrix_cols(*CMO_, Act_idx); + psi::SharedMatrix D_Act = psi::linalg::doublet(C_Act, RDM, false, false); + double value_act = 0.0; for (int w = 0; w < weights_; w++) { - auto COcc = std::make_shared("COcc", nbf_, nbf_); - auto CVir = std::make_shared("COcc", nbf_, nbf_); - for (int mu = 0; mu < nbf_; mu++) + for (int mu = 0; mu < nbf_; mu++) { for (int nu = 0; nu < nbf_; nu++) { - COcc->set(mu, nu, POcc_->get(w, mu * nbf_ + nu)); - CVir->set(mu, nu, PVir_->get(w, mu * nbf_ + nu)); + for (int u = 0; u < nact_; u++) { + value_act += D_Act->get(mu, u) * C_Act->get(nu, u) * Active_Laplace_->get(w, u); + } + PAct_single->set(mu, nu, value_act); + value_act = 0.0; } + } + psi::SharedMatrix LAct = PAct_single->partial_cholesky_factorize(1e-10); + LAct_list_.push_back(LAct); + //n_pseudo_vir_list_.push_back(LAct->coldim()); + } +} + +void AtomicOrbitalHelper::Compute_Cholesky_Density() { + std::vector Occ_idx(nrdocc_); + std::iota(Occ_idx.begin(), Occ_idx.end(), nfrozen_); + psi::SharedMatrix C_Occ = submatrix_cols(*CMO_, Occ_idx); + POcc_real_ = psi::linalg::doublet(C_Occ, C_Occ, false, true); + L_Occ_real_ = POcc_real_->partial_cholesky_factorize(1e-10); +} +/* +void AtomicOrbitalHelper::Householder_QR() { + std::vector Occ_idx(nrdocc_); + std::iota(Occ_idx.begin(), Occ_idx.end(), 0); + psi::SharedMatrix C_Occ = submatrix_cols(*CMO_, Occ_idx); + psi::SharedMatrix C_Occ_transpose = C_Occ->transpose(); + /// Copy C_Occ matrix to R + double** C_T = C_Occ_transpose->pointer(); + auto R = std::make_shared("R", nrdocc_, C_Occ_transpose->coldim()); + double** Rp = R->pointer(); + C_DCOPY(nrdocc_ * C_Occ_transpose->coldim(), C_T[0], 1, Rp[0], 1); + + /// QR decomposition + std::vector tau(nrdocc_); + //std::vector jpvt(C_Occ_transpose->coldim()); + + // First, find out how much workspace to provide + // Optimal size of work vector is written to work_size + double work_size; + //C_DGEQP3(nrdocc_, C_Occ_transpose->coldim(), Rp[0], nrdocc_, jpvt.data(), tau.data(), &work_size, -1); + C_DGEQRF(nrdocc_, C_Occ_transpose->coldim(), Rp[0], nrdocc_, tau.data(), &work_size, -1); + // Now, do the QR decomposition + int lwork = (int)work_size; + std::vector work(lwork); + //C_DGEQP3(nrdocc_, C_Occ_transpose->coldim(), Rp[0], nrdocc_, jpvt.data(), tau.data(), work.data(), lwork); + C_DGEQRF(nrdocc_, C_Occ_transpose->coldim(), Rp[0], nrdocc_, tau.data(), work.data(), lwork); - auto iaia_w = jk.iaia(COcc, CVir); - for (int mu = 0; mu < nbf_; mu++) - for (int nu = 0; nu < nbf_; nu++) - AO_Trans_Screen->set(w, mu * nbf_ + nu, iaia_w->get(mu * nbf_ + nu)); + // Put R in the upper triangle where it belongs + for (int i = 1; i < nrdocc_; i++) { + for (int j = 0; j < i; j++) { + Rp[j][i] = 0.0; + } } - TransAO_Screen_ = AO_Trans_Screen; - TransAO_Screen_->set_name("(u_b {b}^v | u_b {b}^v)"); + + R->transpose()->print(); } -} // namespace forte +*/ +} // namespace forte \ No newline at end of file diff --git a/forte/orbital-helpers/ao_helper.h b/forte/orbital-helpers/ao_helper.h index 74feb0ea8..5745fc331 100644 --- a/forte/orbital-helpers/ao_helper.h +++ b/forte/orbital-helpers/ao_helper.h @@ -35,47 +35,82 @@ namespace forte { class AtomicOrbitalHelper { protected: - std::shared_ptr AO_Screen_; - std::shared_ptr TransAO_Screen_; + psi::SharedMatrix CMO_; + psi::SharedVector eps_rdocc_; + psi::SharedVector eps_virtual_; + psi::SharedVector eps_active_; - std::shared_ptr CMO_; - std::shared_ptr eps_rdocc_; - std::shared_ptr eps_virtual_; + std::vector LOcc_list_; + std::vector LVir_list_; + std::vector LAct_list_; + + std::vector POcc_list_; + std::vector PVir_list_; + + std::vector n_pseudo_occ_list_; + std::vector n_pseudo_vir_list_; + std::vector n_pseudo_act_list_; + + psi::SharedMatrix L_Occ_real_; + //psi::SharedMatrix L_Vir_real_; + psi::SharedMatrix POcc_real_; + // psi::SharedMatrix PVir_real_; - std::shared_ptr POcc_; - std::shared_ptr PVir_; - void Compute_Psuedo_Density(); // LaplaceDenominator Laplace_; - std::shared_ptr Occupied_Laplace_; - std::shared_ptr Virtual_Laplace_; + psi::SharedMatrix Occupied_Laplace_; + psi::SharedMatrix Virtual_Laplace_; + psi::SharedMatrix Active_Laplace_; + double laplace_tolerance_ = 1e-10; int weights_; int nbf_; int nrdocc_; int nvir_; + int nact_; /// How many orbitals does it take to go from occupied to virtual (ie should /// be active) int shift_; + int vir_start_; + int nfrozen_; public: - std::shared_ptr AO_Screen() { return AO_Screen_; } - std::shared_ptr TransAO_Screen() { return TransAO_Screen_; } - std::shared_ptr Occupied_Laplace() { return Occupied_Laplace_; } - std::shared_ptr Virtual_Laplace() { return Virtual_Laplace_; } - std::shared_ptr POcc() { return POcc_; } - std::shared_ptr PVir() { return PVir_; } + psi::SharedMatrix Occupied_Laplace() { return Occupied_Laplace_; } + psi::SharedMatrix Virtual_Laplace() { return Virtual_Laplace_; } + + std::vector POcc_list() { return POcc_list_; } + std::vector PVir_list() { return PVir_list_; } + + std::vector LOcc_list() { return LOcc_list_; } + std::vector LVir_list() { return LVir_list_; } + std::vector LAct_list() { return LAct_list_; } + + std::vector n_pseudo_occ_list() { return n_pseudo_occ_list_; } + std::vector n_pseudo_vir_list() { return n_pseudo_vir_list_; } + std::vector n_pseudo_act_list() { return n_pseudo_act_list_; } + + psi::SharedMatrix L_Occ_real() { return L_Occ_real_; } + + psi::SharedMatrix POcc_real() { return POcc_real_; } + int Weights() { return weights_; } - AtomicOrbitalHelper(std::shared_ptr CMO, std::shared_ptr eps_occ, - std::shared_ptr eps_vir, double laplace_tolerance); - AtomicOrbitalHelper(std::shared_ptr CMO, std::shared_ptr eps_occ, - std::shared_ptr eps_vir, double laplace_tolerance, int shift); - /// Compute (mu nu | mu nu)^{(1/2)} - void Compute_AO_Screen(std::shared_ptr& primary); - void Estimate_TransAO_Screen(std::shared_ptr& primary, - std::shared_ptr& auxiliary); + int vir_start() { return vir_start_; } + + + AtomicOrbitalHelper(psi::SharedMatrix CMO, psi::SharedVector eps_occ, psi::SharedVector eps_vir, + double laplace_tolerance, int shift, int nfrozen, double vir_tol); + + AtomicOrbitalHelper(psi::SharedMatrix CMO, psi::SharedVector eps_occ, psi::SharedVector eps_act, + psi::SharedVector eps_vir, double laplace_tolerance, + int shift, int nfrozen, bool cavv, double vir_tol); + + void Compute_Cholesky_Density(); + void Compute_Cholesky_Pseudo_Density(); + void Compute_Cholesky_Active_Density(psi::SharedMatrix RDM); + //void Householder_QR(); + //void Compute_L_Directly(); ~AtomicOrbitalHelper(); }; diff --git a/forte/register_forte_options.py b/forte/register_forte_options.py index 769f5d76b..58f0d03d6 100644 --- a/forte/register_forte_options.py +++ b/forte/register_forte_options.py @@ -856,11 +856,57 @@ def register_dsrg_options(options): "BATCH_CORE_MPI", "BATCH_CORE_REP", "BATCH_VIRTUAL_REP", + "LT-DSRG" ], "Algorithm to compute the CCVV term in DSRG-MRPT2 (only in three-dsrg-mrpt2 code)", ) options.add_bool("AO_DSRG_MRPT2", False, "Do AO-DSRG-MRPT2 if true (not available)") + + options.add_double("THETA_NB", 0, "theta_NB") + + options.add_double("THETA_NB_IAP", 0, "theta_NB_IAP") + + options.add_double("THETA_IJ", 0, "theta_ij") + + options.add_double("THETA_SCHWARZ", 0, "theta_schwarz") + + options.add_double("OMEGA", 0, "Omega for erfc attenuated coulomb operator") + + options.add_double("LAPLACE_THRESHOLD", 1e-10, "laplace threshold") + + options.add_double("VIR_TOL", 0, "vir_tol") + +# CAVV thresholds + options.add_double("THETA_NB_CAVV", 0, "theta_NB") + + options.add_double("THETA_XNB_CAVV", 0, "theta_XNB") + + options.add_double("THETA_NB_IAP_CAVV", 0, "theta_NB_IAP") + + options.add_double("THETA_NB_XAP_CAVV", 0, "theta_NB_XAP") + + options.add_double("THETA_IJ_CAVV", 0, "theta_ij") + + options.add_double("THETA_SCHWARZ_CAVV", 0, "theta_schwarz") +# CCAV thresholds + options.add_double("THETA_NB_CCAV", 0, "theta_NB") + + options.add_double("THETA_NB_IAP_CCAV", 0, "theta_NB_IAP") + + options.add_double("THETA_IJ_CCAV", 0, "theta_ij") + + options.add_double("THETA_SCHWARZ_CCAV", 0, "theta_schwarz") + + options.add_bool("LAPLACE_CORE", True, "Use core or disk algorithm") + +# options.add_bool("LAPLACE_ONE_ACTIVE", False, "Use Laplace for one active section") + + options.add_bool("LAPLACE_CAVV", False, "Use Laplace for CAVV?") + + options.add_bool("LAPLACE_CCAV", False, "Use Laplace for CCAV? Should only be applied to really large active spaces") + + # options.add_double("THETA_NB_CAVV", 0, "theta_NB_cavv") options.add_int("CCVV_BATCH_NUMBER", -1, "Batches for CCVV_ALGORITHM")