Apply op rewrite (#436)

* First pass at an implementation

* Add missing file

* Base implementation of linear operators

* Fsat version accomplished

* Use templates in the BitArray class to create fast implementations

* Better function documentation

* Final cleanup

forte/CMakeLists.txt

@@ -280,6 +280,7 @@ 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/sparse_state_functions.cc
 sparse_ci/sq_operator_string.cc
 sparse_ci/sq_operator_string_ops.cc
 v2rdm/v2rdm.cc

forte/api/determinant_api.cc

@@ -161,5 +161,11 @@ void export_Determinant(py::module& m) {
         .def("determinants", &DeterminantHashVec::determinants, "Return a vector of Determinants")
         .def("get_det", &DeterminantHashVec::get_det, "Return a specific determinant by reference")
         .def("get_idx", &DeterminantHashVec::get_idx, " Return the index of a determinant");
+
+    m.def(
+        "hilbert_space", &make_hilbert_space, "nmo"_a, "na"_a, "nb"_a, "nirrep"_a = 1,
+        "mo_symmetry"_a = std::vector<int>(), "symmetry"_a = 0,
+        "Generate the Hilbert space for a given number of electrons and orbitals. If information "
+        "about the symmetry of the MOs is not provided, it assumes that all MOs have symmetry 0.");
 }
 } // namespace forte

forte/api/sparse_hamiltonian_api.cc

@@ -40,10 +40,16 @@ 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)
-        .def("apply", &SparseHamiltonian::compute)
-        .def("compute_on_the_fly", &SparseHamiltonian::compute_on_the_fly)
+        .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);
+        .def("to_sparse_operator", &SparseHamiltonian::to_sparse_operator,
+             "Convert the Hamiltonian to a SparseOperator object");
 }
 } // namespace forte

forte/api/sparse_operator_api.cc

@@ -256,7 +256,8 @@ 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");
+          "Create a SparseOperator object from an ActiveSpaceIntegrals object", "as_ints"_a,
+          "screen_thresh"_a = 1.0e-12);
 
     m.def("new_product", [](const SparseOperator A, const SparseOperator B) {
         SparseOperator C;

forte/api/sparse_state_api.cc

@@ -52,6 +52,10 @@ void export_SparseState(py::module& m) {
         .def("norm", &SparseState::norm, "p"_a = 2,
              "Calculate the p-norm of the SparseState (default p = 2, p = -1 for infinity norm)")
         .def("add", &SparseState::add)
+        .def("__add__", &SparseState::operator+, "Add two SparseStates")
+        .def(
+            "__sub__", [](const SparseState& a, const SparseState& b) { return a - b; },
+            "Subtract two SparseStates")
         .def("__iadd__", &SparseState::operator+=, "Add a SparseState to this SparseState")
         .def("__isub__", &SparseState::operator-=, "Subtract a SparseState from this SparseState")
         .def("__imul__", &SparseState::operator*=, "Multiply this SparseState by a scalar")
@@ -74,8 +78,8 @@ void export_SparseState(py::module& m) {
     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
 
+    // there's already a function called spin2, overload the spin2 function
     m.def(
         "spin2",
         [](const SparseState& left_state, const SparseState& right_state) {

forte/forte.cc

@@ -38,7 +38,6 @@
 #include "base_classes/mo_space_info.h"
 #include "helpers/timer.h"
 
-#include "sparse_ci/determinant.h"
 #include "version.h"
 #include "psi4/libpsi4util/PsiOutStream.h"
 #include "psi4/psi4-dec.h"

forte/genci/string_lists_makers.h

@@ -28,14 +28,6 @@
 
 #pragma once
 
-// #include "psi4/libmints/dimension.h"
-
-// #include <map>
-// #include <vector>
-// #include <utility>
-
-// #include "helpers/timer.h"
-// #include "sparse_ci/determinant.h"
 #include "fci/string_list_defs.h"
 #include "genci_string_address.h"
 

forte/helpers/determinant_helpers.cc

@@ -27,6 +27,8 @@
  * @END LICENSE
  */
 
+#include <algorithm>
+
 #include "psi4/libmints/matrix.h"
 
 #include "integrals/active_space_integrals.h"
@@ -83,4 +85,69 @@ make_hamiltonian_matrix(const std::vector<Determinant>& dets,
     }
     return H;
 }
+
+std::vector<std::vector<String>> make_strings(int n, int k, size_t nirrep,
+                                              const std::vector<int>& mo_symmetry) {
+    // n is the number of orbitals
+    // k is the number of electrons
+    std::vector<std::vector<String>> strings(nirrep);
+    if ((k >= 0) and (k <= n)) { // check that (n > 0) makes sense.
+        String I;
+        const auto I_begin = I.begin();
+        const auto I_end = I.begin() + n;
+        // Generate the string 00000001111111
+        //                      {n-k}  { k }
+        I.zero();
+        for (int i = std::max(0, n - k); i < n; ++i)
+            I[i] = true; // 1
+        do {
+            int sym{0};
+            for (int i = 0; i < n; ++i) {
+                if (I[i])
+                    sym ^= mo_symmetry[i];
+            }
+            strings[sym].push_back(I);
+        } while (std::next_permutation(I_begin, I_end));
+    }
+    return strings;
+}
+
+std::vector<Determinant> make_hilbert_space(size_t nmo, size_t na, size_t nb, size_t nirrep,
+                                            std::vector<int> mo_symmetry, int symmetry) {
+    std::vector<Determinant> dets;
+    if (mo_symmetry.size() != nmo) {
+        mo_symmetry = std::vector<int>(nmo, 0);
+    }
+    // find the maximum value in mo_symmetry and check that it is less than nirrep
+    int max_sym = *std::max_element(mo_symmetry.begin(), mo_symmetry.end());
+    if (max_sym >= nirrep) {
+        throw std::runtime_error("The symmetry of the MOs is greater than the number of irreps.");
+    }
+    // implement other sensible checks, like making sure that symmetry is less than nirrep and na <=
+    // nmo, nb <= nmo
+    if (symmetry >= nirrep) {
+        throw std::runtime_error(
+            "The symmetry of the determinants is greater than the number of irreps.");
+    }
+    if (na > nmo) {
+        throw std::runtime_error(
+            "The number of alpha electrons is greater than the number of MOs.");
+    }
+    if (nb > nmo) {
+        throw std::runtime_error("The number of beta electrons is greater than the number of MOs.");
+    }
+
+    auto strings_a = make_strings(nmo, na, nirrep, mo_symmetry);
+    auto strings_b = make_strings(nmo, nb, nirrep, mo_symmetry);
+    for (size_t ha = 0; ha < nirrep; ha++) {
+        int hb = symmetry ^ ha;
+        for (const auto& Ia : strings_a[ha]) {
+            for (const auto& Ib : strings_b[hb]) {
+                dets.push_back(Determinant(Ia, Ib));
+            }
+        }
+    }
+    return dets;
+}
+
 } // namespace forte

forte/helpers/determinant_helpers.h

@@ -52,4 +52,25 @@ std::shared_ptr<psi::Matrix> make_s2_matrix(const std::vector<Determinant>& dets
 std::shared_ptr<psi::Matrix> make_hamiltonian_matrix(const std::vector<Determinant>& dets,
                                                      std::shared_ptr<ActiveSpaceIntegrals> as_ints);
 
+/// @brief Generate all strings of n orbitals and k electrons in each irrep
+/// @param n The number of orbitals
+/// @param k The number of electrons
+/// @param nirrep The number of irreps
+/// @param mo_symmetry The symmetry of the MOs
+/// @return A vector of vectors of strings, one vector for each irrep
+std::vector<std::vector<String>> make_strings(int n, int k, size_t nirrep,
+                                              const std::vector<int>& mo_symmetry);
+
+/// @brief Generate the Hilbert space for a given number of electrons and orbitals
+/// @param nmo The number of orbitals
+/// @param na The number of alpha electrons
+/// @param nb The number of beta electrons
+/// @param nirrep The number of irreps
+/// @param mo_symmetry The symmetry of the MOs
+/// @param symmetry The symmetry of the determinants
+/// @return A vector of determinants
+std::vector<Determinant> make_hilbert_space(size_t nmo, size_t na, size_t nb, size_t nirrep = 1,
+                                            std::vector<int> mo_symmetry = std::vector<int>(),
+                                            int symmetry = 0);
+
 } // namespace forte