Remove SparseHamiltonian, add sign_mask to SparseOperator

fevangelista · fevangelista · commit ec3db2508c3e · 2025-03-27T16:04:03.000-04:00
diff --git a/forte/CMakeLists.txt b/forte/CMakeLists.txt
@@ -102,7 +102,6 @@ api/scf_info_api.cc
 api/sparse_ci_solver_api.cc
 api/sparse_exp_api.cc
 api/sparse_fact_exp_api.cc
-api/sparse_hamiltonian_api.cc
 api/sparse_operator_api.cc
 api/sparse_operator_list_api.cc
 api/sparse_operator_sim_trans_api.cc
@@ -276,7 +275,6 @@ sparse_ci/sorted_string_list.cc
 sparse_ci/sparse_ci_solver.cc
 sparse_ci/sparse_exp.cc
 sparse_ci/sparse_fact_exp.cc
-sparse_ci/sparse_hamiltonian.cc
 sparse_ci/sparse_initial_guess.cc
 sparse_ci/sparse_operator.cc
 sparse_ci/sparse_operator_sim_trans.cc
diff --git a/forte/api/forte_python_module.cc b/forte/api/forte_python_module.cc
@@ -215,7 +215,6 @@ PYBIND11_MODULE(_forte, m) {
     export_SQOperatorString(m);
     export_SparseExp(m);
     export_SparseFactExp(m);
-    export_SparseHamiltonian(m);
     export_SparseOperator(m);
     export_SparseOperatorList(m);
     export_SparseOperatorSimTrans(m);
diff --git a/forte/api/forte_python_module.h b/forte/api/forte_python_module.h
@@ -64,7 +64,6 @@ void export_CI_RDMS(py::module& m);
 void export_SQOperatorString(py::module& m);
 void export_SparseExp(py::module& m);
 void export_SparseFactExp(py::module& m);
-void export_SparseHamiltonian(py::module& m);
 void export_SparseOperator(py::module& m);
 void export_SparseOperatorList(py::module& m);
 void export_SparseOperatorSimTrans(py::module& m);
diff --git a/forte/api/sparse_hamiltonian_api.cc b/forte/api/sparse_hamiltonian_api.cc
@@ -1,55 +0,0 @@
-/*
- * @BEGIN LICENSE
- *
- * Forte: an open-source plugin to Psi4 (https://github.com/psi4/psi4)
- * that implements a variety of quantum chemistry methods for strongly
- * correlated electrons.
- *
- * Copyright (c) 2012-2025 by its authors (see LICENSE, AUTHORS).
- *
- * The copyrights for code used from other parties are included in
- * the corresponding files.
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program.  If not, see http://www.gnu.org/licenses/.
- *
- * @END LICENSE
- */
-
-#include <pybind11/pybind11.h>
-#include <pybind11/stl.h>
-
-#include "sparse_ci/sparse_hamiltonian.h"
-
-namespace py = pybind11;
-using namespace pybind11::literals;
-
-namespace forte {
-
-void export_SparseHamiltonian(py::module& m) {
-    py::class_<SparseHamiltonian>(m, "SparseHamiltonian",
-                                  "A class to represent a sparse Hamiltonian")
-        .def(py::init<std::shared_ptr<ActiveSpaceIntegrals>>())
-        .def("compute", &SparseHamiltonian::compute, "state"_a, "screen_thresh"_a = 1.0e-12,
-             "Compute the state H|state> using an algorithm that caches the elements of H")
-        .def("apply", &SparseHamiltonian::compute, "state"_a, "screen_thresh"_a = 1.0e-12,
-             "Compute the state H|state> using an algorithm that caches the elements of H")
-        .def(
-            "compute_on_the_fly", &SparseHamiltonian::compute_on_the_fly, "state"_a,
-            "screen_thresh"_a = 1.0e-12,
-            "Compute the state H|state> using an on-the-fly algorithm that has no memory footprint")
-        .def("timings", &SparseHamiltonian::timings)
-        .def("to_sparse_operator", &SparseHamiltonian::to_sparse_operator,
-             "Convert the Hamiltonian to a SparseOperator object");
-}
-} // namespace forte
diff --git a/forte/api/sparse_operator_api.cc b/forte/api/sparse_operator_api.cc
@@ -224,6 +224,12 @@ void export_SparseOperator(py::module& m) {
         .def(
             "__str__", [](const SparseOperator& op) { return join(op.str(), "\n"); },
             "Get a string representation of the operator")
+        .def(
+            "apply_to_state",
+            [](const SparseOperator& op, const SparseState& state, double screen_thresh) {
+                return apply_operator_lin(op, state, screen_thresh);
+            },
+            "state"_a, "screen_thresh"_a = 1.0e-12, "Apply the operator to a state")
         .def(
             "fact_trans_lin",
             [](SparseOperator& O, const SparseOperatorList& T, bool reverse, double screen_thresh) {
diff --git a/forte/api/sparse_operator_list_api.cc b/forte/api/sparse_operator_list_api.cc
@@ -179,8 +179,14 @@ void export_SparseOperatorList(py::module& m) {
                 op1 += op2;
                 return op1;
             },
-            "Add (concatenate) a SparseOperatorList object to this SparseOperatorList object");
-        
+            "Add (concatenate) a SparseOperatorList object to this SparseOperatorList object")
+        .def(
+            "apply_to_state",
+            [](const SparseOperatorList& op, const SparseState& state, double screen_thresh) {
+                auto sop = op.to_operator();
+                return apply_operator_lin(sop, state, screen_thresh);
+            },
+            "state"_a, "screen_thresh"_a = 1.0e-12, "Apply the operator to a state");
 
     // Wrapper class that holds a SparseOperator
     // and overloads operator* to apply forte::SparseExp.
diff --git a/forte/attic/sparse_hamiltonian.cc b/forte/attic/sparse_hamiltonian.cc
diff --git a/forte/attic/sparse_hamiltonian.h b/forte/attic/sparse_hamiltonian.h
diff --git a/forte/modules/general_cc.py b/forte/modules/general_cc.py
@@ -16,10 +16,10 @@
     SparseOperatorList,
     Determinant,
     SparseState,
-    SparseHamiltonian,
     SparseFactExp,
     SparseExp,
     get_projection,
+    sparse_operator_hamiltonian,
 )
 
 
@@ -190,8 +190,8 @@ def solve_cc_equations(
         Returns the a tuple containign the converged amplitudes, the energy, the projective energy,
         the number of iterations, and timings information
     """
-    diis = DIIS(t, diis_start, diis_nvec)
-    ham = SparseHamiltonian(as_ints)
+    diis = DIIS(t, diis_start)
+    ham = sparse_operator_hamiltonian(as_ints)
     if cc_type == "cc" or cc_type == "ucc":
         exp = SparseExp(maxk=maxk, screen_thresh=compute_threshold)
     if cc_type == "dcc" or cc_type == "ducc":
@@ -255,19 +255,19 @@ def residual_equations(cc_type, t, op, sop, ref, ham, exp, compute_threshold, li
     c0 = 0.0
     if cc_type == "cc":
         wfn = exp.apply_op(sop, ref)
-        Hwfn = ham.compute_on_the_fly(wfn, compute_threshold)
+        Hwfn = ham.apply_to_state(wfn, compute_threshold)
         R = exp.apply_op(sop, Hwfn, scaling_factor=-1.0)
     elif cc_type == "ucc":
         wfn = exp.apply_antiherm(sop, ref)
-        Hwfn = ham.compute_on_the_fly(wfn, compute_threshold)
+        Hwfn = ham.apply_to_state(wfn, compute_threshold)
         R = exp.apply_antiherm(sop, Hwfn, scaling_factor=-1.0)
     elif cc_type == "dcc":
         wfn = exp.apply_op(sop, ref)
-        Hwfn = ham.compute_on_the_fly(wfn, compute_threshold)
+        Hwfn = ham.apply_to_state(wfn, compute_threshold)
         R = exp.apply_op(sop, Hwfn, inverse=True)
     elif cc_type == "ducc":
         wfn = exp.apply_antiherm(sop, ref)
-        Hwfn = ham.compute_on_the_fly(wfn, compute_threshold)
+        Hwfn = ham.apply_to_state(wfn, compute_threshold)
         R = exp.apply_antiherm(sop, Hwfn, inverse=True)
     else:
         raise ValueError("Incorrect value for cc_type")
@@ -340,7 +340,7 @@ def update(self, t, t_old):
         self.e_diis.append(np.subtract(t, t_old))
 
         diis_dim = len(self.t_diis)
-        if (diis_dim >= self.diis_start):
+        if diis_dim >= self.diis_start:
             # consturct diis B matrix (following Crawford Group github tutorial)
             B = np.ones((diis_dim + 1, diis_dim + 1)) * -1.0
             bsol = np.zeros(diis_dim + 1)
diff --git a/forte/sparse_ci/sparse_state_functions.cc b/forte/sparse_ci/sparse_state_functions.cc
@@ -74,11 +74,10 @@ SparseState apply_operator_impl_naive(bool is_antihermitian, const SparseOperato
     Determinant sign_mask; // a temporary determinant to store the sign mask
     Determinant idx;       // a temporary determinant to store the index of the determinant
     for (const auto& [sqop, t] : sop) {
-        compute_sign_mask(sqop.cre(), sqop.ann(), sign_mask, idx);
         for (const auto& [det, c] : state) {
             if (det.faster_can_apply_operator(sqop.cre(), sqop.ann())) {
-                auto value =
-                    faster_apply_operator_to_det(det, new_det, sqop.cre(), sqop.ann(), sign_mask);
+                auto value = faster_apply_operator_to_det(det, new_det, sqop.cre(), sqop.ann(),
+                                                          sqop.sign_mask());
                 new_terms[new_det] += value * t * c;
             }
         }
@@ -89,11 +88,10 @@ SparseState apply_operator_impl_naive(bool is_antihermitian, const SparseOperato
     }
 
     for (const auto& [sqop, t] : sop) {
-        compute_sign_mask(sqop.ann(), sqop.cre(), sign_mask, idx);
         for (const auto& [det, c] : state) {
             if (det.faster_can_apply_operator(sqop.ann(), sqop.cre())) {
-                auto value =
-                    faster_apply_operator_to_det(det, new_det, sqop.ann(), sqop.cre(), sign_mask);
+                auto value = faster_apply_operator_to_det(det, new_det, sqop.ann(), sqop.cre(),
+                                                          sqop.sign_mask());
                 new_terms[new_det] -= value * t * c;
             }
         }
@@ -109,18 +107,18 @@ template <bool positive>
 void apply_operator_kernel(const auto& sop_groups, const auto& state_sorted,
                            const auto& screen_thresh, auto& new_terms) {
     Determinant new_det;
-    Determinant sign_mask;
+    // Determinant sign_mask;
     Determinant idx;
     for (const auto& [sqop_ann, sqop_group] : sop_groups) {
         for (const auto& [det, c] : state_sorted) {
             if (det.fast_a_and_b_equal_b(sqop_ann)) {
                 // loop over the creation operators in this group
-                for (const auto& [sqop_cre, t] : sqop_group) {
+                for (const auto& [sqop_cre, sqop_sign_mask, t] : sqop_group) {
                     if (det.fast_a_and_b_minus_c_eq_zero(sqop_cre, sqop_ann)) {
                         if (std::abs(c * t) > screen_thresh) {
-                            compute_sign_mask(sqop_cre, sqop_ann, sign_mask, idx);
-                            const auto value = faster_apply_operator_to_det(det, new_det, sqop_cre,
-                                                                            sqop_ann, sign_mask);
+                            // compute_sign_mask(sqop_cre, sqop_ann, sign_mask, idx);
+                            const auto value = faster_apply_operator_to_det(
+                                det, new_det, sqop_cre, sqop_ann, sqop_sign_mask);
                             if constexpr (positive) {
                                 new_terms[new_det] += value * t * c;
                             } else {
@@ -151,11 +149,12 @@ SparseState apply_operator_impl_grouped(bool is_antihermitian, const SparseOpera
               [](const auto& a, const auto& b) { return a.first < b.first; });
 
     // Group the operators by common annihilation strings
-    std::unordered_map<Determinant, std::vector<std::pair<Determinant, sparse_scalar_t>>,
+    std::unordered_map<Determinant,
+                       std::vector<std::tuple<Determinant, Determinant, sparse_scalar_t>>,
                        Determinant::Hash>
         sop_groups;
     for (const auto& [sqop, t] : sop.elements()) {
-        sop_groups[sqop.ann()].emplace_back(sqop.cre(), t);
+        sop_groups[sqop.ann()].emplace_back(sqop.cre(), sqop.sign_mask(), t);
     }
 
     // Call the kernel to apply the operator (adding the result)
@@ -169,7 +168,7 @@ SparseState apply_operator_impl_grouped(bool is_antihermitian, const SparseOpera
     // Here we swap the annihilation and creation operators for the antihermitian case
     sop_groups.clear();
     for (const auto& [sqop, t] : sop.elements()) {
-        sop_groups[sqop.cre()].emplace_back(sqop.ann(), t);
+        sop_groups[sqop.cre()].emplace_back(sqop.ann(), sqop.sign_mask(), t);
     }
 
     // Call the kernel to apply the operator (subtracting the result)
@@ -186,18 +185,15 @@ template <bool positive>
 void apply_operator_kernel_string(const auto& sop_groups, const auto& state_groups,
                                   const auto& screen_thresh, auto& new_terms) {
     Determinant new_det;
-    Determinant sign_mask;
-    Determinant idx;
     for (const auto& [sqop_ann_a, sqop_group] : sop_groups) {
         for (const auto& [det_a, state_group] : state_groups) {
             // can we annihilate the alfa string?
             if (det_a.fast_a_and_b_equal_b(sqop_ann_a)) {
                 // loop over the creation operators in this group
-                for (const auto& [sqop_ann, sqop_cre, t] : sqop_group) {
+                for (const auto& [sqop_ann, sqop_cre, sign_mask, t] : sqop_group) {
                     for (const auto& [det, c] : state_group) {
                         if (det.faster_can_apply_operator(sqop_cre, sqop_ann)) {
                             if (std::abs(c * t) > screen_thresh) {
-                                compute_sign_mask(sqop_cre, sqop_ann, sign_mask, idx);
                                 const auto value = faster_apply_operator_to_det(
                                     det, new_det, sqop_cre, sqop_ann, sign_mask);
                                 if constexpr (positive) {
@@ -232,11 +228,13 @@ SparseState apply_operator_impl_grouped_string(bool is_antihermitian, const Spar
     }
 
     // Group the operators by common alfa annihilation strings
-    std::unordered_map<String, std::vector<std::tuple<Determinant, Determinant, sparse_scalar_t>>,
-                       String::Hash>
+    std::unordered_map<
+        String, std::vector<std::tuple<Determinant, Determinant, Determinant, sparse_scalar_t>>,
+        String::Hash>
         sop_groups;
     for (const auto& [sqop, t] : sop.elements()) {
-        sop_groups[sqop.ann().get_alfa_bits()].emplace_back(sqop.ann(), sqop.cre(), t);
+        sop_groups[sqop.ann().get_alfa_bits()].emplace_back(sqop.ann(), sqop.cre(),
+                                                            sqop.sign_mask(), t);
     }
 
     // Call the kernel to apply the operator (adding the result)
@@ -250,7 +248,8 @@ SparseState apply_operator_impl_grouped_string(bool is_antihermitian, const Spar
     // Here we swap the annihilation and creation operators for the antihermitian case
     sop_groups.clear();
     for (const auto& [sqop, t] : sop.elements()) {
-        sop_groups[sqop.cre().get_alfa_bits()].emplace_back(sqop.cre(), sqop.ann(), t);
+        sop_groups[sqop.cre().get_alfa_bits()].emplace_back(sqop.cre(), sqop.ann(),
+                                                            sqop.sign_mask(), t);
     }
 
     // Call the kernel to apply the operator (subtracting the result)
diff --git a/forte/sparse_ci/sq_operator_string.cc b/forte/sparse_ci/sq_operator_string.cc
@@ -59,13 +59,18 @@ bool compare_ops(const std::tuple<bool, bool, int>& lhs, const std::tuple<bool,
 SQOperatorString::SQOperatorString() {}
 
 SQOperatorString::SQOperatorString(const Determinant& cre, const Determinant& ann)
-    : cre_(cre), ann_(ann) {}
+    : cre_(cre), ann_(ann) {
+    Determinant temp;
+    compute_sign_mask(cre_, ann_, sign_mask_, temp);
+}
 
 SQOperatorString::SQOperatorString(const std::vector<size_t>& acre, const std::vector<size_t>& bcre,
                                    const std::vector<size_t>& aann,
                                    const std::vector<size_t>& bann) {
     cre_ = Determinant(acre, bcre);
     ann_ = Determinant(aann, bann);
+    Determinant temp;
+    compute_sign_mask(cre_, ann_, sign_mask_, temp);
 }
 
 SQOperatorString::SQOperatorString(const std::initializer_list<size_t> acre,
@@ -74,12 +79,16 @@ SQOperatorString::SQOperatorString(const std::initializer_list<size_t> acre,
                                    const std::initializer_list<size_t> bann) {
     cre_ = Determinant(acre, bcre);
     ann_ = Determinant(aann, bann);
+    Determinant temp;
+    compute_sign_mask(cre_, ann_, sign_mask_, temp);
 }
 
 const Determinant& SQOperatorString::cre() const { return cre_; }
 
 const Determinant& SQOperatorString::ann() const { return ann_; }
 
+const Determinant& SQOperatorString::sign_mask() const { return sign_mask_; }
+
 Determinant& SQOperatorString::cre_mod() { return cre_; }
 
 Determinant& SQOperatorString::ann_mod() { return ann_; }
diff --git a/forte/sparse_ci/sq_operator_string.h b/forte/sparse_ci/sq_operator_string.h
@@ -102,6 +102,8 @@ class SQOperatorString {
     const Determinant& cre() const;
     /// @return a Determinant object that represents the annihilation operators
     const Determinant& ann() const;
+    /// @return the sign mask associated with this operator
+    const Determinant& sign_mask() const;
     /// @return a Determinant object that represents the creation operators
     Determinant& cre_mod();
     /// @return a Determinant object that represents the annihilation operators
@@ -161,6 +163,8 @@ class SQOperatorString {
     Determinant cre_;
     /// a Determinant that represents the annihilation operators
     Determinant ann_;
+    /// a Determinant that represents the adjoint of the creation operators
+    Determinant sign_mask_;
 };
 
 class SQOperatorProductComputer {
diff --git a/tests/performance/sparse/sparse_operator_timing.py b/tests/performance/sparse/sparse_operator_timing.py
@@ -39,7 +39,6 @@ def sparse_operator_correctness():
     dets = forte.hilbert_space(wfn.nmo(), na, nb, nirrep, mo_symmetry, wfn_symmetry)
 
     # Build the Hamiltonian
-    spH = forte.SparseHamiltonian(as_ints)
     opH = forte.sparse_operator_hamiltonian(as_ints)
 
     print("\n\n  Number of determinants: ", len(dets))
@@ -49,10 +48,6 @@ def sparse_operator_correctness():
     c = 1 / np.sqrt(len(dets))
     state = forte.SparseState({det: c for det in dets})
     opH_state = forte.apply_op(opH, state)
-    spH_state = spH.apply(state)
-    spH_state -= opH_state
-    assert spH_state.norm() < 1e-9
-    print(f"  Norm of difference: {spH_state.norm()}")
 
 
 def sparse_operator_timing_1():
diff --git a/tests/pytest/sparse_ci/test_sparse_operator2.py b/tests/pytest/sparse_ci/test_sparse_operator2.py
@@ -22,16 +22,11 @@ def test_sparse_operator2():
 
     as_ints = data.as_ints  # forte_objs["as_ints"]
 
-    ham = forte.SparseHamiltonian(as_ints)
+    ham = forte.sparse_operator_hamiltonian(as_ints)
 
     ref = forte.SparseState({det("20"): 1.0})
-    Href1 = ham.compute(ref, 0.0)
-    Href2 = ham.compute_on_the_fly(ref, 0.0)
+    Href1 = ham @ ref
     assert Href1[det("20")] == pytest.approx(-1.094572, abs=1e-6)
-    assert Href2[det("20")] == pytest.approx(-1.094572, abs=1e-6)
-
-    ham_op = ham.to_sparse_operator()
-    assert forte.overlap(ref, forte.apply_op(ham_op, ref)) == pytest.approx(-1.094572, abs=1e-6)
 
     psi4.core.clean()
 
