Merge branch 'sparse_parallel' into eom-dsrg

brianz98 · brianz98 · commit 56f917565fc0 · 2025-02-24T10:54:09.000-05:00
diff --git a/forte/api/sparse_exp_api.cc b/forte/api/sparse_exp_api.cc
@@ -39,6 +39,8 @@ namespace forte {
 void export_SparseExp(py::module& m) {
     py::class_<SparseExp>(m, "SparseExp", "A class to compute the exponential of a sparse operator")
         .def(py::init<int, double>(), "maxk"_a = 19, "screen_thresh"_a = 1.0e-12)
+        .def("set_screen_thresh", &SparseExp::set_screen_thresh, "screen_thresh"_a)
+        .def("set_maxk", &SparseExp::set_maxk, "maxk"_a)
         .def("apply_op",
              py::overload_cast<const SparseOperator&, const SparseState&, double>(
                  &SparseExp::apply_op),
diff --git a/forte/api/sparse_fact_exp_api.cc b/forte/api/sparse_fact_exp_api.cc
@@ -41,9 +41,10 @@ void export_SparseFactExp(py::module& m) {
         m, "SparseFactExp",
         "A class to compute the product exponential of a sparse operator using factorization")
         .def(py::init<double>(), "screen_thresh"_a = 1.0e-12)
-        .def("apply_op", &SparseFactExp::apply_op, "sop"_a, "state"_a, "inverse"_a = false)
+        .def("set_screen_thresh", &SparseFactExp::set_screen_thresh, "screen_thresh"_a)
+        .def("apply_op", &SparseFactExp::apply_op, "sop"_a, "state"_a, "inverse"_a = false, "reverse"_a = false)
         .def("apply_antiherm", &SparseFactExp::apply_antiherm, "sop"_a, "state"_a,
-             "inverse"_a = false);
+             "inverse"_a = false, "reverse"_a = false);
 }
 
 } // namespace forte
