g2g seems to work

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}")

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 <gau2grid/gau2grid.h>
+#include <simde/integration_grids/collocation_matrix.hpp>
+
+namespace scf::xc {
+namespace {
+const auto desc = R"(
+
+CollocationMatrix
+-----------------
+)";
+
+template<typename T>
+std::vector<T> flatten_grid(const chemist::Grid& grid) {
+    std::vector<T> flattened_grid;
+    flattened_grid.reserve(grid.size() * 3);
+    for(const auto& point : grid) {
+        flattened_grid.push_back(static_cast<T>(point.point().x()));
+        flattened_grid.push_back(static_cast<T>(point.point().y()));
+        flattened_grid.push_back(static_cast<T>(point.point().z()));
+    }
+    return flattened_grid;
+}
+
+} // namespace
+
+using pt = simde::CollocationMatrix;
+
+MODULE_CTOR(Gau2Grid) {
+    satisfies_property_type<pt>();
+
+    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<float_type>(grid);
+
+    tensorwrapper::allocator::Eigen<float_type> 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<double> exponents(n_primitives);
+            std::vector<double> coefficients(n_primitives);
+            std::vector<double> 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

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>("Gau2Grid");
     mm.add_module<GridFromFile>("Grid From File");
 }
 

src/scf/xc/xc.hpp

@@ -19,6 +19,7 @@
 
 namespace scf::xc {
 
+DECLARE_MODULE(Gau2Grid);
 DECLARE_MODULE(GridFromFile);
 
 void set_defaults(pluginplay::ModuleManager& mm);

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 <iostream>
+#include <pluginplay/pluginplay.hpp>
+
+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<float_type>;
+    using cg_type        = chemist::basis_set::ContractedGaussian<prim_type>;
+    using shell_type     = chemist::basis_set::Shell<cg_type>;
+    using atomic_bs_type = chemist::basis_set::AtomicBasisSet<shell_type>;
+    using ao_basis_type  = chemist::basis_set::AOBasisSet<atomic_bs_type>;
+
+    // TODO: better testing 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<chemist::GridPoint> grid_points;
+        for(std::size_t i = 0; i < 5; ++i)
+            grid_points.push_back({1.0, 0.0, 0.0, static_cast<double>(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<float_type> coefs{1.0};
+        std::vector<float_type> 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<pt>(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<pt>(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));
+        };
+    }
+}