From 212fe4fbe3e29f78c03e06a3fe3fe8dd1b332c83 Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Tue, 23 Sep 2025 22:26:52 -0500 Subject: [PATCH 1/4] g2g seems to work --- CMakeLists.txt | 12 +++- src/scf/xc/gau2grid.cpp | 104 +++++++++++++++++++++++++++ src/scf/xc/xc.cpp | 1 + src/scf/xc/xc.hpp | 1 + tests/cxx/unit_tests/xc/gau2grid.cpp | 86 ++++++++++++++++++++++ 5 files changed, 203 insertions(+), 1 deletion(-) create mode 100644 src/scf/xc/gau2grid.cpp create mode 100644 tests/cxx/unit_tests/xc/gau2grid.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 9234844..d7128d3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -59,11 +59,21 @@ cmaize_find_or_build_dependency( CMAKE_ARGS EIGEN_BUILD_TESTING=OFF ) +cmaize_find_or_build_dependency( + gau2grid + NAME Gau2grid + URL https://github.com/psi4/gau2grid + VERSION master + BUILD_TARGET gg + FIND_TARGET gg + +) + cmaize_add_library( scf SOURCE_DIR "${CMAKE_CURRENT_LIST_DIR}/src/scf" INCLUDE_DIRS "${CMAKE_CURRENT_LIST_DIR}/include/scf" - DEPENDS "${DEPENDENCIES}" eigen + DEPENDS "${DEPENDENCIES}" eigen gau2grid ) if("${BUILD_TAMM_SCF}") diff --git a/src/scf/xc/gau2grid.cpp b/src/scf/xc/gau2grid.cpp new file mode 100644 index 0000000..389287a --- /dev/null +++ b/src/scf/xc/gau2grid.cpp @@ -0,0 +1,104 @@ +/* + * Copyright 2025 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include "xc.hpp" +#include +#include + +namespace scf::xc { +namespace { +const auto desc = R"( + +CollocationMatrix +----------------- +)"; + +template +std::vector flatten_grid(const chemist::Grid& grid) { + std::vector flattened_grid; + flattened_grid.reserve(grid.size() * 3); + for(const auto& point : grid) { + flattened_grid.push_back(static_cast(point.point().x())); + flattened_grid.push_back(static_cast(point.point().y())); + flattened_grid.push_back(static_cast(point.point().z())); + } + return flattened_grid; +} + +} // namespace + +using pt = simde::CollocationMatrix; + +MODULE_CTOR(Gau2Grid) { + satisfies_property_type(); + + description(desc); +} + +MODULE_RUN(Gau2Grid) { + const auto& [grid, ao_basis] = pt::unwrap_inputs(inputs); + auto n_points = grid.size(); + + using float_type = double; + auto flattened_grid = flatten_grid(grid); + + tensorwrapper::allocator::Eigen allocator(get_runtime()); + tensorwrapper::shape::Smooth matrix_shape{ao_basis.n_aos(), n_points}; + tensorwrapper::layout::Physical layout(matrix_shape); + auto matrix_buffer = allocator.allocate(layout); + auto matrix_data = matrix_buffer->get_mutable_data(); + + std::size_t ao_i = 0; + for(const auto& atomic_basis : ao_basis) { + for(const auto& shell_i : atomic_basis) { + const auto& cg = shell_i.contracted_gaussian(); + const auto L = shell_i.l(); + const auto n_primitives = cg.size(); + const auto n_aos = shell_i.size(); + + // TODO: Expose exponent_data/coefficient_data methods for Shells + std::vector exponents(n_primitives); + std::vector coefficients(n_primitives); + std::vector center(3); + for(std::size_t i = 0; i < n_primitives; ++i) { + exponents[i] = cg.at(i).exponent(); + coefficients[i] = cg.at(i).coefficient(); + } + center[0] = shell_i.center().x(); + center[1] = shell_i.center().y(); + center[2] = shell_i.center().z(); + + auto is_pure = shell_i.pure() == chemist::ShellType::pure; + auto order = is_pure ? GG_SPHERICAL_CCA : GG_CARTESIAN_CCA; + + auto offset = ao_i * n_points; + + auto shell_i_data = matrix_data + offset; + gg_collocation(L, n_points, flattened_grid.data(), 3, n_primitives, + coefficients.data(), exponents.data(), center.data(), + order, shell_i_data); + ao_i += n_aos; + } + } + + simde::type::tensor collocation_matrix(matrix_shape, + std::move(matrix_buffer)); + + auto rv = results(); + return pt::wrap_results(rv, std::move(collocation_matrix)); +} + +} // namespace scf::xc diff --git a/src/scf/xc/xc.cpp b/src/scf/xc/xc.cpp index d09bf11..3378ffe 100644 --- a/src/scf/xc/xc.cpp +++ b/src/scf/xc/xc.cpp @@ -20,6 +20,7 @@ namespace scf::xc { void load_modules(pluginplay::ModuleManager& mm) { gauxc::load_modules(mm); + mm.add_module("Gau2Grid"); mm.add_module("Grid From File"); } diff --git a/src/scf/xc/xc.hpp b/src/scf/xc/xc.hpp index 9eb684e..275b806 100644 --- a/src/scf/xc/xc.hpp +++ b/src/scf/xc/xc.hpp @@ -19,6 +19,7 @@ namespace scf::xc { +DECLARE_MODULE(Gau2Grid); DECLARE_MODULE(GridFromFile); void set_defaults(pluginplay::ModuleManager& mm); diff --git a/tests/cxx/unit_tests/xc/gau2grid.cpp b/tests/cxx/unit_tests/xc/gau2grid.cpp new file mode 100644 index 0000000..d09f518 --- /dev/null +++ b/tests/cxx/unit_tests/xc/gau2grid.cpp @@ -0,0 +1,86 @@ +/* + * Copyright 2025 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include "../../test_scf.hpp" +#include +#include + +using namespace scf; + +using pt = simde::CollocationMatrix; + +using tensorwrapper::operations::approximately_equal; + +TEST_CASE("Gau2Grid") { + pluginplay::ModuleManager mm; + scf::load_modules(mm); + auto& mod = mm.at("Gau2Grid"); + + using float_type = double; + using prim_type = chemist::basis_set::Primitive; + using cg_type = chemist::basis_set::ContractedGaussian; + using shell_type = chemist::basis_set::Shell; + using atomic_bs_type = chemist::basis_set::AtomicBasisSet; + using ao_basis_type = chemist::basis_set::AOBasisSet; + + auto cart = chemist::ShellType::cartesian; + auto pure = chemist::ShellType::pure; + + SECTION("Manual examples") { + // Taken from https://gau2grid.readthedocs.io/en/latest/c_example.html + + std::vector grid_points; + for(std::size_t i = 0; i < 5; ++i) + grid_points.push_back({1.0, 0.0, 0.0, static_cast(i)}); + chemist::Grid grid(grid_points.begin(), grid_points.end()); + + ao_basis_type ao_basis; + atomic_bs_type abs("n/a", 1, {0.0, 0.0, 0.0}); + + std::vector coefs{1.0}; + std::vector exps{1.0}; + cg_type cg(coefs.begin(), coefs.end(), exps.begin(), exps.end(), 0.0, + 0.0, 0.0); + + SECTION("Single Basis Function example") { + abs.add_shell(pure, 0, cg); + ao_basis.add_center(std::move(abs)); + auto rv = mod.run_as(grid, ao_basis); + simde::type::tensor corr{ + {1.0, 0.367879, 0.0183156, 0.00012341, 0.0}}; + REQUIRE(approximately_equal(rv, corr, 1e-6)); + } + + SECTION("Multiple Basis Function example") { + abs.add_shell(pure, 0, cg); + abs.add_shell(pure, 1, cg); + abs.add_shell(pure, 2, cg); + ao_basis.add_center(std::move(abs)); + auto rv = mod.run_as(grid, ao_basis); + simde::type::tensor corr{ + {1, 0.367879, 0.0183156, 0.00012341, 1.12535e-07}, + {0, 0, 0, 0, 0}, + {0, 0.367879, 0.0366313, 0.000370229, 4.50141e-07}, + {0, 0, 0, 0, 0}, + {0, 0, 0, 0, 0}, + {0, 0, 0, 0, 0}, + {0, 0.367879, 0.0732626, 0.00111069, 1.80056e-06}, + {0, 0, 0, 0, 0}, + {0, 0, 0, 0, 0}}; + REQUIRE(approximately_equal(rv, corr, 1e-6)); + }; + } +} From c839021af47137d1433f4a68e6c7087867b3d694 Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Wed, 24 Sep 2025 09:34:48 -0500 Subject: [PATCH 2/4] backup --- src/scf/xc/gau2grid.cpp | 2 +- tests/cxx/unit_tests/xc/gau2grid.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/scf/xc/gau2grid.cpp b/src/scf/xc/gau2grid.cpp index 389287a..381eba2 100644 --- a/src/scf/xc/gau2grid.cpp +++ b/src/scf/xc/gau2grid.cpp @@ -40,7 +40,7 @@ std::vector flatten_grid(const chemist::Grid& grid) { } // namespace -using pt = simde::CollocationMatrix; +using pt = simde::AOCollocationMatrix; MODULE_CTOR(Gau2Grid) { satisfies_property_type(); diff --git a/tests/cxx/unit_tests/xc/gau2grid.cpp b/tests/cxx/unit_tests/xc/gau2grid.cpp index d09f518..615adaf 100644 --- a/tests/cxx/unit_tests/xc/gau2grid.cpp +++ b/tests/cxx/unit_tests/xc/gau2grid.cpp @@ -20,7 +20,7 @@ using namespace scf; -using pt = simde::CollocationMatrix; +using pt = simde::AOCollocationMatrix; using tensorwrapper::operations::approximately_equal; From 2bc90b893e93e7ddba6af178480cb79b471b27d0 Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Fri, 26 Sep 2025 11:15:48 -0500 Subject: [PATCH 3/4] density2grid works (I think) --- src/scf/xc/density2grid.cpp | 98 ++++++++++++++++++++++++ src/scf/xc/xc.cpp | 1 + src/scf/xc/xc.hpp | 1 + tests/cxx/unit_tests/xc/density2grid.cpp | 48 ++++++++++++ 4 files changed, 148 insertions(+) create mode 100644 src/scf/xc/density2grid.cpp create mode 100644 tests/cxx/unit_tests/xc/density2grid.cpp diff --git a/src/scf/xc/density2grid.cpp b/src/scf/xc/density2grid.cpp new file mode 100644 index 0000000..9f99893 --- /dev/null +++ b/src/scf/xc/density2grid.cpp @@ -0,0 +1,98 @@ +/* + * Copyright 2025 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include "xc.hpp" +#include + +namespace scf::xc { +namespace { +const auto desc = R"( + +DensityCollocationMatrix +----------------- +)"; + +struct Kernel { + using buffer_base = tensorwrapper::buffer::BufferBase; + + template + auto run(const buffer_base& aos_on_grid, const buffer_base& X, + parallelzone::runtime::RuntimeView& rv) { + tensorwrapper::allocator::Eigen allocator(rv); + + const auto& eigen_aos_on_grid = allocator.rebind(aos_on_grid); + const auto* paos_on_grid = eigen_aos_on_grid.get_immutable_data(); + const auto& eigen_X = allocator.rebind(X); + const auto* pX = eigen_X.get_immutable_data(); + const auto& shape_X = eigen_X.layout().shape().as_smooth(); + auto n_aos = shape_X.extent(0); + auto n_grid = shape_X.extent(1); + + tensorwrapper::shape::Smooth rv_shape{n_grid}; + tensorwrapper::layout::Physical rv_layout(rv_shape); + auto rv_buffer = allocator.allocate(rv_layout); + + // AOs on rows, grid points on columns + for(std::size_t grid_i = 0; grid_i < n_grid; ++grid_i) { + FloatType sum = 0; + for(std::size_t ao_i = 0; ao_i < n_aos; ++ao_i) { + const auto idx = ao_i * n_grid + grid_i; + sum += paos_on_grid[idx] * pX[idx]; + } + rv_buffer->set_elem(std::vector{grid_i}, sum); + } + return simde::type::tensor(rv_shape, std::move(rv_buffer)); + } +}; +} // namespace + +using pt = simde::EDensityCollocationMatrix; +using ao2grid_pt = simde::AOCollocationMatrix; + +MODULE_CTOR(Density2Grid) { + satisfies_property_type(); + description(desc); + + add_submodule("AOs on a grid"); +} + +MODULE_RUN(Density2Grid) { + const auto& [grid, density] = pt::unwrap_inputs(inputs); + + const auto& rho = density.value(); + const auto& aos = density.basis_set().ao_basis_set(); + + auto& ao2grid_mod = submods.at("AOs on a grid"); + auto aos_on_grid = ao2grid_mod.run_as(grid, aos); + + auto n_grid = grid.size(); + using shape_type = tensorwrapper::shape::Smooth; + simde::type::tensor X; + X("m,i") = rho("m,n") * aos_on_grid("n,i"); + + using tensorwrapper::utilities::floating_point_dispatch; + Kernel k; + auto& runtime = get_runtime(); + const auto& aos_buffer = aos_on_grid.buffer(); + const auto& X_buffer = X.buffer(); + auto rho_on_grid = + floating_point_dispatch(k, aos_buffer, X_buffer, runtime); + + auto rv = results(); + return pt::wrap_results(rv, std::move(rho_on_grid)); +} + +} // namespace scf::xc diff --git a/src/scf/xc/xc.cpp b/src/scf/xc/xc.cpp index 3378ffe..69c9616 100644 --- a/src/scf/xc/xc.cpp +++ b/src/scf/xc/xc.cpp @@ -22,6 +22,7 @@ void load_modules(pluginplay::ModuleManager& mm) { gauxc::load_modules(mm); mm.add_module("Gau2Grid"); mm.add_module("Grid From File"); + mm.add_module("Density2Grid"); } void set_defaults(pluginplay::ModuleManager& mm) { gauxc::set_defaults(mm); } diff --git a/src/scf/xc/xc.hpp b/src/scf/xc/xc.hpp index 275b806..5fcf88d 100644 --- a/src/scf/xc/xc.hpp +++ b/src/scf/xc/xc.hpp @@ -19,6 +19,7 @@ namespace scf::xc { +DECLARE_MODULE(Density2Grid); DECLARE_MODULE(Gau2Grid); DECLARE_MODULE(GridFromFile); diff --git a/tests/cxx/unit_tests/xc/density2grid.cpp b/tests/cxx/unit_tests/xc/density2grid.cpp new file mode 100644 index 0000000..36be9a1 --- /dev/null +++ b/tests/cxx/unit_tests/xc/density2grid.cpp @@ -0,0 +1,48 @@ +/* + * Copyright 2025 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include "../../test_scf.hpp" +#include +#include + +using namespace scf; + +using pt = simde::EDensityCollocationMatrix; +using ao_pt = simde::AOCollocationMatrix; +using tensorwrapper::operations::approximately_equal; + +TEST_CASE("Density2Grid") { + pluginplay::ModuleManager mm; + scf::load_modules(mm); + auto& mod = mm.at("Density2Grid"); + + chemist::Grid grid; // Value doesn't matter for this test b/c of submodule + + SECTION("H2") { + auto rho = test_scf::h2_density(); + simde::type::tensor ao2grid{{1.0, 2.0, 3.0, 4.0}, {5.0, 6.0, 7.0, 8.0}}; + auto ao_mod = + pluginplay::make_lambda([&](auto&& grid_in, auto&& aos) { + REQUIRE(grid_in == grid); + REQUIRE(aos == rho.basis_set().ao_basis_set()); + return ao2grid; + }); + mod.change_submod("AOs on a grid", ao_mod); + auto rv = mod.run_as(grid, rho); + simde::type::tensor corr{11.513088, 20.467712, 31.9808, 46.052352}; + REQUIRE(approximately_equal(rv, corr, 1e-4)); + } +} From 77f24a1d0cb27bf617cca808a05972973827b4ba Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Fri, 26 Sep 2025 11:34:46 -0500 Subject: [PATCH 4/4] fix warnings --- src/scf/xc/density2grid.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/scf/xc/density2grid.cpp b/src/scf/xc/density2grid.cpp index 9f99893..ff3c632 100644 --- a/src/scf/xc/density2grid.cpp +++ b/src/scf/xc/density2grid.cpp @@ -78,8 +78,6 @@ MODULE_RUN(Density2Grid) { auto& ao2grid_mod = submods.at("AOs on a grid"); auto aos_on_grid = ao2grid_mod.run_as(grid, aos); - auto n_grid = grid.size(); - using shape_type = tensorwrapper::shape::Smooth; simde::type::tensor X; X("m,i") = rho("m,n") * aos_on_grid("n,i");