diff --git a/forte/sparse_ci/sparse_exp.h b/forte/sparse_ci/sparse_exp.h
@@ -97,6 +97,10 @@ class SparseExp {
     SparseState apply_antiherm(const SparseOperatorList& sop, const SparseState& state,
                                double scaling_factor = 1.0);
 
+    void set_screen_thresh(double screen_thresh) { screen_thresh_ = screen_thresh; }
+
+    void set_maxk(int maxk) { maxk_ = maxk; }
+
   private:
     int maxk_ = 19;
     double screen_thresh_ = 1e-12;
diff --git a/forte/sparse_ci/sparse_fact_exp.cc b/forte/sparse_ci/sparse_fact_exp.cc
@@ -36,7 +36,7 @@ namespace forte {
 SparseFactExp::SparseFactExp(double screen_thresh) : screen_thresh_(screen_thresh) {}
 
 SparseState SparseFactExp::apply_op(const SparseOperatorList& sop, const SparseState& state,
-                                    bool inverse) {
+                                    bool inverse, bool reverse) {
     // initialize a state object
     SparseState result(state);
 
@@ -47,15 +47,10 @@ SparseState SparseFactExp::apply_op(const SparseOperatorList& sop, const SparseS
     Determinant sign_mask;
     Determinant idx;
     for (size_t m = 0, nterms = sop.size(); m < nterms; m++) {
-        size_t n = inverse ? nterms - m - 1 : m;
+        size_t n = (inverse ^ reverse) ? nterms - m - 1 : m;
         const auto& [sqop, coefficient] = sop(n);
-        if (not sqop.is_nilpotent()) {
-            std::string msg =
-                "compute_on_the_fly_excitation is implemented only for nilpotent operators."
-                "Operator " +
-                sqop.str() + " is not nilpotent";
-            throw std::runtime_error(msg);
-        }
+        bool is_idempotent = !sqop.is_nilpotent();
+
         compute_sign_mask(sqop.cre(), sqop.ann(), sign_mask, idx);
         const auto t = (inverse ? -1.0 : 1.0) * coefficient;
         const auto screen_thresh_div_t = screen_thresh_ / std::abs(t);
@@ -66,9 +61,13 @@ SparseState SparseFactExp::apply_op(const SparseOperatorList& sop, const SparseS
             // test if we can apply this operator to this determinant
             if ((std::abs(c) > screen_thresh_div_t) and
                 det.faster_can_apply_operator(sqop.cre(), sqop.ann())) {
-                const auto sign =
-                    faster_apply_operator_to_det(det, new_det, sqop.cre(), sqop.ann(), sign_mask);
-                new_terms.push_back(std::make_pair(new_det, c * t * sign));
+                if (is_idempotent) {
+                    new_terms.emplace_back(det, c * (std::exp(t) - 1.0));
+                } else {
+                    const auto sign = faster_apply_operator_to_det(det, new_det, sqop.cre(),
+                                                                   sqop.ann(), sign_mask);
+                    new_terms.emplace_back(new_det, c * t * sign);
+                }
             }
         }
         for (const auto& [det, c] : new_terms) {
@@ -82,7 +81,7 @@ SparseState SparseFactExp::apply_op(const SparseOperatorList& sop, const SparseS
 }
 
 SparseState SparseFactExp::apply_antiherm(const SparseOperatorList& sop, const SparseState& state,
-                                          bool inverse) {
+                                          bool inverse, bool reverse) {
 
     // initialize a state object
     SparseState result(state);
@@ -92,16 +91,11 @@ SparseState SparseFactExp::apply_antiherm(const SparseOperatorList& sop, const S
     Determinant sign_mask;
     Determinant idx;
     for (size_t m = 0, nterms = sop.size(); m < nterms; m++) {
-        size_t n = inverse ? nterms - m - 1 : m;
+        size_t n = (inverse ^ reverse) ? nterms - m - 1 : m;
 
         const auto& [sqop, coefficient] = sop(n);
-        if (not sqop.is_nilpotent()) {
-            std::string msg =
-                "compute_on_the_fly_antihermitian is implemented only for nilpotent operators."
-                "Operator " +
-                sqop.str() + " is not nilpotent";
-            throw std::runtime_error(msg);
-        }
+        bool is_idempotent = !sqop.is_nilpotent();
+
         compute_sign_mask(sqop.cre(), sqop.ann(), sign_mask, idx);
         const auto t = (inverse ? -1.0 : 1.0) * coefficient;
         const auto screen_thresh_div_t = screen_thresh_ / std::abs(t);
@@ -111,19 +105,23 @@ SparseState SparseFactExp::apply_antiherm(const SparseOperatorList& sop, const S
             // to have an amplitude less than screen_thresh
             // (here we use the approximation sin(x) ~ x, for x small)
             if (std::abs(c) > screen_thresh_div_t) {
-                if (det.faster_can_apply_operator(sqop.cre(), sqop.ann())) {
-                    const auto theta = t * faster_apply_operator_to_det(det, new_det, sqop.cre(),
-                                                                        sqop.ann(), sign_mask);
-                    new_terms.emplace_back(det, c * (std::cos(std::abs(theta)) - 1.0));
-                    new_terms.emplace_back(new_det, c * std::polar(1.0, std::arg(theta)) *
-                                                        std::sin(std::abs(theta)));
-                } else if (det.faster_can_apply_operator(sqop.ann(), sqop.cre())) {
-                    const auto theta =
-                        -std::conj(t) * faster_apply_operator_to_det(det, new_det, sqop.ann(),
-                                                                     sqop.cre(), sign_mask);
-                    new_terms.emplace_back(det, c * (std::cos(std::abs(theta)) - 1.0));
-                    new_terms.emplace_back(new_det, c * std::polar(1.0, std::arg(theta)) *
-                                                        std::sin(std::abs(theta)));
+                if (is_idempotent and det.faster_can_apply_operator(sqop.cre(), sqop.ann())) {
+                    new_terms.emplace_back(det, c * (std::polar(1.0, 2.0 * std::imag(t)) - 1.0));
+                } else {
+                    if (det.faster_can_apply_operator(sqop.cre(), sqop.ann())) {
+                        const auto theta = t * faster_apply_operator_to_det(
+                                                   det, new_det, sqop.cre(), sqop.ann(), sign_mask);
+                        new_terms.emplace_back(det, c * (std::cos(std::abs(theta)) - 1.0));
+                        new_terms.emplace_back(new_det, c * std::polar(1.0, std::arg(theta)) *
+                                                            std::sin(std::abs(theta)));
+                    } else if (det.faster_can_apply_operator(sqop.ann(), sqop.cre())) {
+                        const auto theta =
+                            -std::conj(t) * faster_apply_operator_to_det(det, new_det, sqop.ann(),
+                                                                         sqop.cre(), sign_mask);
+                        new_terms.emplace_back(det, c * (std::cos(std::abs(theta)) - 1.0));
+                        new_terms.emplace_back(new_det, c * std::polar(1.0, std::arg(theta)) *
+                                                            std::sin(std::abs(theta)));
+                    }
                 }
             }
         }
diff --git a/forte/sparse_ci/sparse_fact_exp.h b/forte/sparse_ci/sparse_fact_exp.h
@@ -64,7 +64,10 @@ class SparseFactExp {
     ///
     ///             exp(-op1) exp(-op2) ... |state>
     ///
-    SparseState apply_op(const SparseOperatorList& sop, const SparseState& state, bool inverse);
+    /// @param reverse If true, apply the operator in reverse order:
+    ///             ... exp(opN-1) exp(opN) |state>
+    SparseState apply_op(const SparseOperatorList& sop, const SparseState& state, 
+                         bool inverse, bool reverse = false);
 
     /// @brief Compute the factorized exponential applied to a state using an exact algorithm
     ///
@@ -82,8 +85,12 @@ class SparseFactExp {
     ///
     ///             exp(-op1 + op1^dagger) exp(-op2 + op2^dagger) ... |state>
     ///
+    /// @param reverse If true, apply the operator in reverse order:
+    ///             ... exp(opN-1) exp(opN) |state>
     SparseState apply_antiherm(const SparseOperatorList& sop, const SparseState& state,
-                               bool inverse);
+                               bool inverse, bool reverse = false);
+
+    void set_screen_thresh(double screen_thresh) { screen_thresh_ = screen_thresh; }
 
   private:
     double screen_thresh_;
diff --git a/tests/pytest/sparse_ci/test_sparse_exp.py b/tests/pytest/sparse_ci/test_sparse_exp.py
@@ -3,11 +3,10 @@
 import numpy as np
 import time
 from forte import det
+import pytest
 
 
 def test_sparse_exp_1():
-    import pytest
-
     ### Test the linear operator ###
     op = forte.SparseOperator()
     ref = forte.SparseState({det("22"): 1.0})
@@ -26,6 +25,8 @@ def test_sparse_exp_1():
     assert wfn[det("0220")] == pytest.approx(0.15, abs=1e-9)
     assert wfn[det("2002")] == pytest.approx(-0.21, abs=1e-9)
 
+
+def test_sparse_exp_2():
     ### Test the exponential operator with excitation operator ###
     op = forte.SparseOperator()
     ref = forte.SparseState({det("22"): 1.0})
@@ -53,6 +54,8 @@ def test_sparse_exp_1():
     assert wfn[det("+0-2")] == pytest.approx(-0.0077, abs=1e-9)
     assert wfn[det("-0+2")] == pytest.approx(-0.0077, abs=1e-9)
 
+
+def test_sparse_exp_3():
     ### Test the exponential operator with antihermitian operator ###
     op = forte.SparseOperator()
     ref = forte.SparseState({det("22"): 1.0})
@@ -75,6 +78,8 @@ def test_sparse_exp_1():
     assert wfn2[det("+2-0")] == pytest.approx(0.0, abs=1e-9)
     assert wfn2[det("-2+0")] == pytest.approx(0.0, abs=1e-9)
 
+
+def test_sparse_exp_4():
     ### Test the factorized exponential operator with an antihermitian operator ###
     op = forte.SparseOperatorList()
     op.add("[2a+ 0a-]", 0.1)
@@ -112,7 +117,6 @@ def test_sparse_exp_1():
     op.add("[1b+ 0b-]", 0.05)
     op.add("[2a+ 2b+ 1b- 1a-]", -0.07)
 
-    dtest = det("20")
     ref = forte.SparseState({det("20"): 0.5, det("02"): 0.8660254038})
     factexp = forte.SparseFactExp()
     wfn = factexp.apply_antiherm(op, ref)
@@ -123,56 +127,114 @@ def test_sparse_exp_1():
     assert wfn[det("020")] == pytest.approx(0.676180171388, abs=1e-9)
     assert wfn[det("-+0")] == pytest.approx(0.016058887563, abs=1e-9)
 
-    ### Test the factorized exponential operator with an antihermitian operator with complex coefficients ###
+    ### Test idempotent operators with complex coefficients ###
     op = forte.SparseOperatorList()
-    op.add("[1a+ 0a-]", 0.1 + 0.2j)
+    op.add("[0a+ 0a-]", np.pi * 0.25j)
+    exp = forte.SparseExp(maxk=100, screen_thresh=1e-15)
+    factexp = forte.SparseFactExp()
+    ref = forte.SparseState({forte.det("20"): 1.0})
+    s1 = exp.apply_op(op, ref)
+    s2 = factexp.apply_op(op, ref)
+    assert s1[det("20")] == pytest.approx(s2[det("20")], abs=1e-9)
+    assert s2[det("20")] == pytest.approx(np.sqrt(2) * (1.0 + 1.0j) / 2, abs=1e-9)
+    s1 = exp.apply_antiherm(op, ref)
+    s2 = factexp.apply_antiherm(op, ref)
+    assert s1[det("20")] == pytest.approx(s2[det("20")], abs=1e-9)
+    assert s2[det("20")] == pytest.approx(1.0j, abs=1e-9)
+    op = forte.SparseOperatorList()
+    op.add("[1a+ 1a-]", np.pi * 0.25j)
+    s1 = exp.apply_antiherm(op, ref)
+    s2 = factexp.apply_antiherm(op, ref)
+    assert s1[det("20")] == pytest.approx(s2[det("20")], abs=1e-9)
+    assert s2[det("20")] == pytest.approx(1.0, abs=1e-9)
 
-    op_inv = forte.SparseOperatorList()
-    op_inv.add("[0a+ 1a-]", 0.1 - 0.2j)
 
-    exp = forte.SparseFactExp()
-    ref = forte.SparseState({forte.det("20"): 0.5, forte.det("02"): 0.8660254038})
+def test_sparse_exp_5():
+    ### Test the reverse argument of factorized exponential operator ###
 
-    s1 = exp.apply_antiherm(op, ref)
-    s2 = exp.apply_antiherm(op_inv, s1)
-    assert s2[det("20")] == pytest.approx(0.5, abs=1e-9)
-    assert s2[det("02")] == pytest.approx(0.8660254038, abs=1e-9)
+    # this is the manually reversed op from test_sparse_exp_4
+    op = forte.SparseOperatorList()
+    op.add("[2a+ 2b+ 0b- 0a-]", 0.15)
+    op.add("[2b+ 0b-]", 0.2)
+    op.add("[2a+ 0a-]", 0.1)
+    ref = forte.SparseState({det("22"): 1.0})
 
-    s1 = exp.apply_antiherm(op, ref, inverse=True)
-    s2 = exp.apply_antiherm(op_inv, ref, inverse=False)
-    assert s1 == s2
+    factexp = forte.SparseFactExp()
+    wfn = factexp.apply_antiherm(op, ref, reverse=True)
+
+    assert wfn[det("+2-0")] == pytest.approx(-0.197676811654, abs=1e-9)
+    assert wfn[det("-2+0")] == pytest.approx(-0.097843395007, abs=1e-9)
+    assert wfn[det("0220")] == pytest.approx(+0.165338757995, abs=1e-9)
+    assert wfn[det("2200")] == pytest.approx(+0.961256283877, abs=1e-9)
+
+    wfn2 = factexp.apply_antiherm(op, wfn, inverse=True, reverse=True)
+
+    assert wfn2[det("2200")] == pytest.approx(1.0, abs=1e-9)
 
-    ### Test the exponential operator with an antihermitian operator with complex coefficients ###
+    op = forte.SparseOperatorList()
+    op.add("[2a+ 2b+ 1b- 1a-]", -0.07)
+    op.add("[1b+ 0b-]", 0.05)
+    op.add("[1a+ 1b+ 0b- 0a-]", -0.3)
+    op.add("[1a+ 0a-]", 0.1)
+
+    ref = forte.SparseState({det("20"): 0.5, det("02"): 0.8660254038})
+    factexp = forte.SparseFactExp()
+    wfn = factexp.apply_antiherm(op, ref, reverse=True)
+
+    assert wfn[det("200")] == pytest.approx(0.733340213919, abs=1e-9)
+    assert wfn[det("+-0")] == pytest.approx(-0.049868863373, abs=1e-9)
+    assert wfn[det("002")] == pytest.approx(-0.047410073759, abs=1e-9)
+    assert wfn[det("020")] == pytest.approx(0.676180171388, abs=1e-9)
+    assert wfn[det("-+0")] == pytest.approx(0.016058887563, abs=1e-9)
+
+    op = forte.SparseOperatorList()
+    ref = forte.SparseState({det("22"): 1.0})
+    op.add("[2a+ 0a-]", 0.1)
+    op.add("[2b+ 0b-]", 0.1)
+    op.add("[2a+ 2b+ 0b- 0a-]", 0.15)
+    op.add("[3a+ 3b+ 1b- 1a-]", -0.077)
+
+    exp = forte.SparseFactExp()
+    wfn = exp.apply_op(op, ref, reverse=True)
+    assert wfn[det("2200")] == pytest.approx(1.0, abs=1e-9)
+    assert wfn[det("0220")] == pytest.approx(0.16, abs=1e-9)
+    assert wfn[det("+2-0")] == pytest.approx(-0.1, abs=1e-9)
+    assert wfn[det("-2+0")] == pytest.approx(-0.1, abs=1e-9)
+    assert wfn[det("2002")] == pytest.approx(-0.077, abs=1e-9)
+    assert wfn[det("+0-2")] == pytest.approx(-0.0077, abs=1e-9)
+    assert wfn[det("-0+2")] == pytest.approx(-0.0077, abs=1e-9)
+
+
+def test_sparse_exp_6():
+    ### Test the factorized exponential operator with an antihermitian operator with complex coefficients ###
     op = forte.SparseOperatorList()
     op.add("[1a+ 0a-]", 0.1 + 0.2j)
-    op_explicit = forte.SparseOperatorList()
-    op_explicit.add("[1a+ 0a-]", 0.1 + 0.2j)
-    op_explicit.add("[0a+ 1a-]", -0.1 + 0.2j)
 
     op_inv = forte.SparseOperatorList()
     op_inv.add("[0a+ 1a-]", 0.1 - 0.2j)
-    op_inv_explicit = forte.SparseOperatorList()
-    op_inv_explicit.add("[0a+ 1a-]", 0.1 - 0.2j)
-    op_inv_explicit.add("[1a+ 0a-]", -0.1 - 0.2j)
 
-    exp = forte.SparseExp()
+    exp = forte.SparseExp(maxk=100, screen_thresh=1e-15)
+    factexp = forte.SparseFactExp()
     ref = forte.SparseState({forte.det("20"): 0.5, forte.det("02"): 0.8660254038})
 
     s1 = exp.apply_antiherm(op, ref)
-    s1_explicit = exp.apply_op(op_explicit, ref)
-    assert s1 == s1_explicit
+    s2 = factexp.apply_antiherm(op, ref)
+    assert s1[det("20")] == pytest.approx(s2[det("20")], abs=1e-9)
+    assert s1[det("02")] == pytest.approx(s2[det("02")], abs=1e-9)
+    assert s1[det("+-")] == pytest.approx(s2[det("+-")], abs=1e-9)
+    assert s1[det("-+")] == pytest.approx(s2[det("-+")], abs=1e-9)
+
+    s1 = exp.apply_antiherm(op, ref)
     s2 = exp.apply_antiherm(op_inv, s1)
     assert s2[det("20")] == pytest.approx(0.5, abs=1e-9)
     assert s2[det("02")] == pytest.approx(0.8660254038, abs=1e-9)
-    s2_explicit = exp.apply_op(op_inv_explicit, s1)
-    assert s2 == s2_explicit
 
-    s1 = exp.apply_antiherm(op, ref)
-    s2 = exp.apply_antiherm(op_inv, ref, scaling_factor=-1.0)
+    s1 = factexp.apply_antiherm(op, ref, inverse=True)
+    s2 = factexp.apply_antiherm(op_inv, ref, inverse=False)
     assert s1 == s2
 
 
-def test_sparse_exp_2():
+def test_sparse_exp_7():
     # Compare the performance of the two methods to apply an operator to a state
     # when the operator all commute with each other
     norb = 10
@@ -239,3 +301,8 @@ def test_sparse_exp_2():
 if __name__ == "__main__":
     test_sparse_exp_1()
     test_sparse_exp_2()
+    test_sparse_exp_3()
+    test_sparse_exp_4()
+    test_sparse_exp_5()
+    test_sparse_exp_6()
+    test_sparse_exp_7()
