fix spin-orbital rdms and prepare_forte_object

forte/api/forte_python_module.cc

@@ -367,9 +367,9 @@ PYBIND11_MODULE(_forte, m) {
         .def("compute_energy", &MRDSRG_SO::compute_energy, "Compute DSRG energy")
         .def("compute_Heff_actv", &MRDSRG_SO::compute_Heff_actv,
              "Return the DSRG dressed ActiveSpaceIntegrals")
-        .def("compute_Heff_full", [](MRDSRG_SO& self) {
-            const auto Heff = self.compute_Heff_full();
-            return py::make_tuple(blockedtensor_to_np(Heff.at(0)), blockedtensor_to_np(Heff.at(1)));
+        .def("save_Heff_full", [](MRDSRG_SO& self) {
+            const auto Heff = self.save_Heff_full();
+            return py::make_tuple(Heff.first, blockedtensor_to_np(Heff.second.at(0)), blockedtensor_to_np(Heff.second.at(1)));
             })
         .def("get_gamma1", [](MRDSRG_SO& self) {
             const auto gamma1 = self.get_gamma1();

forte/mrdsrg-so/mrdsrg_so.cc

@@ -646,6 +646,10 @@ double MRDSRG_SO::compute_energy() {
         ++cycle;
     } while (!converged);
 
+    if (foptions_->get_bool("FULL_HBAR")) {
+        compute_hbar();
+    }
+
     outfile->Printf("\n    "
                     "----------------------------------------------------------"
                     "----------------------------------------");
@@ -1372,9 +1376,13 @@ void MRDSRG_SO::H3_T2_C2(BlockedTensor& H3, BlockedTensor& T2, const double& alp
     temp["mx"] += 0.5 * T2["myuv"] * Lambda2["uvxy"];
     C2["toqr"] -= alpha * temp["mx"] * H3["xtomqr"];
 }
-std::vector<ambit::BlockedTensor> MRDSRG_SO::compute_Heff_full() {
-    compute_hbar();
-    std::vector<ambit::BlockedTensor> Heff = {Hbar1, Hbar2};
-    return Heff;
+std::pair<double, std::vector<BlockedTensor>> MRDSRG_SO::save_Heff_full() {
+    double Edsrg = Eref + Hbar0;
+    ambit::BlockedTensor Hbar1_copy = BTF_->build(tensor_type_, "Hbar1_copy", {"gg"});
+    ambit::BlockedTensor Hbar2_copy = BTF_->build(tensor_type_, "Hbar2_copy", {"gggg"});
+    Hbar1_copy["pq"] = Hbar1["pq"];
+    Hbar2_copy["pqrs"] = Hbar2["pqrs"];
+    std::vector<ambit::BlockedTensor> Heff = {Hbar1_copy, Hbar2_copy};
+    return std::make_pair(Edsrg, Heff);
 }
 } // namespace forte

forte/mrdsrg-so/mrdsrg_so.h

@@ -283,7 +283,7 @@ class MRDSRG_SO : public DynamicCorrelationSolver {
 
     // double compute_eom();
 
-    std::vector<ambit::BlockedTensor> compute_Heff_full();
+    std::pair<double, std::vector<BlockedTensor>> save_Heff_full();
     ambit::BlockedTensor get_gamma1() { return Gamma1; }
     ambit::BlockedTensor get_eta1() { return Eta1; }
     ambit::BlockedTensor get_lambda2() { return Lambda2; }

forte/proc/dsrg.py

@@ -285,8 +285,10 @@ 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 self.solver_type in ["MRDSRG","SA-MRDSRG","SA_MRDSRG","MRDSRG_SO","MRDSRG-SO"] \
-                and n == self.relax_maxiter - 1:
+            if self.options.get_bool("FULL_HBAR") and self.solver_type in ["SA-MRDSRG","SA_MRDSRG","MRDSRG_SO","MRDSRG-SO"]:
+                raise NotImplementedError("Relaxed full Hbar not implemented for SO/SA-DSRG")
+
+            if self.options.get_bool("FULL_HBAR") and self.solver_type in ["MRDSRG"] and n == self.relax_maxiter - 1:
                 psi4.core.print_out("\n  =>** Temporarily saving Full Hbar in de-normal-ordered basis (relaxed) **<=\n")
                 Hbar0, Hbar1, Hbar2 = self.dsrg_solver.save_Heff_full_ambit()
                 psi4.core.print_out(f"\n  The Hbar0 term is: {Hbar0}\n")

forte/utils/helpers.py

@@ -184,7 +184,7 @@ def psi4_cubeprop(wfn, path=".", orbs=[], nocc=0, nvir=0, density=False, frontie
 
 
 def prepare_forte_objects(
-    wfn, mo_spaces=None, active_space="ACTIVE", core_spaces=["RESTRICTED_DOCC"], localize=False, localize_spaces=[]
+    wfn, mo_spaces, active_space="ACTIVE", core_spaces=["RESTRICTED_DOCC"], localize=False, localize_spaces=[]
 ):
     """Take a psi4 wavefunction object and prepare the ForteIntegrals, SCFInfo, and MOSpaceInfo objects
 
@@ -235,15 +235,18 @@ def prepare_forte_objects(
     point_group = wfn.molecule().point_group().symbol()
 
     # create a MOSpaceInfo object
-    if mo_spaces is None:
-        mo_space_info = forte.make_mo_space_info(nmopi, point_group, options)
-    else:
-        mo_space_info = forte.make_mo_space_info_from_map(nmopi, point_group, mo_spaces)
+    mo_space_info = forte.make_mo_space_info_from_map(nmopi, point_group, mo_spaces)
+
+    # These variables are needed in make_state_weights_map
+    nel = wfn.nalpha() + wfn.nbeta()
+    multiplicity = 1
 
+    options.set_int("NEL", nel)
+    options.set_int("MULTIPLICITY", multiplicity)
     state_weights_map = forte.make_state_weights_map(options, mo_space_info)
 
     # make a ForteIntegral object
-    ints = forte.make_ints_from_psi4(wfn, options, mo_space_info)
+    ints = forte.make_ints_from_psi4(wfn, options, scf_info, mo_space_info)
 
     if localize:
         localizer = forte.Localize(forte.forte_options, ints, mo_space_info)