Parallelize similarity transformation code (#434)

* Add more helper functions

* Add tests for new functions

* Cleanup code

* Improved naming

* Add function to create a SparseOperator from an ActiveSpaceIntegral object

* Merge branch 'main' into mcscf_df

* Simplify logic

* Fix compilation

forte/CMakeLists.txt

@@ -278,6 +278,7 @@ sparse_ci/sparse_hamiltonian.cc
 sparse_ci/sparse_initial_guess.cc
 sparse_ci/sparse_operator.cc
 sparse_ci/sparse_operator_sim_trans.cc
+sparse_ci/sparse_operator_hamiltonian.cc
 sparse_ci/sparse_state.cc
 sparse_ci/sq_operator_string.cc
 sparse_ci/sq_operator_string_ops.cc

forte/api/determinant_api.cc

@@ -94,6 +94,7 @@ void export_Determinant(py::module& m) {
                const std::vector<int>& bann,
                const std::vector<int>& bcre) { return gen_excitation(d, aann, acre, bann, bcre); },
             "Apply a generic excitation") // uses gen_excitation() defined in determinant.hpp
+        .def("spin_flip", &Determinant::spin_flip, "Get the spin-flip determinant")
         .def(
             "str", [](const Determinant& a, int n) { return str(a, n); },
             "n"_a = Determinant::norb(),

forte/api/sparse_operator_api.cc

@@ -33,8 +33,11 @@
 
 #include "helpers/string_algorithms.h"
 
+#include "integrals/active_space_integrals.h"
+
 #include "sparse_ci/sparse_operator.h"
 #include "sparse_ci/sq_operator_string_ops.h"
+#include "sparse_ci/sparse_operator_hamiltonian.h"
 
 namespace py = pybind11;
 using namespace pybind11::literals;
@@ -252,6 +255,9 @@ void export_SparseOperator(py::module& m) {
         },
         "list"_a, "Create a SparseOperator object from a list of Tuple[SQOperatorString, complex]");
 
+    m.def("sparse_operator_hamiltonian", &sparse_operator_hamiltonian,
+          "Create a SparseOperator object from an ActiveSpaceIntegrals object");
+
     m.def("new_product", [](const SparseOperator A, const SparseOperator B) {
         SparseOperator C;
         SQOperatorProductComputer computer;

forte/api/sparse_operator_list_api.cc

@@ -49,6 +49,13 @@ void export_SparseOperatorList(py::module& m) {
         .def("add", &SparseOperatorList::add_term_from_str, "str"_a,
              "coefficient"_a = sparse_scalar_t(1), "allow_reordering"_a = false)
         .def("to_operator", &SparseOperatorList::to_operator)
+        .def(
+            "remove",
+            [](SparseOperatorList& op, const std::string& s) {
+                const auto [sqop, _] = make_sq_operator_string(s, false);
+                op.remove(sqop);
+            },
+            "Remove a specific element from the vector space")
         .def("__len__", &SparseOperatorList::size)
         .def(
             "__iter__",
@@ -85,6 +92,7 @@ void export_SparseOperatorList(py::module& m) {
                      op[i] = values[i];
                  }
              })
+        .def("reverse", &SparseOperatorList::reverse, "Reverse the order of the operators")
         .def(
             "__call__",
             [](const SparseOperatorList& op, const size_t n) {
@@ -93,7 +101,7 @@ void export_SparseOperatorList(py::module& m) {
                 }
                 return op(n);
             },
-            "Get the nth element");
+            "Get the nth operator");
 
     m.def(
         "operator_list",

forte/api/sparse_state_api.cc

@@ -50,7 +50,7 @@ void export_SparseState(py::module& m) {
         .def("str", &SparseState::str)
         .def("size", &SparseState::size)
         .def("norm", &SparseState::norm, "p"_a = 2,
-            "Calculate the p-norm of the SparseState (default p = 2, p = -1 for infinity norm)")
+             "Calculate the p-norm of the SparseState (default p = 2, p = -1 for infinity norm)")
         .def("add", &SparseState::add)
         .def("__iadd__", &SparseState::operator+=, "Add a SparseState to this SparseState")
         .def("__isub__", &SparseState::operator-=, "Subtract a SparseState from this SparseState")
@@ -72,12 +72,19 @@ void export_SparseState(py::module& m) {
           "screen_thresh"_a = 1.0e-12);
 
     m.def("apply_number_projector", &apply_number_projector);
+
     m.def("get_projection", &get_projection);
     // there's already a function called spin2, overload the spin2 function
-    m.def("spin2", [](const SparseState& left_state, const SparseState& right_state) {
-        return spin2(left_state, right_state);
-    }, "Calculate the <left_state|S^2|right_state> expectation value");
+
+    m.def(
+        "spin2",
+        [](const SparseState& left_state, const SparseState& right_state) {
+            return spin2(left_state, right_state);
+        },
+        "Calculate the <left_state|S^2|right_state> expectation value");
+
     m.def("overlap", &overlap);
-    m.def("normalize", &normalize, "Returns a normalized version of the input SparseState"); 
+
+    m.def("normalize", &normalize, "Returns a normalized version of the input SparseState");
 }
 } // namespace forte

