From 8a3d63cba733946e9d220dcfddc27942dc0256e6 Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Tue, 4 Mar 2025 11:38:50 -0600 Subject: [PATCH 1/4] fixes for API break --- src/scf/eigen_solver/eigen_generalized.cpp | 28 +++++++--------- src/scf/matrix_builder/density_matrix.cpp | 32 ++++++++++++------- src/scf/matrix_builder/determinant_driver.cpp | 15 ++++----- src/scf/matrix_builder/fock.cpp | 8 ++--- tests/cxx/integration_tests/guess/core.cpp | 6 ++-- .../integration_tests/matrix_builder/fock.cpp | 20 ++++++------ .../matrix_builder/scf_integrals_driver.cpp | 12 +++---- .../update/diagonalization.cpp | 6 ++-- .../eigen_solver/eigen_generalized.cpp | 6 ++-- .../matrix_builder/density_matrix.cpp | 12 +++---- 10 files changed, 70 insertions(+), 75 deletions(-) diff --git a/src/scf/eigen_solver/eigen_generalized.cpp b/src/scf/eigen_solver/eigen_generalized.cpp index 941ed6a..f84cf14 100644 --- a/src/scf/eigen_solver/eigen_generalized.cpp +++ b/src/scf/eigen_solver/eigen_generalized.cpp @@ -38,14 +38,13 @@ MODULE_RUN(EigenGeneralized) { auto&& [A, B] = pt::unwrap_inputs(inputs); // Convert to Eigen buffers - using matrix_alloc_t = tensorwrapper::allocator::Eigen; - using vector_alloc_t = tensorwrapper::allocator::Eigen; - const auto& eigen_A = matrix_alloc_t::rebind(A.buffer()); - const auto& eigen_B = matrix_alloc_t::rebind(B.buffer()); + tensorwrapper::allocator::Eigen allocator(get_runtime()); + const auto& eigen_A = allocator.rebind(A.buffer()); + const auto& eigen_B = allocator.rebind(B.buffer()); // Wrap the tensors in Eigen::Map objects to avoid copy - const auto* pA = eigen_A.value().data(); - const auto* pB = eigen_B.value().data(); + const auto* pA = eigen_A.data(); + const auto* pB = eigen_B.data(); const auto& shape_A = eigen_A.layout().shape().as_smooth(); auto rows = shape_A.extent(0); auto cols = shape_A.extent(1); @@ -63,23 +62,18 @@ MODULE_RUN(EigenGeneralized) { tensorwrapper::layout::Physical vector_layout(vector_shape); tensorwrapper::layout::Physical matrix_layout(matrix_shape); - using matrix_buffer = typename matrix_alloc_t::eigen_buffer_type; - using vector_buffer = typename vector_alloc_t::eigen_buffer_type; + auto pvalues_buffer = allocator.allocate(vector_layout); + auto pvectors_buffer = allocator.allocate(matrix_layout); - typename vector_buffer::data_type vector_tensor(rows); - typename matrix_buffer::data_type matrix_tensor(rows, cols); for(auto i = 0; i < rows; ++i) { - vector_tensor(i) = eigen_values(i); + pvalues_buffer->at(i) = eigen_values(i); for(auto j = 0; j < cols; ++j) { - matrix_tensor(i, j) = eigen_vectors(i, j); + pvectors_buffer->at(i, j) = eigen_vectors(i, j); } } - vector_buffer values_buffer(vector_tensor, vector_layout); - matrix_buffer vectors_buffer(matrix_tensor, matrix_layout); - - simde::type::tensor values(vector_shape, values_buffer); - simde::type::tensor vectors(matrix_shape, vectors_buffer); + simde::type::tensor values(vector_shape, std::move(pvalues_buffer)); + simde::type::tensor vectors(matrix_shape, std::move(pvectors_buffer)); auto rv = results(); return pt::wrap_results(rv, values, vectors); diff --git a/src/scf/matrix_builder/density_matrix.cpp b/src/scf/matrix_builder/density_matrix.cpp index 027af3f..faa33f6 100644 --- a/src/scf/matrix_builder/density_matrix.cpp +++ b/src/scf/matrix_builder/density_matrix.cpp @@ -15,6 +15,7 @@ */ #include "matrix_builder.hpp" +#include namespace scf::matrix_builder { namespace { @@ -62,27 +63,34 @@ MODULE_RUN(DensityMatrix) { "Please shuffle your orbitals so that the ensemble is a " "contiguous slice of the orbitals."); - // Step 2: Grab the orbitals in the ensemble - using allocator_type = tensorwrapper::allocator::Eigen; - using buffer_type = typename allocator_type::eigen_buffer_type; - using eigen_tensor_type = typename buffer_type::data_type; + using allocator_type = tensorwrapper::allocator::Eigen; + using tensor_type = Eigen::Tensor; + using const_tensor_type = Eigen::Tensor; + using map_type = Eigen::TensorMap; + using const_map_type = Eigen::TensorMap; - const auto& c_buffer = allocator_type::rebind(c.buffer()); - const auto& c_eigen = c_buffer.value(); + allocator_type alloc(get_runtime()); + + tensorwrapper::shape::Smooth p_shape(n_aos, n_aos); + tensorwrapper::layout::Physical l(p_shape); + auto pp_buffer = alloc.allocate(l); + + // Step 2: Grab the orbitals in the ensemble + const auto& c_buffer = alloc.rebind(c.buffer()); + const_map_type c_eigen(c_buffer.data(), n_aos, n_aos); + map_type p_eigen(pp_buffer->data(), n_aos, n_aos); using slice_t = Eigen::array; slice_t offsets{0, 0}; slice_t extents{n_aos, Eigen::Index(participants.size())}; - eigen_tensor_type slice = c_eigen.slice(offsets, extents); + auto slice = c_eigen.slice(offsets, extents); // Step 3: CC_dagger using index_pair_t = Eigen::IndexPair; Eigen::array modes{index_pair_t(1, 1)}; - eigen_tensor_type p_eigen = slice.contract(slice, modes); - tensorwrapper::shape::Smooth p_shape(n_aos, n_aos); - tensorwrapper::layout::Physical l(p_shape); - buffer_type p_buffer(p_eigen, l); - simde::type::tensor p(p_shape, p_buffer); + p_eigen = slice.contract(slice, modes); + + simde::type::tensor p(p_shape, std::move(pp_buffer)); auto rv = results(); return pt::wrap_results(rv, p); diff --git a/src/scf/matrix_builder/determinant_driver.cpp b/src/scf/matrix_builder/determinant_driver.cpp index 45ac0d6..e0516d5 100644 --- a/src/scf/matrix_builder/determinant_driver.cpp +++ b/src/scf/matrix_builder/determinant_driver.cpp @@ -103,6 +103,7 @@ MODULE_CTOR(DeterminantDriver) { } MODULE_RUN(DeterminantDriver) { + using float_type = double; using wf_type = simde::type::rscf_wf; const auto&& [braket] = pt::unwrap_inputs(inputs); const auto& bra = braket.bra(); @@ -126,18 +127,14 @@ MODULE_RUN(DeterminantDriver) { const auto& t = visitor.m_pt; // Step 3: Contract - using allocator_type = tensorwrapper::allocator::Eigen; - using scalar_type = tensorwrapper::buffer::Eigen::data_type; + tensorwrapper::Tensor x; + x("") = rho("i,j") * t("i,j"); - const auto& t_eigen = allocator_type::rebind(t.buffer()).value(); - const auto& p_eigen = allocator_type::rebind(rho.buffer()).value(); - - using index_pair_t = Eigen::IndexPair; - Eigen::array modes{index_pair_t(0, 0), index_pair_t(1, 1)}; - scalar_type x = p_eigen.contract(t_eigen, modes); + using allocator_type = tensorwrapper::allocator::Eigen; + auto& x_buffer = allocator_type::rebind(x.buffer()); auto rv = results(); - return pt::wrap_results(rv, x()); + return pt::wrap_results(rv, x_buffer.at()); } } // namespace scf::matrix_builder diff --git a/src/scf/matrix_builder/fock.cpp b/src/scf/matrix_builder/fock.cpp index fc28897..6099489 100644 --- a/src/scf/matrix_builder/fock.cpp +++ b/src/scf/matrix_builder/fock.cpp @@ -49,9 +49,7 @@ MODULE_RUN(Fock) { chemist::braket::BraKet term0(bra_aos, op_0, ket_aos); F = ao_dispatcher.run_as(term0); - using allocator_type = tensorwrapper::allocator::Eigen; - auto& f_buffer = allocator_type::rebind(F.buffer()); - f_buffer.value() = f_buffer.value() * c0; + F("i,j") = F("i,j") * c0; for(decltype(n_terms) i = 1; i < n_terms; ++i) { const auto ci = f.coefficient(i); @@ -59,8 +57,8 @@ MODULE_RUN(Fock) { chemist::braket::BraKet termi(bra_aos, op_i, ket_aos); auto matrix = ao_dispatcher.run_as(termi); - auto& matrix_buffer = allocator_type::rebind(matrix.buffer()); - f_buffer.value() += ci * matrix_buffer.value(); + matrix("i,j") = matrix("i,j") * ci; + F("i,j") = F("i,j") + matrix("i,j"); } } diff --git a/tests/cxx/integration_tests/guess/core.cpp b/tests/cxx/integration_tests/guess/core.cpp index 591c657..8ebb6f2 100644 --- a/tests/cxx/integration_tests/guess/core.cpp +++ b/tests/cxx/integration_tests/guess/core.cpp @@ -32,11 +32,11 @@ TEST_CASE("Core") { REQUIRE(psi.orbital_indices() == occs); REQUIRE(psi.orbitals().from_space() == aos); const auto& evals = psi.orbitals().diagonalized_matrix(); - using allocator_type = tensorwrapper::allocator::Eigen; + using allocator_type = tensorwrapper::allocator::Eigen; const auto& eval_buffer = allocator_type::rebind(evals.buffer()); const auto tol = 1E-6; using Catch::Matchers::WithinAbs; - REQUIRE_THAT(eval_buffer.value()(0), WithinAbs(-1.25330893, tol)); - REQUIRE_THAT(eval_buffer.value()(1), WithinAbs(-0.47506974, tol)); + REQUIRE_THAT(eval_buffer.at(0), WithinAbs(-1.25330893, tol)); + REQUIRE_THAT(eval_buffer.at(1), WithinAbs(-0.47506974, tol)); } \ No newline at end of file diff --git a/tests/cxx/integration_tests/matrix_builder/fock.cpp b/tests/cxx/integration_tests/matrix_builder/fock.cpp index 8283428..1e2b968 100644 --- a/tests/cxx/integration_tests/matrix_builder/fock.cpp +++ b/tests/cxx/integration_tests/matrix_builder/fock.cpp @@ -35,13 +35,13 @@ TEST_CASE("Fock Matrix Builder") { f_e.emplace_back(1.0, std::make_unique(e, h2)); const auto& F = mod.run_as(chemist::braket::BraKet(aos, f_e, aos)); - using alloc_type = tensorwrapper::allocator::Eigen; + using alloc_type = tensorwrapper::allocator::Eigen; const auto& F_eigen = alloc_type::rebind(F.buffer()); using Catch::Matchers::WithinAbs; - REQUIRE_THAT(F_eigen.value()(0, 0), WithinAbs(-1.120958, 1E-6)); - REQUIRE_THAT(F_eigen.value()(0, 1), WithinAbs(-0.959374, 1E-6)); - REQUIRE_THAT(F_eigen.value()(1, 0), WithinAbs(-0.959374, 1E-6)); - REQUIRE_THAT(F_eigen.value()(1, 1), WithinAbs(-1.120958, 1E-6)); + REQUIRE_THAT(F_eigen.at(0, 0), WithinAbs(-1.120958, 1E-6)); + REQUIRE_THAT(F_eigen.at(0, 1), WithinAbs(-0.959374, 1E-6)); + REQUIRE_THAT(F_eigen.at(1, 0), WithinAbs(-0.959374, 1E-6)); + REQUIRE_THAT(F_eigen.at(1, 1), WithinAbs(-1.120958, 1E-6)); } SECTION("With J and K") { @@ -49,13 +49,13 @@ TEST_CASE("Fock Matrix Builder") { const auto& F = mod.run_as(chemist::braket::BraKet(aos, f_e, aos)); - using alloc_type = tensorwrapper::allocator::Eigen; + using alloc_type = tensorwrapper::allocator::Eigen; const auto& F_eigen = alloc_type::rebind(F.buffer()); using Catch::Matchers::WithinAbs; - REQUIRE_THAT(F_eigen.value()(0, 0), WithinAbs(-0.319459, 1E-6)); - REQUIRE_THAT(F_eigen.value()(0, 1), WithinAbs(-0.571781, 1E-6)); - REQUIRE_THAT(F_eigen.value()(1, 0), WithinAbs(-0.571781, 1E-6)); - REQUIRE_THAT(F_eigen.value()(1, 1), WithinAbs(-0.319459, 1E-6)); + REQUIRE_THAT(F_eigen.at(0, 0), WithinAbs(-0.319459, 1E-6)); + REQUIRE_THAT(F_eigen.at(0, 1), WithinAbs(-0.571781, 1E-6)); + REQUIRE_THAT(F_eigen.at(1, 0), WithinAbs(-0.571781, 1E-6)); + REQUIRE_THAT(F_eigen.at(1, 1), WithinAbs(-0.319459, 1E-6)); } } \ No newline at end of file diff --git a/tests/cxx/integration_tests/matrix_builder/scf_integrals_driver.cpp b/tests/cxx/integration_tests/matrix_builder/scf_integrals_driver.cpp index cce4578..b2a61d5 100644 --- a/tests/cxx/integration_tests/matrix_builder/scf_integrals_driver.cpp +++ b/tests/cxx/integration_tests/matrix_builder/scf_integrals_driver.cpp @@ -23,18 +23,16 @@ namespace { void compare_matrices(const tensor& A, const tensor& A_corr) { using Catch::Matchers::WithinAbs; - using alloc_type = tensorwrapper::allocator::Eigen; + using alloc_type = tensorwrapper::allocator::Eigen; const auto& A_buffer = alloc_type::rebind(A.buffer()); const auto& A_corr_buffer = alloc_type::rebind(A_corr.buffer()); - const auto& A_eigen = A_buffer.value(); - const auto& A_corr_eigen = A_corr_buffer.value(); const auto tol = 1E-6; - REQUIRE_THAT(A_eigen(0, 0), WithinAbs(A_corr_eigen(0, 0), 1E-6)); - REQUIRE_THAT(A_eigen(0, 1), WithinAbs(A_corr_eigen(0, 1), 1E-6)); - REQUIRE_THAT(A_eigen(1, 0), WithinAbs(A_corr_eigen(1, 0), 1E-6)); - REQUIRE_THAT(A_eigen(1, 1), WithinAbs(A_corr_eigen(1, 1), 1E-6)); + REQUIRE_THAT(A_buffer.at(0, 0), WithinAbs(A_corr_buffer.at(0, 0), 1E-6)); + REQUIRE_THAT(A_buffer.at(0, 1), WithinAbs(A_corr_buffer.at(0, 1), 1E-6)); + REQUIRE_THAT(A_buffer.at(1, 0), WithinAbs(A_corr_buffer.at(1, 0), 1E-6)); + REQUIRE_THAT(A_buffer.at(1, 1), WithinAbs(A_corr_buffer.at(1, 1), 1E-6)); } } // namespace diff --git a/tests/cxx/integration_tests/update/diagonalization.cpp b/tests/cxx/integration_tests/update/diagonalization.cpp index 4b46ec0..c174044 100644 --- a/tests/cxx/integration_tests/update/diagonalization.cpp +++ b/tests/cxx/integration_tests/update/diagonalization.cpp @@ -47,12 +47,12 @@ TEST_CASE("Diagaonalization") { // Check orbital energies const auto& evals = psi.orbitals().diagonalized_matrix(); - using vector_allocator = tensorwrapper::allocator::Eigen; + using vector_allocator = tensorwrapper::allocator::Eigen; const auto& eval_buffer = vector_allocator::rebind(evals.buffer()); const auto tol = 1E-6; using Catch::Matchers::WithinAbs; - REQUIRE_THAT(eval_buffer.value()(0), WithinAbs(-1.25330893, tol)); - REQUIRE_THAT(eval_buffer.value()(1), WithinAbs(-0.47506974, tol)); + REQUIRE_THAT(eval_buffer.at(0), WithinAbs(-1.25330893, tol)); + REQUIRE_THAT(eval_buffer.at(1), WithinAbs(-0.47506974, tol)); } } \ No newline at end of file diff --git a/tests/cxx/unit_tests/eigen_solver/eigen_generalized.cpp b/tests/cxx/unit_tests/eigen_solver/eigen_generalized.cpp index 32e1ad7..577cdff 100644 --- a/tests/cxx/unit_tests/eigen_solver/eigen_generalized.cpp +++ b/tests/cxx/unit_tests/eigen_solver/eigen_generalized.cpp @@ -31,10 +31,10 @@ TEST_CASE("EigenGeneralized") { auto&& [values, vector] = mod.run_as(A, B); - using value_alloc_t = tensorwrapper::allocator::Eigen; + using value_alloc_t = tensorwrapper::allocator::Eigen; const auto& eigen_values = value_alloc_t::rebind(values.buffer()); using Catch::Matchers::WithinAbs; - REQUIRE_THAT(eigen_values.value()(0), WithinAbs(-0.236068, 1E-6)); - REQUIRE_THAT(eigen_values.value()(1), WithinAbs(4.236068, 1E-6)); + REQUIRE_THAT(eigen_values.at(0), WithinAbs(-0.236068, 1E-6)); + REQUIRE_THAT(eigen_values.at(1), WithinAbs(4.236068, 1E-6)); } \ No newline at end of file diff --git a/tests/cxx/unit_tests/matrix_builder/density_matrix.cpp b/tests/cxx/unit_tests/matrix_builder/density_matrix.cpp index bad620d..b050194 100644 --- a/tests/cxx/unit_tests/matrix_builder/density_matrix.cpp +++ b/tests/cxx/unit_tests/matrix_builder/density_matrix.cpp @@ -29,12 +29,12 @@ TEST_CASE("Density Matrix Builder") { chemist::braket::BraKet p_mn(aos, rho_hat, aos); const auto& P = mod.run_as(p_mn); - using allocator_type = tensorwrapper::allocator::Eigen; - const auto& P_eigen = allocator_type::rebind(P.buffer()).value(); + using allocator_type = tensorwrapper::allocator::Eigen; + const auto& P_eigen = allocator_type::rebind(P.buffer()); using Catch::Matchers::WithinAbs; - REQUIRE_THAT(P_eigen(0, 0), WithinAbs(0.31980835, 1E-6)); - REQUIRE_THAT(P_eigen(0, 1), WithinAbs(0.31980835, 1E-6)); - REQUIRE_THAT(P_eigen(1, 0), WithinAbs(0.31980835, 1E-6)); - REQUIRE_THAT(P_eigen(1, 1), WithinAbs(0.31980835, 1E-6)); + REQUIRE_THAT(P_eigen.at(0, 0), WithinAbs(0.31980835, 1E-6)); + REQUIRE_THAT(P_eigen.at(0, 1), WithinAbs(0.31980835, 1E-6)); + REQUIRE_THAT(P_eigen.at(1, 0), WithinAbs(0.31980835, 1E-6)); + REQUIRE_THAT(P_eigen.at(1, 1), WithinAbs(0.31980835, 1E-6)); } \ No newline at end of file From ec66623f721066fa3828068d4dab1dabaa456c97 Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Tue, 4 Mar 2025 22:28:58 -0600 Subject: [PATCH 2/4] most of SCF works with UQ [skip ci] --- src/scf/eigen_solver/eigen_generalized.cpp | 59 ++++++++++------ src/scf/matrix_builder/density_matrix.cpp | 68 +++++++++++-------- tests/cxx/integration_tests/coulombs_law.cpp | 5 +- .../integration_tests/driver/scf_driver.cpp | 7 +- .../cxx/integration_tests/driver/scf_loop.cpp | 9 +-- tests/cxx/integration_tests/guess/core.cpp | 16 +++-- .../integration_tests/integration_tests.hpp | 15 +++- .../matrix_builder/determinant_driver.cpp | 11 +-- .../matrix_builder/electronic_energy.cpp | 9 +-- .../integration_tests/matrix_builder/fock.cpp | 41 ++++++----- .../matrix_builder/scf_integrals_driver.cpp | 61 +++++++++-------- .../update/diagonalization.cpp | 17 +++-- tests/cxx/test_scf.hpp | 52 ++++++++++---- .../unit_tests/fock_operator/restricted.cpp | 16 ++--- .../matrix_builder/density_matrix.cpp | 3 +- 15 files changed, 241 insertions(+), 148 deletions(-) diff --git a/src/scf/eigen_solver/eigen_generalized.cpp b/src/scf/eigen_solver/eigen_generalized.cpp index f84cf14..fffab6b 100644 --- a/src/scf/eigen_solver/eigen_generalized.cpp +++ b/src/scf/eigen_solver/eigen_generalized.cpp @@ -22,23 +22,12 @@ namespace scf::eigen_solver { using pt = simde::GeneralizedEigenSolve; -const auto desc = R"( -Generalized Eigen Solve via Eigen ---------------------------------- - -TODO: Write me!!! -)"; - -MODULE_CTOR(EigenGeneralized) { - description(desc); - satisfies_property_type(); -} - -MODULE_RUN(EigenGeneralized) { - auto&& [A, B] = pt::unwrap_inputs(inputs); - +template +auto eigen_kernel(const tensorwrapper::Tensor& A, + const tensorwrapper::Tensor& B, + parallelzone::runtime::RuntimeView& rv) { // Convert to Eigen buffers - tensorwrapper::allocator::Eigen allocator(get_runtime()); + tensorwrapper::allocator::Eigen allocator(rv); const auto& eigen_A = allocator.rebind(A.buffer()); const auto& eigen_B = allocator.rebind(B.buffer()); @@ -48,11 +37,17 @@ MODULE_RUN(EigenGeneralized) { const auto& shape_A = eigen_A.layout().shape().as_smooth(); auto rows = shape_A.extent(0); auto cols = shape_A.extent(1); - Eigen::Map A_map(pA, rows, cols); - Eigen::Map B_map(pB, rows, cols); + + constexpr auto rmajor = Eigen::RowMajor; + constexpr auto edynam = Eigen::Dynamic; + using matrix_type = Eigen::Matrix; + using map_type = Eigen::Map; + + map_type A_map(pA, rows, cols); + map_type B_map(pB, rows, cols); // Compute - Eigen::GeneralizedSelfAdjointEigenSolver ges(A_map, B_map); + Eigen::GeneralizedSelfAdjointEigenSolver ges(A_map, B_map); auto eigen_values = ges.eigenvalues(); auto eigen_vectors = ges.eigenvectors(); @@ -74,6 +69,32 @@ MODULE_RUN(EigenGeneralized) { simde::type::tensor values(vector_shape, std::move(pvalues_buffer)); simde::type::tensor vectors(matrix_shape, std::move(pvectors_buffer)); + return std::make_pair(values, vectors); +} + +const auto desc = R"( +Generalized Eigen Solve via Eigen +--------------------------------- + +TODO: Write me!!! +)"; + +MODULE_CTOR(EigenGeneralized) { + description(desc); + satisfies_property_type(); +} + +MODULE_RUN(EigenGeneralized) { + auto&& [A, B] = pt::unwrap_inputs(inputs); + + simde::type::tensor values; + simde::type::tensor vectors; + if(tensorwrapper::allocator::Eigen::can_rebind(A.buffer())) { + std::tie(values, vectors) = eigen_kernel(A, B, get_runtime()); + } else { + using udouble = tensorwrapper::types::udouble; + std::tie(values, vectors) = eigen_kernel(A, B, get_runtime()); + } auto rv = results(); return pt::wrap_results(rv, values, vectors); diff --git a/src/scf/matrix_builder/density_matrix.cpp b/src/scf/matrix_builder/density_matrix.cpp index faa33f6..4e33fd5 100644 --- a/src/scf/matrix_builder/density_matrix.cpp +++ b/src/scf/matrix_builder/density_matrix.cpp @@ -27,6 +27,37 @@ const auto desc = R"( using pt = simde::aos_rho_e_aos; +template +auto build_density(std::size_t n_aos, std::size_t n_occ, + const tensorwrapper::Tensor& c) { + constexpr auto rmajor = Eigen::RowMajor; + constexpr auto edynam = Eigen::Dynamic; + using allocator_type = tensorwrapper::allocator::Eigen; + using tensor_type = Eigen::Matrix; + using map_type = Eigen::Map; + using const_map_type = Eigen::Map; + auto rv = c.buffer().allocator().runtime(); + allocator_type alloc(rv); + + tensorwrapper::shape::Smooth p_shape(n_aos, n_aos); + tensorwrapper::layout::Physical l(p_shape); + auto pp_buffer = alloc.allocate(l); + + // Step 2: Grab the orbitals in the ensemble + auto& c_buffer = alloc.rebind(c.buffer()); + + const_map_type c_eigen(c_buffer.data(), n_aos, n_aos); + map_type p_eigen(pp_buffer->data(), n_aos, n_aos); + auto slice = c_eigen.block(0, 0, n_aos, n_occ); + + // Step 3: CC_dagger + using index_pair_t = Eigen::IndexPair; + Eigen::array modes{index_pair_t(1, 1)}; + p_eigen = slice * slice.transpose(); + + return simde::type::tensor(p_shape, std::move(pp_buffer)); +} + MODULE_CTOR(DensityMatrix) { description(desc); satisfies_property_type(); @@ -63,34 +94,15 @@ MODULE_RUN(DensityMatrix) { "Please shuffle your orbitals so that the ensemble is a " "contiguous slice of the orbitals."); - using allocator_type = tensorwrapper::allocator::Eigen; - using tensor_type = Eigen::Tensor; - using const_tensor_type = Eigen::Tensor; - using map_type = Eigen::TensorMap; - using const_map_type = Eigen::TensorMap; - - allocator_type alloc(get_runtime()); - - tensorwrapper::shape::Smooth p_shape(n_aos, n_aos); - tensorwrapper::layout::Physical l(p_shape); - auto pp_buffer = alloc.allocate(l); - - // Step 2: Grab the orbitals in the ensemble - const auto& c_buffer = alloc.rebind(c.buffer()); - const_map_type c_eigen(c_buffer.data(), n_aos, n_aos); - map_type p_eigen(pp_buffer->data(), n_aos, n_aos); - - using slice_t = Eigen::array; - slice_t offsets{0, 0}; - slice_t extents{n_aos, Eigen::Index(participants.size())}; - auto slice = c_eigen.slice(offsets, extents); - - // Step 3: CC_dagger - using index_pair_t = Eigen::IndexPair; - Eigen::array modes{index_pair_t(1, 1)}; - p_eigen = slice.contract(slice, modes); - - simde::type::tensor p(p_shape, std::move(pp_buffer)); + // TODO: The need to dispatch like this goes away when TW supports slicing + simde::type::tensor p; + using double_alloc = tensorwrapper::allocator::Eigen; + if(double_alloc::can_rebind(c.buffer())) { + p = build_density(n_aos, participants.size(), c); + } else { + using udouble = tensorwrapper::types::udouble; + p = build_density(n_aos, participants.size(), c); + } auto rv = results(); return pt::wrap_results(rv, p); diff --git a/tests/cxx/integration_tests/coulombs_law.cpp b/tests/cxx/integration_tests/coulombs_law.cpp index bf7525c..c9cc280 100644 --- a/tests/cxx/integration_tests/coulombs_law.cpp +++ b/tests/cxx/integration_tests/coulombs_law.cpp @@ -20,8 +20,9 @@ using pt = simde::charge_charge_interaction; using Catch::Matchers::WithinAbs; TEST_CASE("CoulombsLaw") { - auto mm = test_scf::load_modules(); - auto& mod = mm.at("Coulomb's Law"); + using float_type = double; + auto mm = test_scf::load_modules(); + auto& mod = mm.at("Coulomb's Law"); auto h2_nuclei = test_scf::make_h2(); // TODO: Conversions are missing in Chemist. Use those when they're in place diff --git a/tests/cxx/integration_tests/driver/scf_driver.cpp b/tests/cxx/integration_tests/driver/scf_driver.cpp index 417f187..dd69d9e 100644 --- a/tests/cxx/integration_tests/driver/scf_driver.cpp +++ b/tests/cxx/integration_tests/driver/scf_driver.cpp @@ -21,9 +21,10 @@ using Catch::Matchers::WithinAbs; using pt = simde::AOEnergy; TEST_CASE("SCFDriver") { - auto mm = test_scf::load_modules(); - auto h2 = test_scf::make_h2(); - auto aos = test_scf::h2_aos().ao_basis_set(); + using float_type = double; + auto mm = test_scf::load_modules(); + auto h2 = test_scf::make_h2(); + auto aos = test_scf::h2_aos().ao_basis_set(); const auto e = mm.run_as("SCF Driver", aos, h2); REQUIRE_THAT(e, WithinAbs(-1.1167592336, 1E-6)); diff --git a/tests/cxx/integration_tests/driver/scf_loop.cpp b/tests/cxx/integration_tests/driver/scf_loop.cpp index 5c9cbfd..f4ac385 100644 --- a/tests/cxx/integration_tests/driver/scf_loop.cpp +++ b/tests/cxx/integration_tests/driver/scf_loop.cpp @@ -25,12 +25,13 @@ template using pt = simde::Optimize, WFType>; TEST_CASE("SCFLoop") { - using wf_type = simde::type::rscf_wf; - auto mm = test_scf::load_modules(); - auto& mod = mm.at("Loop"); + using float_type = double; + using wf_type = simde::type::rscf_wf; + auto mm = test_scf::load_modules(); + auto& mod = mm.at("Loop"); using index_set = typename wf_type::orbital_index_set_type; - wf_type psi0(index_set{0}, test_scf::h2_cmos()); + wf_type psi0(index_set{0}, test_scf::h2_cmos()); auto H = test_scf::h2_hamiltonian(); diff --git a/tests/cxx/integration_tests/guess/core.cpp b/tests/cxx/integration_tests/guess/core.cpp index 8ebb6f2..66a7766 100644 --- a/tests/cxx/integration_tests/guess/core.cpp +++ b/tests/cxx/integration_tests/guess/core.cpp @@ -20,23 +20,27 @@ using rscf_wf = simde::type::rscf_wf; using pt = simde::InitialGuess; using simde::type::tensor; -TEST_CASE("Core") { - auto mm = test_scf::load_modules(); +TEMPLATE_LIST_TEST_CASE("Core", "", test_scf::float_types) { + using float_type = TestType; + + auto mm = test_scf::load_modules(); auto aos = test_scf::h2_aos(); auto H = test_scf::h2_hamiltonian(); auto mod = mm.at("Core guess"); - auto psi = mod.run_as(H, aos); + auto psi = mod.template run_as(H, aos); typename rscf_wf::orbital_index_set_type occs{0}; REQUIRE(psi.orbital_indices() == occs); REQUIRE(psi.orbitals().from_space() == aos); const auto& evals = psi.orbitals().diagonalized_matrix(); - using allocator_type = tensorwrapper::allocator::Eigen; + using allocator_type = tensorwrapper::allocator::Eigen; const auto& eval_buffer = allocator_type::rebind(evals.buffer()); const auto tol = 1E-6; using Catch::Matchers::WithinAbs; - REQUIRE_THAT(eval_buffer.at(0), WithinAbs(-1.25330893, tol)); - REQUIRE_THAT(eval_buffer.at(1), WithinAbs(-0.47506974, tol)); + if constexpr(std::is_same_v) { + REQUIRE_THAT(eval_buffer.at(0), WithinAbs(-1.25330893, tol)); + REQUIRE_THAT(eval_buffer.at(1), WithinAbs(-0.47506974, tol)); + } } \ No newline at end of file diff --git a/tests/cxx/integration_tests/integration_tests.hpp b/tests/cxx/integration_tests/integration_tests.hpp index 36ed094..1f879b2 100644 --- a/tests/cxx/integration_tests/integration_tests.hpp +++ b/tests/cxx/integration_tests/integration_tests.hpp @@ -23,8 +23,12 @@ namespace test_scf { +/// Floating point types to test +using float_types = std::tuple; + /// Factors out setting submodules for SCF plugin from other plugins -inline auto load_modules() { +template +pluginplay::ModuleManager load_modules() { auto rv = std::make_shared(); pluginplay::ModuleManager mm(rv, nullptr); scf::load_modules(mm); @@ -44,6 +48,15 @@ inline auto load_modules() { mm.change_submod("Diagonalization Fock update", "Overlap matrix builder", "Overlap"); + if constexpr(!std::is_same_v) { + mm.change_input("Evaluate 2-Index BraKet", "With UQ?", true); + mm.change_input("Evaluate 4-Index BraKet", "With UQ?", true); + mm.change_input("Overlap", "With UQ?", true); + mm.change_input("ERI4", "With UQ?", true); + mm.change_input("Kinetic", "With UQ?", true); + mm.change_input("Nuclear", "With UQ?", true); + } + return mm; } diff --git a/tests/cxx/integration_tests/matrix_builder/determinant_driver.cpp b/tests/cxx/integration_tests/matrix_builder/determinant_driver.cpp index 3655ed8..5ab85e4 100644 --- a/tests/cxx/integration_tests/matrix_builder/determinant_driver.cpp +++ b/tests/cxx/integration_tests/matrix_builder/determinant_driver.cpp @@ -27,13 +27,14 @@ using erased_type = chemist::braket::BraKet; TEST_CASE("DeterminantDriver") { - auto mm = test_scf::load_modules(); - auto mod = mm.at("Determinant driver"); + using float_type = double; + auto mm = test_scf::load_modules(); + auto mod = mm.at("Determinant driver"); using wf_type = simde::type::rscf_wf; using index_set = typename wf_type::orbital_index_set_type; - wf_type psi(index_set{0}, test_scf::h2_cmos()); + wf_type psi(index_set{0}, test_scf::h2_cmos()); simde::type::many_electrons es{2}; SECTION("Calling Kinetic") { @@ -54,7 +55,7 @@ TEST_CASE("DeterminantDriver") { } SECTION("Calling J") { - auto rho = test_scf::h2_density(); + auto rho = test_scf::h2_density(); simde::type::J_e_type J_e(es, rho); chemist::braket::BraKet braket(psi, J_e, psi); erased_type copy_braket(braket); @@ -63,7 +64,7 @@ TEST_CASE("DeterminantDriver") { } SECTION("Calling K") { - auto rho = test_scf::h2_density(); + auto rho = test_scf::h2_density(); simde::type::K_e_type K_e(es, rho); chemist::braket::BraKet braket(psi, K_e, psi); erased_type copy_braket(braket); diff --git a/tests/cxx/integration_tests/matrix_builder/electronic_energy.cpp b/tests/cxx/integration_tests/matrix_builder/electronic_energy.cpp index 28406f2..7fd7e75 100644 --- a/tests/cxx/integration_tests/matrix_builder/electronic_energy.cpp +++ b/tests/cxx/integration_tests/matrix_builder/electronic_energy.cpp @@ -22,13 +22,14 @@ using pt = simde::eval_braket; TEST_CASE("ElectronicEnergy") { - auto mm = test_scf::load_modules(); - auto mod = mm.at("Electronic energy"); + using float_type = double; + auto mm = test_scf::load_modules(); + auto mod = mm.at("Electronic energy"); using wf_type = simde::type::rscf_wf; using index_set = typename wf_type::orbital_index_set_type; - wf_type psi(index_set{0}, test_scf::h2_cmos()); + wf_type psi(index_set{0}, test_scf::h2_cmos()); simde::type::many_electrons es{2}; simde::type::T_e_type T_e(es); @@ -36,7 +37,7 @@ TEST_CASE("ElectronicEnergy") { auto h2_nuclei = test_scf::make_h2(); simde::type::V_en_type V_en(es, h2_nuclei); - auto rho = test_scf::h2_density(); + auto rho = test_scf::h2_density(); simde::type::J_e_type J_e(es, rho); simde::type::K_e_type K_e(es, rho); simde::type::electronic_hamiltonian H_e(T_e * 2.0 + V_en * 2.0 + J_e * 2.0 - diff --git a/tests/cxx/integration_tests/matrix_builder/fock.cpp b/tests/cxx/integration_tests/matrix_builder/fock.cpp index 1e2b968..6054627 100644 --- a/tests/cxx/integration_tests/matrix_builder/fock.cpp +++ b/tests/cxx/integration_tests/matrix_builder/fock.cpp @@ -21,8 +21,10 @@ using pt = simde::aos_f_e_aos; using simde::type::t_e_type; using simde::type::v_en_type; -TEST_CASE("Fock Matrix Builder") { - auto mm = test_scf::load_modules(); +TEMPLATE_LIST_TEST_CASE("Fock Matrix Builder", "", test_scf::float_types) { + using float_type = TestType; + auto mm = test_scf::load_modules(); + auto& mod = mm.at("Fock Matrix Builder"); auto aos = test_scf::h2_aos(); @@ -33,29 +35,34 @@ TEST_CASE("Fock Matrix Builder") { simde::type::fock f_e; f_e.emplace_back(1.0, std::make_unique(e)); f_e.emplace_back(1.0, std::make_unique(e, h2)); - const auto& F = mod.run_as(chemist::braket::BraKet(aos, f_e, aos)); + chemist::braket::BraKet f_mn(aos, f_e, aos); + const auto& F = mod.template run_as(f_mn); - using alloc_type = tensorwrapper::allocator::Eigen; + using alloc_type = tensorwrapper::allocator::Eigen; const auto& F_eigen = alloc_type::rebind(F.buffer()); using Catch::Matchers::WithinAbs; - REQUIRE_THAT(F_eigen.at(0, 0), WithinAbs(-1.120958, 1E-6)); - REQUIRE_THAT(F_eigen.at(0, 1), WithinAbs(-0.959374, 1E-6)); - REQUIRE_THAT(F_eigen.at(1, 0), WithinAbs(-0.959374, 1E-6)); - REQUIRE_THAT(F_eigen.at(1, 1), WithinAbs(-1.120958, 1E-6)); + if constexpr(std::is_same_v) { + REQUIRE_THAT(F_eigen.at(0, 0), WithinAbs(-1.120958, 1E-6)); + REQUIRE_THAT(F_eigen.at(0, 1), WithinAbs(-0.959374, 1E-6)); + REQUIRE_THAT(F_eigen.at(1, 0), WithinAbs(-0.959374, 1E-6)); + REQUIRE_THAT(F_eigen.at(1, 1), WithinAbs(-1.120958, 1E-6)); + } } SECTION("With J and K") { - auto f_e = test_scf::h2_fock(); - - const auto& F = mod.run_as(chemist::braket::BraKet(aos, f_e, aos)); + auto f_e = test_scf::h2_fock(); + chemist::braket::BraKet f_mn(aos, f_e, aos); + const auto& F = mod.template run_as(f_mn); - using alloc_type = tensorwrapper::allocator::Eigen; + using alloc_type = tensorwrapper::allocator::Eigen; const auto& F_eigen = alloc_type::rebind(F.buffer()); using Catch::Matchers::WithinAbs; - REQUIRE_THAT(F_eigen.at(0, 0), WithinAbs(-0.319459, 1E-6)); - REQUIRE_THAT(F_eigen.at(0, 1), WithinAbs(-0.571781, 1E-6)); - REQUIRE_THAT(F_eigen.at(1, 0), WithinAbs(-0.571781, 1E-6)); - REQUIRE_THAT(F_eigen.at(1, 1), WithinAbs(-0.319459, 1E-6)); + if constexpr(std::is_same_v) { + REQUIRE_THAT(F_eigen.at(0, 0), WithinAbs(-0.319459, 1E-6)); + REQUIRE_THAT(F_eigen.at(0, 1), WithinAbs(-0.571781, 1E-6)); + REQUIRE_THAT(F_eigen.at(1, 0), WithinAbs(-0.571781, 1E-6)); + REQUIRE_THAT(F_eigen.at(1, 1), WithinAbs(-0.319459, 1E-6)); + } } -} \ No newline at end of file +} diff --git a/tests/cxx/integration_tests/matrix_builder/scf_integrals_driver.cpp b/tests/cxx/integration_tests/matrix_builder/scf_integrals_driver.cpp index b2a61d5..90dbcd7 100644 --- a/tests/cxx/integration_tests/matrix_builder/scf_integrals_driver.cpp +++ b/tests/cxx/integration_tests/matrix_builder/scf_integrals_driver.cpp @@ -21,18 +21,24 @@ using simde::type::tensor; namespace { +template void compare_matrices(const tensor& A, const tensor& A_corr) { using Catch::Matchers::WithinAbs; - using alloc_type = tensorwrapper::allocator::Eigen; + using alloc_type = tensorwrapper::allocator::Eigen; const auto& A_buffer = alloc_type::rebind(A.buffer()); const auto& A_corr_buffer = alloc_type::rebind(A_corr.buffer()); const auto tol = 1E-6; - - REQUIRE_THAT(A_buffer.at(0, 0), WithinAbs(A_corr_buffer.at(0, 0), 1E-6)); - REQUIRE_THAT(A_buffer.at(0, 1), WithinAbs(A_corr_buffer.at(0, 1), 1E-6)); - REQUIRE_THAT(A_buffer.at(1, 0), WithinAbs(A_corr_buffer.at(1, 0), 1E-6)); - REQUIRE_THAT(A_buffer.at(1, 1), WithinAbs(A_corr_buffer.at(1, 1), 1E-6)); + if constexpr(std::is_same_v) { + REQUIRE_THAT(A_buffer.at(0, 0), + WithinAbs(A_corr_buffer.at(0, 0), 1E-6)); + REQUIRE_THAT(A_buffer.at(0, 1), + WithinAbs(A_corr_buffer.at(0, 1), 1E-6)); + REQUIRE_THAT(A_buffer.at(1, 0), + WithinAbs(A_corr_buffer.at(1, 0), 1E-6)); + REQUIRE_THAT(A_buffer.at(1, 1), + WithinAbs(A_corr_buffer.at(1, 1), 1E-6)); + } } } // namespace @@ -41,21 +47,22 @@ using erased_type = chemist::braket::BraKet; -TEST_CASE("SCFIntegralsDriver") { - auto mm = test_scf::load_modules(); - auto aos = test_scf::h2_aos(); - auto mod = mm.at("SCF integral driver"); +TEMPLATE_LIST_TEST_CASE("SCFIntegralsDriver", "", test_scf::float_types) { + using float_type = TestType; + auto mm = test_scf::load_modules(); + auto aos = test_scf::h2_aos(); + auto mod = mm.at("SCF integral driver"); simde::type::electron e; - auto rho = test_scf::h2_density(); + auto rho = test_scf::h2_density(); SECTION("Calling Kinetic") { auto& tmod = mm.at("Kinetic"); simde::type::t_e_type t_e(e); chemist::braket::BraKet braket(aos, t_e, aos); erased_type copy_braket(braket); - const auto& T = mod.run_as(copy_braket); - const auto& T_corr = tmod.run_as(braket); - compare_matrices(T, T_corr); + const auto& T = mod.template run_as(copy_braket); + const auto& T_corr = tmod.template run_as(braket); + compare_matrices(T, T_corr); } SECTION("Calling Electron-Nuclear Attraction") { @@ -64,9 +71,9 @@ TEST_CASE("SCFIntegralsDriver") { simde::type::v_en_type v_en(e, h2_nuclei); chemist::braket::BraKet braket(aos, v_en, aos); erased_type copy_braket(braket); - const auto& V = mod.run_as(copy_braket); - const auto& V_corr = tmod.run_as(braket); - compare_matrices(V, V_corr); + const auto& V = mod.template run_as(copy_braket); + const auto& V_corr = tmod.template run_as(braket); + compare_matrices(V, V_corr); } SECTION("Calling J Matrix") { @@ -74,9 +81,9 @@ TEST_CASE("SCFIntegralsDriver") { simde::type::j_e_type j_e(e, rho); chemist::braket::BraKet braket(aos, j_e, aos); erased_type copy_braket(braket); - const auto& J = mod.run_as(copy_braket); - const auto& J_corr = jmod.run_as(braket); - compare_matrices(J, J_corr); + const auto& J = mod.template run_as(copy_braket); + const auto& J_corr = jmod.template run_as(braket); + compare_matrices(J, J_corr); } SECTION("Calling K Matrix") { @@ -84,22 +91,22 @@ TEST_CASE("SCFIntegralsDriver") { simde::type::k_e_type k_e(e, rho); chemist::braket::BraKet braket(aos, k_e, aos); erased_type copy_braket(braket); - const auto& K = mod.run_as(copy_braket); - const auto& K_corr = kmod.run_as(braket); - compare_matrices(K, K_corr); + const auto& K = mod.template run_as(copy_braket); + const auto& K_corr = kmod.template run_as(braket); + compare_matrices(K, K_corr); } SECTION("Calling density matrix") { auto& pmod = mm.at("Density matrix builder"); - auto cmos = test_scf::h2_cmos(); + auto cmos = test_scf::h2_cmos(); std::vector occs{1, 0}; simde::type::rho_e rho_hat(cmos, occs); chemist::braket::BraKet braket(aos, rho_hat, aos); erased_type copy_braket(braket); - const auto& P = mod.run_as(copy_braket); + const auto& P = mod.template run_as(copy_braket); using op_pt = simde::aos_rho_e_aos; - const auto& P_corr = pmod.run_as(braket); - compare_matrices(P, P_corr); + const auto& P_corr = pmod.template run_as(braket); + compare_matrices(P, P_corr); } // SECTION("Calling Fock Matrix") { diff --git a/tests/cxx/integration_tests/update/diagonalization.cpp b/tests/cxx/integration_tests/update/diagonalization.cpp index c174044..4ba00d3 100644 --- a/tests/cxx/integration_tests/update/diagonalization.cpp +++ b/tests/cxx/integration_tests/update/diagonalization.cpp @@ -23,9 +23,10 @@ using orbital_index_set = typename wf_t::orbital_index_set_type; using simde::type::t_e_type; using simde::type::v_en_type; -TEST_CASE("Diagaonalization") { - auto mm = test_scf::load_modules(); - auto& mod = mm.at("Diagonalization Fock update"); +TEMPLATE_LIST_TEST_CASE("Diagaonalization", "", test_scf::float_types) { + using float_type = TestType; + auto mm = test_scf::load_modules(); + auto& mod = mm.at("Diagonalization Fock update"); auto aos = test_scf::h2_aos(); auto h2 = test_scf::make_h2(); @@ -41,18 +42,20 @@ TEST_CASE("Diagaonalization") { simde::type::cmos cmos(empty, aos, empty); simde::type::rscf_wf core_guess(occs, cmos); - const auto& psi = mod.run_as(f_e, core_guess); + const auto& psi = mod.template run_as(f_e, core_guess); REQUIRE(psi.orbital_indices() == occs); REQUIRE(psi.orbitals().from_space() == aos); // Check orbital energies const auto& evals = psi.orbitals().diagonalized_matrix(); - using vector_allocator = tensorwrapper::allocator::Eigen; + using vector_allocator = tensorwrapper::allocator::Eigen; const auto& eval_buffer = vector_allocator::rebind(evals.buffer()); const auto tol = 1E-6; using Catch::Matchers::WithinAbs; - REQUIRE_THAT(eval_buffer.at(0), WithinAbs(-1.25330893, tol)); - REQUIRE_THAT(eval_buffer.at(1), WithinAbs(-0.47506974, tol)); + if constexpr(std::is_same_v) { + REQUIRE_THAT(eval_buffer.at(0), WithinAbs(-1.25330893, tol)); + REQUIRE_THAT(eval_buffer.at(1), WithinAbs(-0.47506974, tol)); + } } } \ No newline at end of file diff --git a/tests/cxx/test_scf.hpp b/tests/cxx/test_scf.hpp index 62f8c31..fb7891b 100644 --- a/tests/cxx/test_scf.hpp +++ b/tests/cxx/test_scf.hpp @@ -91,29 +91,53 @@ inline auto h2_aos() { return simde::type::aos(h_basis(h2)); } +template inline auto h2_mos() { - using mos_type = simde::type::mos; - using tensor_type = typename mos_type::transform_type; - tensor_type c({{-0.565516, -1.07019}, {-0.565516, 1.07019}}); - return mos_type(h2_aos(), std::move(c)); + using mos_type = simde::type::mos; + using tensor_type = typename mos_type::transform_type; + using allocator_type = tensorwrapper::allocator::Eigen; + allocator_type alloc(parallelzone::runtime::RuntimeView{}); + tensorwrapper::shape::Smooth shape{2, 2}; + tensorwrapper::layout::Physical l(shape); + auto c_buffer = alloc.allocate(l); + c_buffer->at(0, 0) = -0.565516; + c_buffer->at(0, 1) = -1.07019; + c_buffer->at(1, 0) = -0.565516; + c_buffer->at(1, 1) = 1.07019; + tensor_type t(shape, std::move(c_buffer)); + return mos_type(h2_aos(), std::move(t)); } +template inline auto h2_cmos() { - using cmos_type = simde::type::cmos; - using tensor_type = typename cmos_type::transform_type; - tensor_type e({-1.25330893, -0.47506974}); - return cmos_type(e, h2_aos(), h2_mos().transform()); + using cmos_type = simde::type::cmos; + using tensor_type = typename cmos_type::transform_type; + using allocator_type = tensorwrapper::allocator::Eigen; + allocator_type alloc(parallelzone::runtime::RuntimeView{}); + tensorwrapper::shape::Smooth shape{2}; + tensorwrapper::layout::Physical l(shape); + auto e_buffer = alloc.allocate(l); + e_buffer->at(0) = -1.25330893; + e_buffer->at(1) = -0.47506974; + tensor_type e(shape, std::move(e_buffer)); + return cmos_type(std::move(e), h2_aos(), h2_mos().transform()); } +template inline auto h2_density() { - using density_type = simde::type::decomposable_e_density; - typename density_type::value_type rho( - {{0.31980835, 0.31980835}, {0.31980835, 0.31980835}}); - return density_type(rho, h2_mos()); + using density_type = simde::type::decomposable_e_density; + using tensor_type = typename density_type::value_type; + using allocator_type = tensorwrapper::allocator::Eigen; + allocator_type alloc(parallelzone::runtime::RuntimeView{}); + tensorwrapper::shape::Smooth shape{2, 2}; + tensorwrapper::layout::Physical l(shape); + auto pbuffer = alloc.construct(l, 0.31980835); + tensor_type t(shape, std::move(pbuffer)); + return density_type(std::move(t), h2_mos()); } /// The Fock matrix consistent with h2_hamiltonian and h2_density -template +template inline auto h2_fock() { ElectronType es; if constexpr(std::is_same_v) { @@ -121,7 +145,7 @@ inline auto h2_fock() { } auto h2 = make_h2(); - auto rho = h2_density(); + auto rho = h2_density(); simde::type::fock F; using namespace chemist::qm_operator; using t_type = Kinetic; diff --git a/tests/cxx/unit_tests/fock_operator/restricted.cpp b/tests/cxx/unit_tests/fock_operator/restricted.cpp index 26f7a17..7d0a522 100644 --- a/tests/cxx/unit_tests/fock_operator/restricted.cpp +++ b/tests/cxx/unit_tests/fock_operator/restricted.cpp @@ -21,19 +21,19 @@ TEST_CASE("Restricted") { pluginplay::ModuleManager mm; scf::load_modules(mm); - + using float_type = double; using density_type = simde::type::decomposable_e_density; using pt = simde::FockOperator; auto H = test_scf::h2_hamiltonian(); auto h2 = test_scf::make_h2(); - auto rho = test_scf::h2_density(); + auto rho = test_scf::h2_density(); SECTION("Many-electron version") { auto& mod = mm.at("Restricted Fock Op"); SECTION("No density") { - density_type rho; - auto F = mod.run_as(H, rho); + density_type rho_empty; + auto F = mod.run_as(H, rho_empty); simde::type::many_electrons es(2); simde::type::fock F_corr; @@ -45,8 +45,6 @@ TEST_CASE("Restricted") { } SECTION("Density") { - auto rho = test_scf::h2_density(); - auto F = mod.run_as(H, rho); auto F_corr = test_scf::h2_fock(); REQUIRE(F == F_corr); @@ -56,8 +54,8 @@ TEST_CASE("Restricted") { auto& mod = mm.at("Restricted One-Electron Fock Op"); SECTION(" No density") { - density_type rho; - auto f = mod.run_as(H, rho); + density_type rho_empty; + auto f = mod.run_as(H, rho_empty); simde::type::fock f_corr; simde::type::electron e; using simde::type::t_e_type; @@ -68,8 +66,6 @@ TEST_CASE("Restricted") { } SECTION("Density") { - auto rho = test_scf::h2_density(); - auto f = mod.run_as(H, rho); auto f_corr = test_scf::h2_fock(); REQUIRE(f == f_corr); diff --git a/tests/cxx/unit_tests/matrix_builder/density_matrix.cpp b/tests/cxx/unit_tests/matrix_builder/density_matrix.cpp index b050194..cbe232b 100644 --- a/tests/cxx/unit_tests/matrix_builder/density_matrix.cpp +++ b/tests/cxx/unit_tests/matrix_builder/density_matrix.cpp @@ -19,11 +19,12 @@ using pt = simde::aos_rho_e_aos; TEST_CASE("Density Matrix Builder") { + using float_type = double; pluginplay::ModuleManager mm; scf::load_modules(mm); auto& mod = mm.at("Density matrix builder"); auto aos = test_scf::h2_aos(); - auto cmos = test_scf::h2_cmos(); + auto cmos = test_scf::h2_cmos(); std::vector occs{1, 0}; simde::type::rho_e rho_hat(cmos, occs); From b0d02511a47b28080f654ce610e5bbbc4c7f32a9 Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Wed, 5 Mar 2025 22:33:21 -0600 Subject: [PATCH 3/4] end-to-end UQ!!! --- src/scf/driver/scf_loop.cpp | 23 +++- src/scf/eigen_solver/eigen_generalized.cpp | 115 +++++++++--------- src/scf/matrix_builder/density_matrix.cpp | 78 ++++++------ src/scf/matrix_builder/determinant_driver.cpp | 6 +- src/scf/matrix_builder/electronic_energy.cpp | 8 +- tests/cxx/integration_tests/coulombs_law.cpp | 12 +- .../integration_tests/driver/scf_driver.cpp | 17 ++- .../cxx/integration_tests/driver/scf_loop.cpp | 16 ++- tests/cxx/integration_tests/guess/core.cpp | 16 ++- .../matrix_builder/determinant_driver.cpp | 37 ++++-- .../matrix_builder/electronic_energy.cpp | 16 ++- .../integration_tests/matrix_builder/fock.cpp | 39 +++--- .../matrix_builder/scf_integrals_driver.cpp | 36 ++---- .../update/diagonalization.cpp | 18 ++- tests/cxx/test_scf.hpp | 2 + tests/cxx/unit_tests/coulombs_law.cpp | 24 +++- .../eigen_solver/eigen_generalized.cpp | 33 +++-- .../unit_tests/fock_operator/restricted.cpp | 2 +- .../matrix_builder/density_matrix.cpp | 20 +-- 19 files changed, 293 insertions(+), 225 deletions(-) diff --git a/src/scf/driver/scf_loop.cpp b/src/scf/driver/scf_loop.cpp index c005c96..f4e97d8 100644 --- a/src/scf/driver/scf_loop.cpp +++ b/src/scf/driver/scf_loop.cpp @@ -107,8 +107,8 @@ MODULE_RUN(SCFLoop) { } auto e_nuclear = V_nn_mod.run_as(qs_lhs, qs_rhs); - wf_type psi_old = psi0; - double e_old = 0; + wf_type psi_old = psi0; + simde::type::tensor e_old; const unsigned int max_iter = 3; unsigned int iter = 0; @@ -154,8 +154,23 @@ MODULE_RUN(SCFLoop) { psi_old = new_psi; ++iter; } - auto rv = results(); - return pt::wrap_results(rv, e_old + e_nuclear, psi_old); + simde::type::tensor e_total; + + // e_nuclear is a double. This hack converts it to udouble (if needed) + tensorwrapper::allocator::Eigen dalloc(get_runtime()); + using tensorwrapper::types::udouble; + tensorwrapper::allocator::Eigen ualloc(get_runtime()); + + if(ualloc.can_rebind(e_old.buffer())) { + simde::type::tensor temp(e_old); + auto val = dalloc.rebind(e_nuclear.buffer()).at(); + ualloc.rebind(temp.buffer()).at() = val; + e_nuclear = temp; + } + + e_total("") = e_old("") + e_nuclear(""); + auto rv = results(); + return pt::wrap_results(rv, e_total, psi_old); } } // namespace scf::driver \ No newline at end of file diff --git a/src/scf/eigen_solver/eigen_generalized.cpp b/src/scf/eigen_solver/eigen_generalized.cpp index fffab6b..31014dc 100644 --- a/src/scf/eigen_solver/eigen_generalized.cpp +++ b/src/scf/eigen_solver/eigen_generalized.cpp @@ -20,57 +20,61 @@ namespace scf::eigen_solver { -using pt = simde::GeneralizedEigenSolve; - -template -auto eigen_kernel(const tensorwrapper::Tensor& A, - const tensorwrapper::Tensor& B, - parallelzone::runtime::RuntimeView& rv) { - // Convert to Eigen buffers - tensorwrapper::allocator::Eigen allocator(rv); - const auto& eigen_A = allocator.rebind(A.buffer()); - const auto& eigen_B = allocator.rebind(B.buffer()); - - // Wrap the tensors in Eigen::Map objects to avoid copy - const auto* pA = eigen_A.data(); - const auto* pB = eigen_B.data(); - const auto& shape_A = eigen_A.layout().shape().as_smooth(); - auto rows = shape_A.extent(0); - auto cols = shape_A.extent(1); - - constexpr auto rmajor = Eigen::RowMajor; - constexpr auto edynam = Eigen::Dynamic; - using matrix_type = Eigen::Matrix; - using map_type = Eigen::Map; - - map_type A_map(pA, rows, cols); - map_type B_map(pB, rows, cols); - - // Compute - Eigen::GeneralizedSelfAdjointEigenSolver ges(A_map, B_map); - auto eigen_values = ges.eigenvalues(); - auto eigen_vectors = ges.eigenvectors(); - - // Wrap in TensorWrapper Tensor - tensorwrapper::shape::Smooth vector_shape{rows}; - tensorwrapper::shape::Smooth matrix_shape{rows, cols}; - tensorwrapper::layout::Physical vector_layout(vector_shape); - tensorwrapper::layout::Physical matrix_layout(matrix_shape); - - auto pvalues_buffer = allocator.allocate(vector_layout); - auto pvectors_buffer = allocator.allocate(matrix_layout); - - for(auto i = 0; i < rows; ++i) { - pvalues_buffer->at(i) = eigen_values(i); - for(auto j = 0; j < cols; ++j) { - pvectors_buffer->at(i, j) = eigen_vectors(i, j); +namespace { +struct Kernel { + template + auto run(const tensorwrapper::buffer::BufferBase& A, + const tensorwrapper::buffer::BufferBase& B, + parallelzone::runtime::RuntimeView& rv) { + // Convert to Eigen buffers + tensorwrapper::allocator::Eigen allocator(rv); + const auto& eigen_A = allocator.rebind(A); + const auto& eigen_B = allocator.rebind(B); + + // Wrap the tensors in Eigen::Map objects to avoid copy + const auto* pA = eigen_A.data(); + const auto* pB = eigen_B.data(); + const auto& shape_A = eigen_A.layout().shape().as_smooth(); + auto rows = shape_A.extent(0); + auto cols = shape_A.extent(1); + + constexpr auto rmajor = Eigen::RowMajor; + constexpr auto edynam = Eigen::Dynamic; + using matrix_type = Eigen::Matrix; + using map_type = Eigen::Map; + + map_type A_map(pA, rows, cols); + map_type B_map(pB, rows, cols); + + // Compute + Eigen::GeneralizedSelfAdjointEigenSolver ges(A_map, B_map); + auto eigen_values = ges.eigenvalues(); + auto eigen_vectors = ges.eigenvectors(); + + // Wrap in TensorWrapper Tensor + tensorwrapper::shape::Smooth vector_shape{rows}; + tensorwrapper::shape::Smooth matrix_shape{rows, cols}; + tensorwrapper::layout::Physical vector_layout(vector_shape); + tensorwrapper::layout::Physical matrix_layout(matrix_shape); + + auto pvalues_buffer = allocator.allocate(vector_layout); + auto pvectors_buffer = allocator.allocate(matrix_layout); + + for(auto i = 0; i < rows; ++i) { + pvalues_buffer->at(i) = eigen_values(i); + for(auto j = 0; j < cols; ++j) { + pvectors_buffer->at(i, j) = eigen_vectors(i, j); + } } + + simde::type::tensor values(vector_shape, std::move(pvalues_buffer)); + simde::type::tensor vectors(matrix_shape, std::move(pvectors_buffer)); + return std::make_pair(values, vectors); } +}; +} // namespace - simde::type::tensor values(vector_shape, std::move(pvalues_buffer)); - simde::type::tensor vectors(matrix_shape, std::move(pvectors_buffer)); - return std::make_pair(values, vectors); -} +using pt = simde::GeneralizedEigenSolve; const auto desc = R"( Generalized Eigen Solve via Eigen @@ -87,14 +91,13 @@ MODULE_CTOR(EigenGeneralized) { MODULE_RUN(EigenGeneralized) { auto&& [A, B] = pt::unwrap_inputs(inputs); - simde::type::tensor values; - simde::type::tensor vectors; - if(tensorwrapper::allocator::Eigen::can_rebind(A.buffer())) { - std::tie(values, vectors) = eigen_kernel(A, B, get_runtime()); - } else { - using udouble = tensorwrapper::types::udouble; - std::tie(values, vectors) = eigen_kernel(A, B, get_runtime()); - } + using tensorwrapper::utilities::floating_point_dispatch; + + auto r = get_runtime(); + Kernel k; + const auto& A_buffer = A.buffer(); + const auto& B_buffer = B.buffer(); + auto [values, vectors] = floating_point_dispatch(k, A_buffer, B_buffer, r); auto rv = results(); return pt::wrap_results(rv, values, vectors); diff --git a/src/scf/matrix_builder/density_matrix.cpp b/src/scf/matrix_builder/density_matrix.cpp index 4e33fd5..ecf9fab 100644 --- a/src/scf/matrix_builder/density_matrix.cpp +++ b/src/scf/matrix_builder/density_matrix.cpp @@ -23,40 +23,42 @@ namespace { const auto desc = R"( )"; -} +struct Kernel { + template + auto run(const tensorwrapper::buffer::BufferBase& c, std::size_t n_aos, + std::size_t n_occ) { + constexpr auto rmajor = Eigen::RowMajor; + constexpr auto edynam = Eigen::Dynamic; + using allocator_type = tensorwrapper::allocator::Eigen; + using tensor_type = Eigen::Matrix; + using map_type = Eigen::Map; + using const_map_type = Eigen::Map; + auto rv = c.allocator().runtime(); + allocator_type alloc(rv); + + tensorwrapper::shape::Smooth p_shape(n_aos, n_aos); + tensorwrapper::layout::Physical l(p_shape); + auto pp_buffer = alloc.allocate(l); + + // Step 2: Grab the orbitals in the ensemble + auto& c_buffer = alloc.rebind(c); + + const_map_type c_eigen(c_buffer.data(), n_aos, n_aos); + map_type p_eigen(pp_buffer->data(), n_aos, n_aos); + auto slice = c_eigen.block(0, 0, n_aos, n_occ); + + // Step 3: CC_dagger + using index_pair_t = Eigen::IndexPair; + Eigen::array modes{index_pair_t(1, 1)}; + p_eigen = slice * slice.transpose(); + + return simde::type::tensor(p_shape, std::move(pp_buffer)); + } +}; -using pt = simde::aos_rho_e_aos; +} // namespace -template -auto build_density(std::size_t n_aos, std::size_t n_occ, - const tensorwrapper::Tensor& c) { - constexpr auto rmajor = Eigen::RowMajor; - constexpr auto edynam = Eigen::Dynamic; - using allocator_type = tensorwrapper::allocator::Eigen; - using tensor_type = Eigen::Matrix; - using map_type = Eigen::Map; - using const_map_type = Eigen::Map; - auto rv = c.buffer().allocator().runtime(); - allocator_type alloc(rv); - - tensorwrapper::shape::Smooth p_shape(n_aos, n_aos); - tensorwrapper::layout::Physical l(p_shape); - auto pp_buffer = alloc.allocate(l); - - // Step 2: Grab the orbitals in the ensemble - auto& c_buffer = alloc.rebind(c.buffer()); - - const_map_type c_eigen(c_buffer.data(), n_aos, n_aos); - map_type p_eigen(pp_buffer->data(), n_aos, n_aos); - auto slice = c_eigen.block(0, 0, n_aos, n_occ); - - // Step 3: CC_dagger - using index_pair_t = Eigen::IndexPair; - Eigen::array modes{index_pair_t(1, 1)}; - p_eigen = slice * slice.transpose(); - - return simde::type::tensor(p_shape, std::move(pp_buffer)); -} +using pt = simde::aos_rho_e_aos; MODULE_CTOR(DensityMatrix) { description(desc); @@ -95,15 +97,9 @@ MODULE_RUN(DensityMatrix) { "contiguous slice of the orbitals."); // TODO: The need to dispatch like this goes away when TW supports slicing - simde::type::tensor p; - using double_alloc = tensorwrapper::allocator::Eigen; - if(double_alloc::can_rebind(c.buffer())) { - p = build_density(n_aos, participants.size(), c); - } else { - using udouble = tensorwrapper::types::udouble; - p = build_density(n_aos, participants.size(), c); - } - + using tensorwrapper::utilities::floating_point_dispatch; + Kernel k; + auto p = floating_point_dispatch(k, c.buffer(), n_aos, participants.size()); auto rv = results(); return pt::wrap_results(rv, p); } diff --git a/src/scf/matrix_builder/determinant_driver.cpp b/src/scf/matrix_builder/determinant_driver.cpp index e0516d5..fd7c276 100644 --- a/src/scf/matrix_builder/determinant_driver.cpp +++ b/src/scf/matrix_builder/determinant_driver.cpp @@ -103,7 +103,6 @@ MODULE_CTOR(DeterminantDriver) { } MODULE_RUN(DeterminantDriver) { - using float_type = double; using wf_type = simde::type::rscf_wf; const auto&& [braket] = pt::unwrap_inputs(inputs); const auto& bra = braket.bra(); @@ -130,11 +129,8 @@ MODULE_RUN(DeterminantDriver) { tensorwrapper::Tensor x; x("") = rho("i,j") * t("i,j"); - using allocator_type = tensorwrapper::allocator::Eigen; - auto& x_buffer = allocator_type::rebind(x.buffer()); - auto rv = results(); - return pt::wrap_results(rv, x_buffer.at()); + return pt::wrap_results(rv, x); } } // namespace scf::matrix_builder diff --git a/src/scf/matrix_builder/electronic_energy.cpp b/src/scf/matrix_builder/electronic_energy.cpp index 66b0fc6..b456549 100644 --- a/src/scf/matrix_builder/electronic_energy.cpp +++ b/src/scf/matrix_builder/electronic_energy.cpp @@ -48,7 +48,7 @@ MODULE_RUN(ElectronicEnergy) { const auto& ket = H_ij.ket(); auto& mod = submods.at("determinant driver"); - double e = 0.0; + simde::type::tensor e; auto n_ops = H.size(); for(decltype(n_ops) i = 0; i < n_ops; ++i) { @@ -57,7 +57,11 @@ MODULE_RUN(ElectronicEnergy) { chemist::braket::BraKet O_ij(bra, O_i, ket); const auto o = mod.run_as>(O_ij); - e += ci * o; + simde::type::tensor temp; + temp("") = o("") * ci; + if(i == 0) + e = temp; + else { e("") = e("") + temp(""); } } auto rv = results(); diff --git a/tests/cxx/integration_tests/coulombs_law.cpp b/tests/cxx/integration_tests/coulombs_law.cpp index c9cc280..d8c73fd 100644 --- a/tests/cxx/integration_tests/coulombs_law.cpp +++ b/tests/cxx/integration_tests/coulombs_law.cpp @@ -19,11 +19,16 @@ using pt = simde::charge_charge_interaction; using Catch::Matchers::WithinAbs; -TEST_CASE("CoulombsLaw") { +TEMPLATE_LIST_TEST_CASE("CoulombsLaw", "", test_scf::float_types) { using float_type = double; auto mm = test_scf::load_modules(); auto& mod = mm.at("Coulomb's Law"); + tensorwrapper::allocator::Eigen alloc(mm.get_runtime()); + tensorwrapper::shape::Smooth shape_corr{}; + auto pcorr = alloc.allocate(tensorwrapper::layout::Physical(shape_corr)); + using tensorwrapper::operations::approximately_equal; + auto h2_nuclei = test_scf::make_h2(); // TODO: Conversions are missing in Chemist. Use those when they're in place simde::type::charges qs; @@ -31,5 +36,8 @@ TEST_CASE("CoulombsLaw") { qs.push_back(nucleus.as_point_charge()); auto e_nuclear = mod.run_as(qs, qs); - REQUIRE_THAT(e_nuclear, WithinAbs(0.71510297482837526, 1E-6)); + + pcorr->at() = 0.71510297482837526; + tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); + REQUIRE(approximately_equal(corr, e_nuclear, 1E-6)); } \ No newline at end of file diff --git a/tests/cxx/integration_tests/driver/scf_driver.cpp b/tests/cxx/integration_tests/driver/scf_driver.cpp index dd69d9e..ae80976 100644 --- a/tests/cxx/integration_tests/driver/scf_driver.cpp +++ b/tests/cxx/integration_tests/driver/scf_driver.cpp @@ -16,16 +16,21 @@ #include "../integration_tests.hpp" -using Catch::Matchers::WithinAbs; - using pt = simde::AOEnergy; -TEST_CASE("SCFDriver") { - using float_type = double; +TEMPLATE_LIST_TEST_CASE("SCFDriver", "", test_scf::float_types) { + using float_type = TestType; auto mm = test_scf::load_modules(); auto h2 = test_scf::make_h2(); auto aos = test_scf::h2_aos().ao_basis_set(); - const auto e = mm.run_as("SCF Driver", aos, h2); - REQUIRE_THAT(e, WithinAbs(-1.1167592336, 1E-6)); + tensorwrapper::allocator::Eigen alloc(mm.get_runtime()); + tensorwrapper::shape::Smooth shape_corr{}; + auto pcorr = alloc.allocate(tensorwrapper::layout::Physical(shape_corr)); + using tensorwrapper::operations::approximately_equal; + pcorr->at() = -1.1167592336; + tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); + + const auto e = mm.template run_as("SCF Driver", aos, h2); + REQUIRE(approximately_equal(corr, e, 1E-6)); } \ No newline at end of file diff --git a/tests/cxx/integration_tests/driver/scf_loop.cpp b/tests/cxx/integration_tests/driver/scf_loop.cpp index f4ac385..51870ce 100644 --- a/tests/cxx/integration_tests/driver/scf_loop.cpp +++ b/tests/cxx/integration_tests/driver/scf_loop.cpp @@ -24,18 +24,26 @@ using egy_pt = simde::eval_braket; template using pt = simde::Optimize, WFType>; -TEST_CASE("SCFLoop") { - using float_type = double; +TEMPLATE_LIST_TEST_CASE("SCFLoop", "", test_scf::float_types) { + using float_type = TestType; using wf_type = simde::type::rscf_wf; auto mm = test_scf::load_modules(); auto& mod = mm.at("Loop"); + tensorwrapper::allocator::Eigen alloc(mm.get_runtime()); + tensorwrapper::shape::Smooth shape_corr{}; + auto pcorr = alloc.allocate(tensorwrapper::layout::Physical(shape_corr)); + using tensorwrapper::operations::approximately_equal; + using index_set = typename wf_type::orbital_index_set_type; wf_type psi0(index_set{0}, test_scf::h2_cmos()); auto H = test_scf::h2_hamiltonian(); chemist::braket::BraKet H_00(psi0, H, psi0); - const auto& [e, psi] = mod.run_as>(H_00, psi0); - REQUIRE_THAT(e, WithinAbs(-1.1167592336, 1E-6)); + const auto& [e, psi] = mod.template run_as>(H_00, psi0); + + pcorr->at() = -1.1167592336; + tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); + REQUIRE(approximately_equal(corr, e, 1E-6)); } \ No newline at end of file diff --git a/tests/cxx/integration_tests/guess/core.cpp b/tests/cxx/integration_tests/guess/core.cpp index 66a7766..7c12cf6 100644 --- a/tests/cxx/integration_tests/guess/core.cpp +++ b/tests/cxx/integration_tests/guess/core.cpp @@ -33,14 +33,12 @@ TEMPLATE_LIST_TEST_CASE("Core", "", test_scf::float_types) { typename rscf_wf::orbital_index_set_type occs{0}; REQUIRE(psi.orbital_indices() == occs); REQUIRE(psi.orbitals().from_space() == aos); - const auto& evals = psi.orbitals().diagonalized_matrix(); - using allocator_type = tensorwrapper::allocator::Eigen; - const auto& eval_buffer = allocator_type::rebind(evals.buffer()); + const auto& evals = psi.orbitals().diagonalized_matrix(); + tensorwrapper::allocator::Eigen alloc(mm.get_runtime()); + tensorwrapper::shape::Smooth shape_corr{2}; + auto pbuffer = alloc.construct({-1.25330893, -0.47506974}); + tensorwrapper::Tensor corr(shape_corr, std::move(pbuffer)); - const auto tol = 1E-6; - using Catch::Matchers::WithinAbs; - if constexpr(std::is_same_v) { - REQUIRE_THAT(eval_buffer.at(0), WithinAbs(-1.25330893, tol)); - REQUIRE_THAT(eval_buffer.at(1), WithinAbs(-0.47506974, tol)); - } + using tensorwrapper::operations::approximately_equal; + REQUIRE(approximately_equal(corr, evals, 1E-6)); } \ No newline at end of file diff --git a/tests/cxx/integration_tests/matrix_builder/determinant_driver.cpp b/tests/cxx/integration_tests/matrix_builder/determinant_driver.cpp index 5ab85e4..cc9a1cc 100644 --- a/tests/cxx/integration_tests/matrix_builder/determinant_driver.cpp +++ b/tests/cxx/integration_tests/matrix_builder/determinant_driver.cpp @@ -26,8 +26,8 @@ template using erased_type = chemist::braket::BraKet; -TEST_CASE("DeterminantDriver") { - using float_type = double; +TEMPLATE_LIST_TEST_CASE("DeterminantDriver", "", test_scf::float_types) { + using float_type = TestType; auto mm = test_scf::load_modules(); auto mod = mm.at("Determinant driver"); @@ -37,12 +37,20 @@ TEST_CASE("DeterminantDriver") { wf_type psi(index_set{0}, test_scf::h2_cmos()); simde::type::many_electrons es{2}; + tensorwrapper::allocator::Eigen alloc(mm.get_runtime()); + tensorwrapper::shape::Smooth shape_corr{}; + auto pcorr = alloc.allocate(tensorwrapper::layout::Physical(shape_corr)); + using tensorwrapper::operations::approximately_equal; + SECTION("Calling Kinetic") { simde::type::T_e_type T_e(es); chemist::braket::BraKet braket(psi, T_e, psi); erased_type copy_braket(braket); - const auto& T = mod.run_as>(copy_braket); - REQUIRE_THAT(T, WithinAbs(0.637692, 1E-6)); + const auto& T = mod.template run_as>(copy_braket); + + pcorr->at() = 0.637692; + tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); + REQUIRE(approximately_equal(corr, T, 1E-6)); } SECTION("Calling Electron-Nuclear Attraction") { @@ -50,8 +58,11 @@ TEST_CASE("DeterminantDriver") { simde::type::V_en_type V_en(es, h2_nuclei); chemist::braket::BraKet braket(psi, V_en, psi); erased_type copy_braket(braket); - const auto& V = mod.run_as>(copy_braket); - REQUIRE_THAT(V, WithinAbs(-1.96830777255516853, 1E-6)); + const auto& V = mod.template run_as>(copy_braket); + + pcorr->at() = -1.96830777255516853; + tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); + REQUIRE(approximately_equal(corr, V, 1E-6)); } SECTION("Calling J") { @@ -59,8 +70,11 @@ TEST_CASE("DeterminantDriver") { simde::type::J_e_type J_e(es, rho); chemist::braket::BraKet braket(psi, J_e, psi); erased_type copy_braket(braket); - const auto& J = mod.run_as>(copy_braket); - REQUIRE_THAT(J, WithinAbs(0.76056339681664897, 1E-6)); + const auto& J = mod.template run_as>(copy_braket); + + pcorr->at() = 0.76056339681664897; + tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); + REQUIRE(approximately_equal(corr, J, 1E-6)); } SECTION("Calling K") { @@ -68,7 +82,10 @@ TEST_CASE("DeterminantDriver") { simde::type::K_e_type K_e(es, rho); chemist::braket::BraKet braket(psi, K_e, psi); erased_type copy_braket(braket); - const auto& K = mod.run_as>(copy_braket); - REQUIRE_THAT(K, WithinAbs(0.76056339681664897, 1E-6)); + const auto& K = mod.template run_as>(copy_braket); + + pcorr->at() = 0.76056339681664897; + tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); + REQUIRE(approximately_equal(corr, K, 1E-6)); } } \ No newline at end of file diff --git a/tests/cxx/integration_tests/matrix_builder/electronic_energy.cpp b/tests/cxx/integration_tests/matrix_builder/electronic_energy.cpp index 7fd7e75..88b81ac 100644 --- a/tests/cxx/integration_tests/matrix_builder/electronic_energy.cpp +++ b/tests/cxx/integration_tests/matrix_builder/electronic_energy.cpp @@ -21,14 +21,19 @@ template using pt = simde::eval_braket; -TEST_CASE("ElectronicEnergy") { - using float_type = double; +TEMPLATE_LIST_TEST_CASE("ElectronicEnergy", "", test_scf::float_types) { + using float_type = TestType; auto mm = test_scf::load_modules(); auto mod = mm.at("Electronic energy"); using wf_type = simde::type::rscf_wf; using index_set = typename wf_type::orbital_index_set_type; + tensorwrapper::allocator::Eigen alloc(mm.get_runtime()); + tensorwrapper::shape::Smooth shape_corr{}; + auto pcorr = alloc.allocate(tensorwrapper::layout::Physical(shape_corr)); + using tensorwrapper::operations::approximately_equal; + wf_type psi(index_set{0}, test_scf::h2_cmos()); simde::type::many_electrons es{2}; @@ -44,6 +49,9 @@ TEST_CASE("ElectronicEnergy") { K_e); chemist::braket::BraKet braket(psi, H_e, psi); - const auto& E_elec = mod.run_as>(braket); - REQUIRE_THAT(E_elec, WithinAbs(-1.90066758625308307, 1E-6)); + const auto& E_elec = mod.template run_as>(braket); + + pcorr->at() = -1.90066758625308307; + tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); + REQUIRE(approximately_equal(corr, E_elec, 1E-6)); } \ No newline at end of file diff --git a/tests/cxx/integration_tests/matrix_builder/fock.cpp b/tests/cxx/integration_tests/matrix_builder/fock.cpp index 6054627..b228c7d 100644 --- a/tests/cxx/integration_tests/matrix_builder/fock.cpp +++ b/tests/cxx/integration_tests/matrix_builder/fock.cpp @@ -28,6 +28,12 @@ TEMPLATE_LIST_TEST_CASE("Fock Matrix Builder", "", test_scf::float_types) { auto& mod = mm.at("Fock Matrix Builder"); auto aos = test_scf::h2_aos(); + using tensorwrapper::operations::approximately_equal; + tensorwrapper::allocator::Eigen alloc(mm.get_runtime()); + tensorwrapper::shape::Smooth shape_corr{2, 2}; + tensorwrapper::layout::Physical l(shape_corr); + auto pcorr = alloc.allocate(l); + SECTION("No J or K") { auto h2 = test_scf::make_h2(); simde::type::electron e; @@ -38,15 +44,14 @@ TEMPLATE_LIST_TEST_CASE("Fock Matrix Builder", "", test_scf::float_types) { chemist::braket::BraKet f_mn(aos, f_e, aos); const auto& F = mod.template run_as(f_mn); - using alloc_type = tensorwrapper::allocator::Eigen; - const auto& F_eigen = alloc_type::rebind(F.buffer()); - using Catch::Matchers::WithinAbs; - if constexpr(std::is_same_v) { - REQUIRE_THAT(F_eigen.at(0, 0), WithinAbs(-1.120958, 1E-6)); - REQUIRE_THAT(F_eigen.at(0, 1), WithinAbs(-0.959374, 1E-6)); - REQUIRE_THAT(F_eigen.at(1, 0), WithinAbs(-0.959374, 1E-6)); - REQUIRE_THAT(F_eigen.at(1, 1), WithinAbs(-1.120958, 1E-6)); - } + pcorr->at(0, 0) = -1.120958; + pcorr->at(0, 1) = -0.959374; + pcorr->at(1, 0) = -0.959374; + pcorr->at(1, 1) = -1.120958; + + tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); + + REQUIRE(approximately_equal(F, corr, 1E-6)); } SECTION("With J and K") { @@ -54,15 +59,13 @@ TEMPLATE_LIST_TEST_CASE("Fock Matrix Builder", "", test_scf::float_types) { chemist::braket::BraKet f_mn(aos, f_e, aos); const auto& F = mod.template run_as(f_mn); - using alloc_type = tensorwrapper::allocator::Eigen; - const auto& F_eigen = alloc_type::rebind(F.buffer()); + pcorr->at(0, 0) = -0.319459; + pcorr->at(0, 1) = -0.571781; + pcorr->at(1, 0) = -0.571781; + pcorr->at(1, 1) = -0.319459; + + tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); - using Catch::Matchers::WithinAbs; - if constexpr(std::is_same_v) { - REQUIRE_THAT(F_eigen.at(0, 0), WithinAbs(-0.319459, 1E-6)); - REQUIRE_THAT(F_eigen.at(0, 1), WithinAbs(-0.571781, 1E-6)); - REQUIRE_THAT(F_eigen.at(1, 0), WithinAbs(-0.571781, 1E-6)); - REQUIRE_THAT(F_eigen.at(1, 1), WithinAbs(-0.319459, 1E-6)); - } + REQUIRE(approximately_equal(F, corr, 1E-6)); } } diff --git a/tests/cxx/integration_tests/matrix_builder/scf_integrals_driver.cpp b/tests/cxx/integration_tests/matrix_builder/scf_integrals_driver.cpp index 90dbcd7..d3a0b46 100644 --- a/tests/cxx/integration_tests/matrix_builder/scf_integrals_driver.cpp +++ b/tests/cxx/integration_tests/matrix_builder/scf_integrals_driver.cpp @@ -19,30 +19,6 @@ using pt = simde::aos_op_base_aos; using simde::type::tensor; -namespace { - -template -void compare_matrices(const tensor& A, const tensor& A_corr) { - using Catch::Matchers::WithinAbs; - using alloc_type = tensorwrapper::allocator::Eigen; - const auto& A_buffer = alloc_type::rebind(A.buffer()); - const auto& A_corr_buffer = alloc_type::rebind(A_corr.buffer()); - - const auto tol = 1E-6; - if constexpr(std::is_same_v) { - REQUIRE_THAT(A_buffer.at(0, 0), - WithinAbs(A_corr_buffer.at(0, 0), 1E-6)); - REQUIRE_THAT(A_buffer.at(0, 1), - WithinAbs(A_corr_buffer.at(0, 1), 1E-6)); - REQUIRE_THAT(A_buffer.at(1, 0), - WithinAbs(A_corr_buffer.at(1, 0), 1E-6)); - REQUIRE_THAT(A_buffer.at(1, 1), - WithinAbs(A_corr_buffer.at(1, 1), 1E-6)); - } -} - -} // namespace - using erased_type = chemist::braket::BraKet; @@ -55,6 +31,8 @@ TEMPLATE_LIST_TEST_CASE("SCFIntegralsDriver", "", test_scf::float_types) { simde::type::electron e; auto rho = test_scf::h2_density(); + using tensorwrapper::operations::approximately_equal; + SECTION("Calling Kinetic") { auto& tmod = mm.at("Kinetic"); simde::type::t_e_type t_e(e); @@ -62,7 +40,7 @@ TEMPLATE_LIST_TEST_CASE("SCFIntegralsDriver", "", test_scf::float_types) { erased_type copy_braket(braket); const auto& T = mod.template run_as(copy_braket); const auto& T_corr = tmod.template run_as(braket); - compare_matrices(T, T_corr); + REQUIRE(approximately_equal(T, T_corr, 1E-6)); } SECTION("Calling Electron-Nuclear Attraction") { @@ -73,7 +51,7 @@ TEMPLATE_LIST_TEST_CASE("SCFIntegralsDriver", "", test_scf::float_types) { erased_type copy_braket(braket); const auto& V = mod.template run_as(copy_braket); const auto& V_corr = tmod.template run_as(braket); - compare_matrices(V, V_corr); + REQUIRE(approximately_equal(V, V_corr, 1E-6)); } SECTION("Calling J Matrix") { @@ -83,7 +61,7 @@ TEMPLATE_LIST_TEST_CASE("SCFIntegralsDriver", "", test_scf::float_types) { erased_type copy_braket(braket); const auto& J = mod.template run_as(copy_braket); const auto& J_corr = jmod.template run_as(braket); - compare_matrices(J, J_corr); + REQUIRE(approximately_equal(J, J_corr, 1E-6)); } SECTION("Calling K Matrix") { @@ -93,7 +71,7 @@ TEMPLATE_LIST_TEST_CASE("SCFIntegralsDriver", "", test_scf::float_types) { erased_type copy_braket(braket); const auto& K = mod.template run_as(copy_braket); const auto& K_corr = kmod.template run_as(braket); - compare_matrices(K, K_corr); + REQUIRE(approximately_equal(K, K_corr, 1E-6)); } SECTION("Calling density matrix") { @@ -106,7 +84,7 @@ TEMPLATE_LIST_TEST_CASE("SCFIntegralsDriver", "", test_scf::float_types) { const auto& P = mod.template run_as(copy_braket); using op_pt = simde::aos_rho_e_aos; const auto& P_corr = pmod.template run_as(braket); - compare_matrices(P, P_corr); + REQUIRE(approximately_equal(P, P_corr, 1E-6)); } // SECTION("Calling Fock Matrix") { diff --git a/tests/cxx/integration_tests/update/diagonalization.cpp b/tests/cxx/integration_tests/update/diagonalization.cpp index 4ba00d3..7a10d3f 100644 --- a/tests/cxx/integration_tests/update/diagonalization.cpp +++ b/tests/cxx/integration_tests/update/diagonalization.cpp @@ -47,15 +47,13 @@ TEMPLATE_LIST_TEST_CASE("Diagaonalization", "", test_scf::float_types) { REQUIRE(psi.orbitals().from_space() == aos); // Check orbital energies - const auto& evals = psi.orbitals().diagonalized_matrix(); - using vector_allocator = tensorwrapper::allocator::Eigen; - const auto& eval_buffer = vector_allocator::rebind(evals.buffer()); - - const auto tol = 1E-6; - using Catch::Matchers::WithinAbs; - if constexpr(std::is_same_v) { - REQUIRE_THAT(eval_buffer.at(0), WithinAbs(-1.25330893, tol)); - REQUIRE_THAT(eval_buffer.at(1), WithinAbs(-0.47506974, tol)); - } + const auto& evals = psi.orbitals().diagonalized_matrix(); + tensorwrapper::allocator::Eigen alloc(mm.get_runtime()); + auto corr_buffer = alloc.construct({-1.25330893, -0.47506974}); + tensorwrapper::shape::Smooth corr_shape{2}; + tensorwrapper::Tensor corr(corr_shape, std::move(corr_buffer)); + + using tensorwrapper::operations::approximately_equal; + REQUIRE(approximately_equal(evals, corr, 1E-6)); } } \ No newline at end of file diff --git a/tests/cxx/test_scf.hpp b/tests/cxx/test_scf.hpp index fb7891b..08719cb 100644 --- a/tests/cxx/test_scf.hpp +++ b/tests/cxx/test_scf.hpp @@ -23,6 +23,8 @@ namespace test_scf { +using float_types = std::tuple; + /// Makes a H nucleus at the point @p x, @p y, @p z inline auto h_nucleus(double x, double y, double z) { return simde::type::nucleus("H", 1ul, 1836.15, x, y, z); diff --git a/tests/cxx/unit_tests/coulombs_law.cpp b/tests/cxx/unit_tests/coulombs_law.cpp index 80038bf..e144366 100644 --- a/tests/cxx/unit_tests/coulombs_law.cpp +++ b/tests/cxx/unit_tests/coulombs_law.cpp @@ -23,10 +23,16 @@ using pt = simde::charge_charge_interaction; using Catch::Matchers::WithinAbs; TEST_CASE("CoulombsLaw") { + using float_type = double; pluginplay::ModuleManager mm; scf::load_modules(mm); auto& mod = mm.at("Coulomb's Law"); + tensorwrapper::allocator::Eigen alloc(mm.get_runtime()); + tensorwrapper::shape::Smooth shape_corr{}; + auto pcorr = alloc.allocate(tensorwrapper::layout::Physical(shape_corr)); + using tensorwrapper::operations::approximately_equal; + using charge_type = typename simde::type::charges::value_type; charge_type q0(1.0, 2.0, 3.0, 4.0); charge_type q1(-1.0, 3.0, 4.0, 5.0); @@ -41,19 +47,25 @@ TEST_CASE("CoulombsLaw") { } SECTION("charges w/ itself") { - auto e = mod.run_as(qs, qs); - REQUIRE_THAT(e, WithinAbs(-1.0103629710818451, 1E-6)); + auto e = mod.run_as(qs, qs); + pcorr->at() = -1.0103629710818451; + tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); + REQUIRE(approximately_equal(corr, e, 1E-6)); } SECTION("charges w/ empty") { - auto e = mod.run_as(qs, empty); - REQUIRE_THAT(e, WithinAbs(0.0, 1E-6)); + auto e = mod.run_as(qs, empty); + pcorr->at() = 0.0; + tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); + REQUIRE(approximately_equal(corr, e, 1E-6)); } SECTION("charges w/ different charges") { simde::type::charges qs0{q0}; simde::type::charges qs12{q1, q2}; - auto e = mod.run_as(qs0, qs12); - REQUIRE_THAT(e, WithinAbs(-0.1443375672974065, 1E-6)); + auto e = mod.run_as(qs0, qs12); + pcorr->at() = -0.1443375672974065; + tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); + REQUIRE(approximately_equal(corr, e, 1E-6)); } } \ No newline at end of file diff --git a/tests/cxx/unit_tests/eigen_solver/eigen_generalized.cpp b/tests/cxx/unit_tests/eigen_solver/eigen_generalized.cpp index 577cdff..fcbcbad 100644 --- a/tests/cxx/unit_tests/eigen_solver/eigen_generalized.cpp +++ b/tests/cxx/unit_tests/eigen_solver/eigen_generalized.cpp @@ -18,7 +18,8 @@ #include #include -TEST_CASE("EigenGeneralized") { +TEMPLATE_LIST_TEST_CASE("EigenGeneralized", "", test_scf::float_types) { + using float_type = TestType; pluginplay::ModuleManager mm; scf::load_modules(mm); @@ -26,15 +27,31 @@ TEST_CASE("EigenGeneralized") { auto& mod = mm.at("Generalized eigensolve via Eigen"); - simde::type::tensor A({{1.0, 2.0}, {2.0, 3.0}}); - simde::type::tensor B({{1.0, 0.0}, {0.0, 1.0}}); + tensorwrapper::allocator::Eigen alloc(mm.get_runtime()); + tensorwrapper::shape::Smooth shape{2, 2}; + tensorwrapper::layout::Physical l(shape); + auto A_buffer = alloc.allocate(l); + A_buffer->at(0, 0) = 1.0; + A_buffer->at(0, 1) = 2.0; + A_buffer->at(1, 0) = 2.0; + A_buffer->at(1, 1) = 3.0; + + auto B_buffer = alloc.allocate(l); + B_buffer->at(0, 0) = 1.0; + B_buffer->at(0, 1) = 0.0; + B_buffer->at(1, 0) = 0.0; + B_buffer->at(1, 1) = 1.0; + + simde::type::tensor A(shape, std::move(A_buffer)); + simde::type::tensor B(shape, std::move(B_buffer)); auto&& [values, vector] = mod.run_as(A, B); - using value_alloc_t = tensorwrapper::allocator::Eigen; - const auto& eigen_values = value_alloc_t::rebind(values.buffer()); + const auto& eigen_values = alloc.rebind(values.buffer()); + + auto corr_buffer = alloc.construct({-0.236068, 4.236068}); + tensorwrapper::shape::Smooth corr_shape{2}; + simde::type::tensor corr(corr_shape, std::move(corr_buffer)); - using Catch::Matchers::WithinAbs; - REQUIRE_THAT(eigen_values.at(0), WithinAbs(-0.236068, 1E-6)); - REQUIRE_THAT(eigen_values.at(1), WithinAbs(4.236068, 1E-6)); + REQUIRE(tensorwrapper::operations::approximately_equal(corr, values, 1E-6)); } \ No newline at end of file diff --git a/tests/cxx/unit_tests/fock_operator/restricted.cpp b/tests/cxx/unit_tests/fock_operator/restricted.cpp index 7d0a522..45cdcf1 100644 --- a/tests/cxx/unit_tests/fock_operator/restricted.cpp +++ b/tests/cxx/unit_tests/fock_operator/restricted.cpp @@ -18,7 +18,7 @@ #include #include -TEST_CASE("Restricted") { +TEMPLATE_LIST_TEST_CASE("Restricted", "", test_scf::float_types) { pluginplay::ModuleManager mm; scf::load_modules(mm); using float_type = double; diff --git a/tests/cxx/unit_tests/matrix_builder/density_matrix.cpp b/tests/cxx/unit_tests/matrix_builder/density_matrix.cpp index cbe232b..1b5b558 100644 --- a/tests/cxx/unit_tests/matrix_builder/density_matrix.cpp +++ b/tests/cxx/unit_tests/matrix_builder/density_matrix.cpp @@ -18,8 +18,8 @@ using pt = simde::aos_rho_e_aos; -TEST_CASE("Density Matrix Builder") { - using float_type = double; +TEMPLATE_LIST_TEST_CASE("Density Matrix Builder", "", test_scf::float_types) { + using float_type = TestType; pluginplay::ModuleManager mm; scf::load_modules(mm); auto& mod = mm.at("Density matrix builder"); @@ -29,13 +29,13 @@ TEST_CASE("Density Matrix Builder") { simde::type::rho_e rho_hat(cmos, occs); chemist::braket::BraKet p_mn(aos, rho_hat, aos); - const auto& P = mod.run_as(p_mn); - using allocator_type = tensorwrapper::allocator::Eigen; - const auto& P_eigen = allocator_type::rebind(P.buffer()); + const auto& P = mod.run_as(p_mn); + tensorwrapper::allocator::Eigen alloc(mm.get_runtime()); + tensorwrapper::shape::Smooth corr_shape{2, 2}; + tensorwrapper::layout::Physical l(corr_shape); + auto corr_buffer = alloc.construct(l, 0.31980835); + tensorwrapper::Tensor corr(corr_shape, std::move(corr_buffer)); - using Catch::Matchers::WithinAbs; - REQUIRE_THAT(P_eigen.at(0, 0), WithinAbs(0.31980835, 1E-6)); - REQUIRE_THAT(P_eigen.at(0, 1), WithinAbs(0.31980835, 1E-6)); - REQUIRE_THAT(P_eigen.at(1, 0), WithinAbs(0.31980835, 1E-6)); - REQUIRE_THAT(P_eigen.at(1, 1), WithinAbs(0.31980835, 1E-6)); + using tensorwrapper::operations::approximately_equal; + REQUIRE(approximately_equal(P, corr, 1E-6)); } \ No newline at end of file From 6ef373f372ea4b37862a6397e9a4d862d6752681 Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Thu, 6 Mar 2025 11:21:01 -0600 Subject: [PATCH 4/4] swap to approx equal --- tests/cxx/unit_tests/coulombs_law.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/tests/cxx/unit_tests/coulombs_law.cpp b/tests/cxx/unit_tests/coulombs_law.cpp index e144366..9300d26 100644 --- a/tests/cxx/unit_tests/coulombs_law.cpp +++ b/tests/cxx/unit_tests/coulombs_law.cpp @@ -42,8 +42,10 @@ TEST_CASE("CoulombsLaw") { simde::type::charges qs{q0, q1, q2}; SECTION("empty points") { - auto e = mod.run_as(empty, empty); - REQUIRE(e == 0.0); + auto e = mod.run_as(empty, empty); + pcorr->at() = 0.0; + tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); + REQUIRE(approximately_equal(corr, e, 1E-6)); } SECTION("charges w/ itself") {