From 11eb4ddbb62f4150a1458623cdd9fab08c831ed1 Mon Sep 17 00:00:00 2001 From: Shuhang Li Date: Tue, 12 Nov 2024 13:50:51 -0500 Subject: [PATCH] ref_relax --- forte/api/orbital_api.cc | 4 + forte/orbital-helpers/semi_canonicalize.cc | 74 ++++++++++++ forte/orbital-helpers/semi_canonicalize.h | 49 ++++++-- forte/proc/dsrg.py | 132 +++++++++------------ 4 files changed, 176 insertions(+), 83 deletions(-) diff --git a/forte/api/orbital_api.cc b/forte/api/orbital_api.cc index bfaeda282..15f3db4e9 100644 --- a/forte/api/orbital_api.cc +++ b/forte/api/orbital_api.cc @@ -78,8 +78,12 @@ void export_SemiCanonical(py::module& m) { "Semicanonicalize the orbitals and transform the integrals and reference") .def("Ua", &SemiCanonical::Ua, "Return the alpha rotation matrix") .def("Ub", &SemiCanonical::Ub, "Return the alpha rotation matrix") + .def("Ua_c", &SemiCanonical::Ua_c, "Return the alpha rotation matrix in the core space") + .def("Ub_c", &SemiCanonical::Ub_c, "Return the beta rotation matrix in the core space") .def("Ua_t", &SemiCanonical::Ua_t, "Return the alpha rotation matrix in the active space") .def("Ub_t", &SemiCanonical::Ub_t, "Return the beta rotation matrix in the active space") + .def("Ua_v", &SemiCanonical::Ua_v, "Return the alpha rotation matrix in the virtual space") + .def("Ub_v", &SemiCanonical::Ub_v, "Return the beta rotation matrix in the virtual space") .def("fix_orbital_success", &SemiCanonical::fix_orbital_success, "Return if the orbital ordering and phases are fixed successfully"); } diff --git a/forte/orbital-helpers/semi_canonicalize.cc b/forte/orbital-helpers/semi_canonicalize.cc index 462282982..fbb1a1159 100644 --- a/forte/orbital-helpers/semi_canonicalize.cc +++ b/forte/orbital-helpers/semi_canonicalize.cc @@ -74,16 +74,26 @@ void SemiCanonical::read_options(const std::shared_ptr& options) { void SemiCanonical::startup() { nirrep_ = mo_space_info_->nirrep(); nmopi_ = mo_space_info_->dimension("ALL"); + ncore_ = mo_space_info_->size("CORE"); nact_ = mo_space_info_->size("ACTIVE"); + nvirt_ = mo_space_info_->size("VIRTUAL"); // Prepare orbital rotation matrix, which transforms all MOs Ua_ = std::make_shared("Ua", nmopi_, nmopi_); Ub_ = std::make_shared("Ub", nmopi_, nmopi_); + // Prepare orbital rotation matrix, which transforms only core MOs + Ua_c_ = ambit::Tensor::build(ambit::CoreTensor, "Ua", {ncore_, ncore_}); + Ub_c_ = ambit::Tensor::build(ambit::CoreTensor, "Ub", {ncore_, ncore_}); + // Prepare orbital rotation matrix, which transforms only active MOs Ua_t_ = ambit::Tensor::build(ambit::CoreTensor, "Ua", {nact_, nact_}); Ub_t_ = ambit::Tensor::build(ambit::CoreTensor, "Ub", {nact_, nact_}); + // Prepare orbital rotation matrix, which transforms only virtual MOs + Ua_v_ = ambit::Tensor::build(ambit::CoreTensor, "Ua", {nvirt_, nvirt_}); + Ub_v_ = ambit::Tensor::build(ambit::CoreTensor, "Ub", {nvirt_, nvirt_}); + // Initialize U to identity set_U_to_identity(); @@ -120,11 +130,23 @@ void SemiCanonical::set_U_to_identity() { Ua_->identity(); Ub_->identity(); + Ua_c_.iterate( + [&](const std::vector& i, double& value) { value = (i[0] == i[1]) ? 1.0 : 0.0; }); + + Ub_c_.iterate( + [&](const std::vector& i, double& value) { value = (i[0] == i[1]) ? 1.0 : 0.0; }); + Ua_t_.iterate( [&](const std::vector& i, double& value) { value = (i[0] == i[1]) ? 1.0 : 0.0; }); Ub_t_.iterate( [&](const std::vector& i, double& value) { value = (i[0] == i[1]) ? 1.0 : 0.0; }); + + Ua_v_.iterate( + [&](const std::vector& i, double& value) { value = (i[0] == i[1]) ? 1.0 : 0.0; }); + + Ub_v_.iterate( + [&](const std::vector& i, double& value) { value = (i[0] == i[1]) ? 1.0 : 0.0; }); } void SemiCanonical::semicanonicalize(std::shared_ptr rdms, const bool& build_fock, @@ -274,12 +296,64 @@ void SemiCanonical::build_transformation_matrices(const bool& semi) { // keep phase and order unchanged fix_orbital_success_ = ints_->fix_orbital_phases(Ua_, true); + // fill in Ua_c_ + fill_Ucore(Ua_, Ua_c_); // fill in Ua_t_ fill_Uactv(Ua_, Ua_t_); + // fill in Ua_v_ + fill_Uvirt(Ua_, Ua_v_); // pass to Ub Ub_->copy(Ua_); + Ub_c_.copy(Ua_c_); Ub_t_.copy(Ua_t_); + Ub_v_.copy(Ua_v_); +} + +void SemiCanonical::fill_Ucore(const std::shared_ptr& U, ambit::Tensor& Ut) { + auto core_names = mo_space_info_->composite_space_names()["CORE"]; + auto& Ut_data = Ut.data(); + + for (const std::string& name : core_names) { + auto size = mo_space_info_->size(name); + if (size == 0) + continue; + + auto pos = mo_space_info_->pos_in_space(name, "CORE"); + auto relative_mos = mo_space_info_->relative_mo(name); + for (size_t p = 0; p < size; ++p) { + const auto& [hp, np] = relative_mos[p]; + for (size_t q = 0; q < size; ++q) { + const auto& [hq, nq] = relative_mos[q]; + if (hp != hq) + continue; + Ut_data[pos[p] * ncore_ + pos[q]] = U->get(hp, np, nq); + } + } + } +} + +void SemiCanonical::fill_Uvirt(const std::shared_ptr& U, ambit::Tensor& Ut) { + auto virtual_names = mo_space_info_->composite_space_names()["VIRTUAL"]; + auto& Ut_data = Ut.data(); + + for (const std::string& name : virtual_names) { + auto size = mo_space_info_->size(name); + if (size == 0) + continue; + + auto pos = mo_space_info_->pos_in_space(name, "VIRTUAL"); + auto relative_mos = mo_space_info_->relative_mo(name); + for (size_t p = 0; p < size; ++p) { + const auto& [hp, np] = relative_mos[p]; + for (size_t q = 0; q < size; ++q) { + const auto& [hq, nq] = relative_mos[q]; + if (hp != hq) + continue; + Ut_data[pos[p] * nvirt_ + pos[q]] = U->get(hp, np, nq); + } + } + } } void SemiCanonical::fill_Uactv(const std::shared_ptr& U, ambit::Tensor& Ut) { diff --git a/forte/orbital-helpers/semi_canonicalize.h b/forte/orbital-helpers/semi_canonicalize.h index 5fc779537..9e325b2f5 100644 --- a/forte/orbital-helpers/semi_canonicalize.h +++ b/forte/orbital-helpers/semi_canonicalize.h @@ -42,16 +42,16 @@ namespace forte { class ForteOptions; /// @brief The SemiCanonical class -/// This class computes semi-canonical orbitals from the 1RDM and optionally transforms the integrals and RDMs -/// Semi-canonical orbitals are obtained by diagonalizing the Fock matrix in each orbital space separately -/// The class can also produce natural orbitals. These differ by the semi-canonical orbital only in the active space -/// where they are defined to be eigenvectors of the 1RDM -/// -/// The final orbitals are ordered by increasing energy within each irrep and space. Natural orbitals are ordered -/// by decreasing occupation number +/// This class computes semi-canonical orbitals from the 1RDM and optionally transforms the +/// integrals and RDMs Semi-canonical orbitals are obtained by diagonalizing the Fock matrix in each +/// orbital space separately The class can also produce natural orbitals. These differ by the +/// semi-canonical orbital only in the active space where they are defined to be eigenvectors of the +/// 1RDM +/// +/// The final orbitals are ordered by increasing energy within each irrep and space. Natural +/// orbitals are ordered by decreasing occupation number class SemiCanonical { public: - /// @brief SemiCanonical Constructor /// @param mo_space_info The MOSpaceInfo object /// @param ints The ForteIntegrals object @@ -78,12 +78,22 @@ class SemiCanonical { /// @return the beta rotation matrix std::shared_ptr Ub() { return Ub_; } + /// @return the alpha rotation matrix in the core space + ambit::Tensor Ua_c() const { return Ua_c_.clone(); } + /// @return the beta rotation matrix in the core space + ambit::Tensor Ub_c() const { return Ub_c_.clone(); } + /// @return the alpha rotation matrix in the active space ambit::Tensor Ua_t() const { return Ua_t_.clone(); } /// @return the beta rotation matrix in the active space ambit::Tensor Ub_t() const { return Ub_t_.clone(); } + /// @return the alpha rotation matrix in the virtual space + ambit::Tensor Ua_v() const { return Ua_v_.clone(); } + /// @return the beta rotation matrix in the virtual space + ambit::Tensor Ub_v() const { return Ub_v_.clone(); } + /// @return if the orbital ordering and phases are fixed successfully bool fix_orbital_success() const { return fix_orbital_success_; } @@ -117,9 +127,15 @@ class SemiCanonical { /// Offset of GAS orbitals within ACTIVE std::map actv_offsets_; + /// Number of core MOs + size_t ncore_; + /// Number of active MOs size_t nact_; + /// Number of virtual MOs + size_t nvirt_; + /// Number of irreps size_t nirrep_; @@ -127,11 +143,22 @@ class SemiCanonical { std::shared_ptr Ua_; /// Unitary matrix for beta orbital rotation std::shared_ptr Ub_; + + /// Unitary matrix for alpha orbital rotation in the core space + ambit::Tensor Ua_c_; + /// Unitary matrix for beta orbital rotation in the core space + ambit::Tensor Ub_c_; + /// Unitary matrix for alpha orbital rotation in the active space ambit::Tensor Ua_t_; /// Unitary matrix for beta orbital rotation in the active space ambit::Tensor Ub_t_; + /// Unitary matrix for alpha orbital rotation in the virtual space + ambit::Tensor Ua_v_; + /// Unitary matrix for beta orbital rotation in the virtual space + ambit::Tensor Ub_v_; + /// Set Ua_, Ub_, Ua_t_, and Ub_t_ to identity void set_U_to_identity(); @@ -153,9 +180,15 @@ class SemiCanonical { /// Builds unitary matrices used to diagonalize diagonal blocks of Fock void build_transformation_matrices(const bool& semi); + /// Fill ambit::Tensor Ua_c_ (Ub_c_) using std::shared_ptr Ua_ (Ub_) + void fill_Ucore(const std::shared_ptr& U, ambit::Tensor& Ut); + /// Fill ambit::Tensor Ua_t_ (Ub_t_) using std::shared_ptr Ua_ (Ub_) void fill_Uactv(const std::shared_ptr& U, ambit::Tensor& Ut); + /// Fill ambit::Tensor Ua_v_ (Ub_v_) using std::shared_ptr Ua_ (Ub_) + void fill_Uvirt(const std::shared_ptr& U, ambit::Tensor& Ut); + /// Successfully fix the orbital ordering and phases bool fix_orbital_success_; }; diff --git a/forte/proc/dsrg.py b/forte/proc/dsrg.py index 4aae8191d..04088a644 100644 --- a/forte/proc/dsrg.py +++ b/forte/proc/dsrg.py @@ -246,41 +246,11 @@ def compute_energy(self): if not self.Heff_implemented: self.relax_maxiter = 0 - # if self.options.get_bool('FULL_HBAR'): - # self.relax_maxiter = 0 - if self.options.get_bool("FULL_HBAR") and self.relax_maxiter == 0: - psi4.core.print_out("\n =>** Saving Full Hbar **<=\n") + psi4.core.print_out("\n =>** Saving Full Hbar (unrelax) **<=\n") Heff = self.dsrg_solver.compute_Heff_full() Heff_dict = forte.Heff_dict(Heff) np.savez("save_Hbar", **Heff_dict) - psi4.core.print_out("\n =>** Getting gamma1 **<=\n") - gamma1 = self.dsrg_solver.get_gamma1() - psi4.core.print_out("\n =>** Getting eta1 **<=\n") - eta1 = self.dsrg_solver.get_eta1() - psi4.core.print_out("\n =>** Getting lambda2 **<=\n") - lambda2 = self.dsrg_solver.get_lambda2() - psi4.core.print_out("\n =>** Getting lambda3 **<=\n") - lambda3 = self.dsrg_solver.get_lambda3() - gamma1_dict = forte.blocktensor_to_dict(gamma1) - eta_1_dict = forte.blocktensor_to_dict(eta1) - lambda2_dict = forte.blocktensor_to_dict(lambda2) - if self.solver_type in ["MRDSRG_SO", "MRDSRG-SO"]: - lambda3_dict = forte.blocktensor_to_dict(lambda3) - else: - lambda3_dict = forte.L3_dict(lambda3) - - np.savez("save_gamma1", **gamma1_dict) - np.savez("save_eta1", **eta_1_dict) - np.savez("save_lambda2", **lambda2_dict) - np.savez("save_lambda3", **lambda3_dict) - del gamma1, eta1, lambda2, lambda3, gamma1_dict, eta_1_dict, lambda2_dict, lambda3_dict - if self.options.get_str("FOURPDC") != "ZERO": - psi4.core.print_out("\n =>** Getting lambda4 **<=\n") - lambda4 = self.dsrg_solver.get_lambda4() - lambda4_dict = forte.L4_dict(lambda4) - np.savez("save_lambda4", **lambda4_dict) - del lambda4, lambda4_dict if self.solver_type not in ["MRDSRG_SO", "MRDSRG-SO"]: psi4.core.print_out("\n =>** Getting dipole integral **<=\n") @@ -305,6 +275,27 @@ def compute_energy(self): # so that the CI vectors are comparable before and after DSRG dressing. # However, the ForteIntegrals object and the dipole integrals always refer to the current semi-canonical basis. # so to compute the dipole moment correctly, we need to make the RDMs and orbital basis consistent + + if self.options.get_bool("FULL_HBAR") and n == self.relax_maxiter - 1: + psi4.core.print_out("\n =>** Saving Full Hbar (relax) **<=\n") + Heff = self.dsrg_solver.compute_Heff_full() + Heff_dict = forte.Heff_dict(Heff) + np.savez("save_Hbar", **Heff_dict) + + if self.solver_type not in ["MRDSRG_SO", "MRDSRG-SO"]: + psi4.core.print_out("\n =>** Getting dipole integral **<=\n") + Mbar0 = self.dsrg_solver.compute_Mbar0_full() + print(Mbar0) + np.save("Mbar0", Mbar0) + Mbar1 = self.dsrg_solver.compute_Mbar1_full() + Mbar2 = self.dsrg_solver.compute_Mbar2_full() + + for i in range(3): + np.savez(f"Mbar1_{i}", **forte.blocktensor_to_dict(Mbar1[i])) + np.savez(f"Mbar2_{i}", **forte.blocktensor_to_dict(Mbar2[i])) + + del Mbar0, Mbar1, Mbar2 + ints_dressed = self.dsrg_solver.compute_Heff_actv() if self.fno_pt2_Heff_shift is not None: ints_dressed.add(self.fno_pt2_Heff_shift, 1.0) @@ -434,52 +425,43 @@ def compute_energy(self): psi4.core.print_out(f"\n\n DSRG-MRPT2 FNO energy correction: {self.fno_pt2_energy_shift:20.15f}") psi4.core.print_out(f"\n DSRG-MRPT2 FNO corrected energy: {e_dsrg:20.15f}") - if self.options.get_bool("FULL_HBAR"): - self.rdms.rotate(self.Ua, self.Ub) - psi4.core.print_out("\n =>** Saving Full Hbar **<=\n") - Heff = self.dsrg_solver.compute_Heff_full() - Heff_dict = forte.Heff_dict(Heff) - np.savez("save_Hbar", **Heff_dict) - psi4.core.print_out("\n =>** Getting gamma1 **<=\n") - gamma1 = self.dsrg_solver.get_gamma1() - psi4.core.print_out("\n =>** Getting eta1 **<=\n") - eta1 = self.dsrg_solver.get_eta1() - psi4.core.print_out("\n =>** Getting lambda2 **<=\n") - lambda2 = self.dsrg_solver.get_lambda2() - psi4.core.print_out("\n =>** Getting lambda3 **<=\n") - lambda3 = self.dsrg_solver.get_lambda3() - gamma1_dict = forte.blocktensor_to_dict(gamma1) - eta_1_dict = forte.blocktensor_to_dict(eta1) - lambda2_dict = forte.blocktensor_to_dict(lambda2) - if self.solver_type in ["MRDSRG_SO", "MRDSRG-SO"]: - lambda3_dict = forte.blocktensor_to_dict(lambda3) - else: - lambda3_dict = forte.L3_dict(lambda3) - - np.savez("save_gamma1", **gamma1_dict) - np.savez("save_eta1", **eta_1_dict) - np.savez("save_lambda2", **lambda2_dict) - np.savez("save_lambda3", **lambda3_dict) - del gamma1, eta1, lambda2, lambda3, gamma1_dict, eta_1_dict, lambda2_dict, lambda3_dict - if self.options.get_str("FOURPDC") != "ZERO": - psi4.core.print_out("\n =>** Getting lambda4 **<=\n") - lambda4 = self.dsrg_solver.get_lambda4() - lambda4_dict = forte.L4_dict(lambda4) - np.savez("save_lambda4", **lambda4_dict) - del lambda4, lambda4_dict - if self.solver_type not in ["MRDSRG_SO", "MRDSRG-SO"]: - psi4.core.print_out("\n =>** Getting dipole integral **<=\n") - Mbar0 = self.dsrg_solver.compute_Mbar0_full() - print(Mbar0) - np.save("Mbar0", Mbar0) - Mbar1 = self.dsrg_solver.compute_Mbar1_full() - Mbar2 = self.dsrg_solver.compute_Mbar2_full() + if self.options.get_bool("FULL_HBAR"): + if self.relax_maxiter != 0: + self.rdms = self.active_space_solver.compute_average_rdms( + self.state_weights_map, self.max_rdm_level, self.rdm_type + ) + self.rdms.rotate(self.Ua, self.Ub) # To previous semi-canonical basis - for i in range(3): - np.savez(f"Mbar1_{i}", **forte.blocktensor_to_dict(Mbar1[i])) - np.savez(f"Mbar2_{i}", **forte.blocktensor_to_dict(Mbar2[i])) + self.make_dsrg_solver() - del Mbar0, Mbar1, Mbar2 + psi4.core.print_out("\n =>** Getting gamma1 **<=\n") + gamma1 = self.dsrg_solver.get_gamma1() + psi4.core.print_out("\n =>** Getting eta1 **<=\n") + eta1 = self.dsrg_solver.get_eta1() + psi4.core.print_out("\n =>** Getting lambda2 **<=\n") + lambda2 = self.dsrg_solver.get_lambda2() + psi4.core.print_out("\n =>** Getting lambda3 **<=\n") + lambda3 = self.dsrg_solver.get_lambda3() + + gamma1_dict = forte.blocktensor_to_dict(gamma1) + eta1_dict = forte.blocktensor_to_dict(eta1) + lambda2_dict = forte.blocktensor_to_dict(lambda2) + if self.solver_type in ["MRDSRG_SO", "MRDSRG-SO"]: + lambda3_dict = forte.blocktensor_to_dict(lambda3) + else: + lambda3_dict = forte.L3_dict(lambda3) + + np.savez("save_gamma1", **gamma1_dict) + np.savez("save_eta1", **eta1_dict) + np.savez("save_lambda2", **lambda2_dict) + np.savez("save_lambda3", **lambda3_dict) + del gamma1, eta1, lambda2, lambda3, gamma1_dict, eta1_dict, lambda2_dict, lambda3_dict + if self.options.get_str("FOURPDC") != "ZERO": + psi4.core.print_out("\n =>** Getting lambda4 **<=\n") + lambda4 = self.dsrg_solver.get_lambda4() + lambda4_dict = forte.L4_dict(lambda4) + np.savez("save_lambda4", **lambda4_dict) + del lambda4, lambda4_dict self.dsrg_cleanup()