forte/api/sq_operator_string_api.cc

@@ -48,6 +48,7 @@ void export_SQOperatorString(py::module& m) {
         .def("str", &SQOperatorString::str, "Get the string representation of the operator string")
         .def("count", &SQOperatorString::count, "Get the number of operators")
         .def("adjoint", &SQOperatorString::adjoint, "Get the adjoint operator string")
+        .def("spin_flip", &SQOperatorString::spin_flip, "Get the spin-flipped operator string")
         .def("number_component", &SQOperatorString::number_component,
              "Get the number component of the operator string")
         .def("non_number_component", &SQOperatorString::non_number_component,

forte/helpers/math_structures.h

@@ -345,6 +345,11 @@ template <typename Derived, typename T, typename F> class VectorSpaceList {
     /// @return the number of elements in the vector
     size_t size() const { return elements_.size(); }
 
+    /// @brief Remove a specific element from the vector space
+    void remove(const T& e) {
+        std::erase_if(elements_, [&e](const auto& p) { return p.first == e; });
+    }
+
     /// @return an element of the vector
     const F& operator[](size_t n) const { return elements_[n].second; }
 
@@ -408,6 +413,14 @@ template <typename Derived, typename T, typename F> class VectorSpaceList {
         return result;
     }
 
+    /// @brief Return a reversed copy of the vector
+    Derived reverse() {
+        // avoid issues with const
+        Derived result = static_cast<Derived&>(*this);
+        std::reverse(result.elements_.begin(), result.elements_.end());
+        return result;
+    }
+
   private:
     // Using an unordered_map with a custom hash function
     container elements_;

forte/sparse_ci/determinant.hpp

@@ -537,7 +537,7 @@ template <size_t N> class DeterminantImpl : public BitArray<N> {
         }
     }
 
-    BitArray<nbits_half> get_bits(DetSpinType spin_type) {
+    BitArray<nbits_half> get_bits(DetSpinType spin_type) const {
         return (spin_type == DetSpinType::Alpha ? get_alfa_bits() : get_beta_bits());
     }
 
@@ -552,6 +552,15 @@ template <size_t N> class DeterminantImpl : public BitArray<N> {
         for (size_t n = nwords_half; n < nwords_; n++)
             words_[n] = u_int64_t(0);
     }
+
+    /// Swap the alpha and beta bits of a determinant
+    DeterminantImpl<N> spin_flip() const {
+        DeterminantImpl<N> d(*this);
+        for (size_t n = 0; n < nwords_half; n++) {
+            std::swap(d.words_[n], d.words_[n + nwords_half]);
+        }
+        return d;
+    }
 };
 
 // Functions

forte/sparse_ci/sparse.h

@@ -31,12 +31,14 @@
 #include <complex>
 
 namespace forte {
+
+// Define the scalar type used in the SparseOperator and SparseState objects
 using sparse_scalar_t = std::complex<double>;
 
-// // For double
-// double to_double(const double& input) { return input; }
+// For double
+inline double to_double(const double& input) { return input; }
 
 // For std::complex<double>
-template <typename T> double to_double(const T& input) { return std::real(input); }
+inline double to_double(const std::complex<double>& input) { return std::real(input); }
 
 } // namespace forte

forte/sparse_ci/sparse_hamiltonian.cc

@@ -29,6 +29,7 @@
 #include <cmath>
 
 #include "sparse_ci/sparse_hamiltonian.h"
+#include "sparse_ci/sparse_operator_hamiltonian.h"
 
 namespace forte {
 
@@ -342,43 +343,7 @@ SparseState SparseHamiltonian::compute_on_the_fly(const SparseState& state, doub
 std::map<std::string, double> SparseHamiltonian::timings() const { return timings_; }
 
 SparseOperator SparseHamiltonian::to_sparse_operator() const {
-    SparseOperator H;
-    size_t nmo = as_ints_->nmo();
-
-    H.add_term_from_str("[]", as_ints_->nuclear_repulsion_energy() + as_ints_->scalar_energy() +
-                                  as_ints_->frozen_core_energy());
-    for (size_t p = 0; p < nmo; p++) {
-        for (size_t q = 0; q < nmo; q++) {
-            H.add_term_from_str(std::format("{}a+ {}a-", p, q), as_ints_->oei_a(p, q));
-            H.add_term_from_str(std::format("{}b+ {}b-", p, q), as_ints_->oei_b(p, q));
-        }
-    }
-
-    for (size_t p = 0; p < nmo; p++) {
-        for (size_t q = p + 1; q < nmo; q++) {
-            for (size_t r = 0; r < nmo; r++) {
-                for (size_t s = r + 1; s < nmo; s++) {
-                    H.add_term_from_str(std::format("{}a+ {}a+ {}a- {}a-", p, q, s, r),
-                                        as_ints_->tei_aa(p, q, r, s));
-                    H.add_term_from_str(std::format("{}b+ {}b+ {}b- {}b-", p, q, s, r),
-                                        as_ints_->tei_bb(p, q, r, s));
-                }
-            }
-        }
-    }
-
-    for (size_t p = 0; p < nmo; p++) {
-        for (size_t q = 0; q < nmo; q++) {
-            for (size_t r = 0; r < nmo; r++) {
-                for (size_t s = 0; s < nmo; s++) {
-                    H.add_term_from_str(std::format("{}a+ {}b+ {}b- {}a-", p, q, s, r),
-                                        as_ints_->tei_ab(p, q, r, s));
-                }
-            }
-        }
-    }
-
-    return H;
+    return sparse_operator_hamiltonian(as_ints_);
 }
 
 } // namespace forte

forte/sparse_ci/sparse_operator_hamiltonian.cc

@@ -0,0 +1,78 @@
+/*
+ * @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-2024 by its authors (see COPYING, COPYING.LESSER, 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 Lesser 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 Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see http://www.gnu.org/licenses/.
+ *
+ * @END LICENSE
+ */
+
+#include <format>
+
+#include "sparse_ci/sparse_operator_hamiltonian.h"
+
+#include "integrals/active_space_integrals.h"
+
+namespace forte {
+
+SparseOperator sparse_operator_hamiltonian(std::shared_ptr<ActiveSpaceIntegrals> as_ints) {
+    SparseOperator H;
+    size_t nmo = as_ints->nmo();
+
+    H.add_term_from_str("[]", as_ints->nuclear_repulsion_energy() + as_ints->scalar_energy() +
+                                  as_ints->frozen_core_energy());
+
+    for (size_t p = 0; p < nmo; p++) {
+        for (size_t q = 0; q < nmo; q++) {
+            H.add_term_from_str(std::format("{}a+ {}a-", p, q), as_ints->oei_a(p, q));
+            H.add_term_from_str(std::format("{}b+ {}b-", p, q), as_ints->oei_b(p, q));
+        }
+    }
+
+    for (size_t p = 0; p < nmo; p++) {
+        for (size_t q = p + 1; q < nmo; q++) {
+            for (size_t r = 0; r < nmo; r++) {
+                for (size_t s = r + 1; s < nmo; s++) {
+                    H.add_term_from_str(std::format("{}a+ {}a+ {}a- {}a-", p, q, s, r),
+                                        as_ints->tei_aa(p, q, r, s));
+                    H.add_term_from_str(std::format("{}b+ {}b+ {}b- {}b-", p, q, s, r),
+                                        as_ints->tei_bb(p, q, r, s));
+                }
+            }
+        }
+    }
+
+    for (size_t p = 0; p < nmo; p++) {
+        for (size_t q = 0; q < nmo; q++) {
+            for (size_t r = 0; r < nmo; r++) {
+                for (size_t s = 0; s < nmo; s++) {
+                    H.add_term_from_str(std::format("{}a+ {}b+ {}b- {}a-", p, q, s, r),
+                                        as_ints->tei_ab(p, q, r, s));
+                }
+            }
+        }
+    }
+
+    return H;
+}
+
+} // namespace forte

forte/sparse_ci/sparse_operator_hamiltonian.h

@@ -0,0 +1,43 @@
+/*
+ * @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-2024 by its authors (see COPYING, COPYING.LESSER, 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 Lesser 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 Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see http://www.gnu.org/licenses/.
+ *
+ * @END LICENSE
+ */
+
+#pragma once
+
+#include <memory>
+
+#include "sparse_ci/sparse_operator.h"
+
+namespace forte {
+
+class ActiveSpaceIntegrals;
+
+/// @brief Generate the a SparseOperator representation of the Hamiltonian using integrals from an
+/// ActiveSpaceIntegrals object
+SparseOperator sparse_operator_hamiltonian(std::shared_ptr<ActiveSpaceIntegrals> as_ints);
+
+} // namespace forte