diff --git a/forte/api/orbital_api.cc b/forte/api/orbital_api.cc index 9856917c7..8e614541e 100644 --- a/forte/api/orbital_api.cc +++ b/forte/api/orbital_api.cc @@ -88,12 +88,8 @@ 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/mrdsrg-spin-integrated/master_mrdsrg.cc b/forte/mrdsrg-spin-integrated/master_mrdsrg.cc index 9d8c0c229..b3be3ac7f 100644 --- a/forte/mrdsrg-spin-integrated/master_mrdsrg.cc +++ b/forte/mrdsrg-spin-integrated/master_mrdsrg.cc @@ -707,7 +707,6 @@ void MASTER_DSRG::deGNO_ints_full(const std::string& name, double& H0, BlockedTe L1h.block("CC").iterate( [&](const std::vector& i, double& value) { value = i[0] == i[1] ? 1.0 : 0.0; }); - // scalar from H1 double scalar1 = 0.0; scalar1 -= H1["ji"] * L1h["ij"]; diff --git a/forte/orbital-helpers/semi_canonicalize.cc b/forte/orbital-helpers/semi_canonicalize.cc index 66b8efc32..b8224c693 100644 --- a/forte/orbital-helpers/semi_canonicalize.cc +++ b/forte/orbital-helpers/semi_canonicalize.cc @@ -64,26 +64,16 @@ SemiCanonical::SemiCanonical(std::shared_ptr mo_space_info, 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(); @@ -122,23 +112,11 @@ 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, bool build_fock, @@ -299,64 +277,12 @@ void SemiCanonical::build_transformation_matrices(const bool& semi) { if (!inactive_mix_) 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_spaces_def().at("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_spaces_def().at("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 76bf9a84c..90b58c6d9 100644 --- a/forte/orbital-helpers/semi_canonicalize.h +++ b/forte/orbital-helpers/semi_canonicalize.h @@ -109,22 +109,12 @@ 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_; } @@ -158,15 +148,9 @@ 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_; @@ -175,21 +159,11 @@ class SemiCanonical { /// 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(); @@ -211,15 +185,9 @@ 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 43ac44bb6..1cb6c3e67 100644 --- a/forte/proc/dsrg.py +++ b/forte/proc/dsrg.py @@ -249,7 +249,7 @@ def compute_energy(self): if not self.Heff_implemented: self.relax_maxiter = 0 - if self.options.get_bool("FULL_HBAR") and self.solver_type in ["MRDSRG","SA-MRDSRG","SA_MRDSRG"]: + if self.options.get_bool("FULL_HBAR") and self.solver_type in ["MRDSRG","SA-MRDSRG","SA_MRDSRG"] and self.relax_maxiter == 0: psi4.core.print_out("\n =>** Saving Full Hbar (unrelaxed) **<=\n") Hbar1, Hbar2 = self.dsrg_solver.compute_Heff_full() np.savez("save_Hbar", **Hbar1, **Hbar2) @@ -267,7 +267,7 @@ def compute_energy(self): del Mbar0, Mbar1, Mbar2 - if self.options.get_bool("FULL_HBAR_DEGNO") and self.solver_type in ["SA-MRDSRG","SA_MRDSRG","MRDSRG"]: + if self.options.get_bool("FULL_HBAR_DEGNO") and self.solver_type in ["MRDSRG","SA-MRDSRG","SA_MRDSRG"] and self.relax_maxiter == 0: psi4.core.print_out("\n =>** Saving Full Hbar in de-normal-ordered basis (unrelaxed) **<=\n") Hbar1, Hbar2 = self.dsrg_solver.compute_Heff_full_degno() np.savez("save_Hbar_degno", **Hbar1, **Hbar2) @@ -283,16 +283,14 @@ def compute_energy(self): # 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.options.get_bool("FULL_HBAR") and self.solver_type in ["MRDSRG","SA-MRDSRG","SA_MRDSRG"] and n == self.relax_maxiter - 1: + psi4.core.print_out("\n =>** Saving Full Hbar (relaxed) **<=\n") + Hbar1, Hbar2 = self.dsrg_solver.compute_Heff_full() + np.savez("save_Hbar", **Hbar1, **Hbar2) 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() @@ -302,6 +300,11 @@ def compute_energy(self): np.savez(f"Mbar2_{i}", **Mbar2[i]) del Mbar0, Mbar1, Mbar2 + + if self.options.get_bool("FULL_HBAR_DEGNO") and self.solver_type in ["MRDSRG","SA-MRDSRG","SA_MRDSRG"] and n == self.relax_maxiter - 1: + psi4.core.print_out("\n =>** Saving Full Hbar in de-normal-ordered basis (relaxed) **<=\n") + Hbar1, Hbar2 = self.dsrg_solver.compute_Heff_full_degno() + np.savez("save_Hbar_degno", **Hbar1, **Hbar2) ints_dressed = self.dsrg_solver.compute_Heff_actv() if self.fno_pt2_Heff_shift is not None: