diff --git a/forte/CMakeLists.txt b/forte/CMakeLists.txt index 892724f04..cc43e2751 100644 --- 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 @@ -275,7 +274,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 index 22dd5d29e..b62b31a12 100644 --- a/forte/api/forte_python_module.cc +++ b/forte/api/forte_python_module.cc @@ -214,7 +214,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 index fe587a49a..37c4d477e 100644 --- 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 index 1c61131a2..e69de29bb 100644 --- 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 -#include - -#include "sparse_ci/sparse_hamiltonian.h" - -namespace py = pybind11; -using namespace pybind11::literals; - -namespace forte { - -void export_SparseHamiltonian(py::module& m) { - py::class_(m, "SparseHamiltonian", - "A class to represent a sparse Hamiltonian") - .def(py::init>()) - .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 index 73e3a3e73..bbdf111e1 100644 --- 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 index 756d50bc4..bf0232703 100644 --- a/forte/api/sparse_operator_list_api.cc +++ b/forte/api/sparse_operator_list_api.cc @@ -164,7 +164,14 @@ void export_SparseOperatorList(py::module& m) { auto sop = op.to_operator(); return apply_operator_lin(sop, st); }, - "Multiply a SparseOperator and a SparseState"); + "Multiply a SparseOperator and a SparseState") + .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/sparse_ci/sparse_hamiltonian.cc b/forte/attic/sparse_hamiltonian.cc similarity index 100% rename from forte/sparse_ci/sparse_hamiltonian.cc rename to forte/attic/sparse_hamiltonian.cc diff --git a/forte/sparse_ci/sparse_hamiltonian.h b/forte/attic/sparse_hamiltonian.h similarity index 100% rename from forte/sparse_ci/sparse_hamiltonian.h rename to forte/attic/sparse_hamiltonian.h diff --git a/forte/modules/general_cc.py b/forte/modules/general_cc.py index 7eb24f9c1..2b41db26f 100644 --- a/forte/modules/general_cc.py +++ b/forte/modules/general_cc.py @@ -15,10 +15,10 @@ SparseOperatorList, Determinant, SparseState, - SparseHamiltonian, SparseFactExp, SparseExp, get_projection, + sparse_operator_hamiltonian, ) @@ -189,7 +189,7 @@ def solve_cc_equations( the number of iterations, and timings information """ diis = DIIS(t, diis_start) - ham = SparseHamiltonian(as_ints) + 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": @@ -253,19 +253,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") diff --git a/forte/sparse_ci/sparse_state_functions.cc b/forte/sparse_ci/sparse_state_functions.cc index a36e96c84..996957340 100644 --- 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 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>, + std::unordered_map>, 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 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::Hash> + std::unordered_map< + String, std::vector>, + 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 index 468b20545..22174ad8b 100644 --- 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& lhs, const std::tuple& acre, const std::vector& bcre, const std::vector& aann, const std::vector& 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 acre, @@ -74,12 +79,16 @@ SQOperatorString::SQOperatorString(const std::initializer_list acre, const std::initializer_list 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 index 2979ea8d3..09de397e3 100644 --- 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 index 399c358e0..05575cacf 100644 --- 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 index 6ced54a3c..535d9a939 100644 --- 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()