Adds building E_xc and V_xc via LibXC (#50)

* g2g seems to work

* backup

* density2grid works (I think)

* fix warnings

* it freaking works!!!!

* vxc done too!!!

* fix warnings

* fix warning

* Update CMakeLists.txt

Co-authored-by: Jonathan M. Waldrop <relik107@gmail.com>

* make libxc optional

* jonathan's suggestions

* formatting

---------

Co-authored-by: Jonathan M. Waldrop <relik107@gmail.com>

Author: ryanmrichard

CMakeLists.txt

@@ -32,12 +32,19 @@ set(SCF_TESTS_DIR "${CMAKE_CURRENT_SOURCE_DIR}/tests")
 nwx_cxx_api_docs("${SCF_SOURCE_DIR}" "${SCF_INCLUDE_DIR}")
 
 ### Options ###
+
 cmaize_option_list(
     BUILD_TESTING OFF "Should we build the tests?"
     BUILD_PYBIND11_PYBINDINGS ON "Build Python bindings with pybind11?"
     BUILD_TAMM_SCF OFF "Should we build modules that rely on TAMM/Exachem?"
+    BUILD_LIBXC ON "Should we build libxc?"
     INTEGRATION_TESTING OFF "Should we build the integration tests?"
 )
+
+if(CMAKE_CXX_COMPILER_ID STREQUAL "Clang")
+    set(BUILD_LIBXC OFF)
+endif()
+
 if("${BUILD_TAMM_SCF}")
     set(DEPENDENCIES simde gauxc tamm exachem chemcache)
     include(get_libint2)
@@ -69,11 +76,25 @@ cmaize_find_or_build_dependency(
 
 )
 
+set(BUILD_TESTING_BACKUP ${BUILD_TESTING})
+set(BUILD_TESTING OFF)
+cmaize_find_or_build_optional_dependency(
+    libxc
+    BUILD_LIBXC
+    NAME libxc
+    URL https://gitlab.com/libxc/libxc
+    VERSION devel
+    BUILD_TARGET xc
+    FIND_TARGET xc
+    CMAKE_ARGS BUILD_TESTING=OFF
+)
+set(BUILD_TESTING ${BUILD_TESTING_BACKUP})
+
 cmaize_add_library(
     scf
     SOURCE_DIR "${CMAKE_CURRENT_LIST_DIR}/src/scf"
     INCLUDE_DIRS "${CMAKE_CURRENT_LIST_DIR}/include/scf"
-    DEPENDS "${DEPENDENCIES}" eigen gau2grid
+    DEPENDS "${DEPENDENCIES}" eigen gau2grid libxc
 )
 
 if("${BUILD_TAMM_SCF}")

src/scf/xc/aos_on_grid.cpp

@@ -0,0 +1,73 @@
+/*
+ * 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 <simde/integration_grids/collocation_matrix.hpp>
+
+namespace scf::xc {
+namespace {
+const auto desc = R"(
+
+AO CollocationMatrix
+-----------------
+)";
+
+} // namespace
+
+using pt = simde::AOCollocationMatrix;
+
+MODULE_CTOR(AOsOnGrid) {
+    satisfies_property_type<pt>();
+    description(desc);
+}
+
+MODULE_RUN(AOsOnGrid) {
+    const auto& [grid, ao_basis] = pt::unwrap_inputs(inputs);
+    auto n_points                = grid.size();
+    auto n_aos                   = ao_basis.n_aos();
+    using float_type             = double;
+
+    tensorwrapper::allocator::Eigen<float_type> allocator(get_runtime());
+    tensorwrapper::shape::Smooth matrix_shape{n_aos, n_points};
+    tensorwrapper::layout::Physical layout(matrix_shape);
+    auto buffer = allocator.allocate(layout);
+
+    std::vector<std::size_t> idx{0, 0};
+    for(const auto& atomic_basis : ao_basis) {
+        for(const auto& shell_i : atomic_basis) {
+            assert(shell_i.l() == 0); // only s is supported for now
+            const auto& cg = shell_i.contracted_gaussian();
+            for(; idx[1] < n_points; ++idx[1]) {
+                float_type ao_value = 0.0;
+                for(const auto& prim : cg) {
+                    // TODO: update when eval accounts for normalization
+                    const auto val = prim.evaluate(grid.at(idx[1]).point());
+                    const auto exponent = prim.exponent();
+                    auto norm = std::sqrt(std::pow(2.0 * exponent / M_PI, 1.5));
+                    ao_value += norm * val;
+                }
+                buffer->set_elem(idx, ao_value);
+            }
+            idx[0] += shell_i.size();
+            idx[1] = 0;
+        }
+    }
+    simde::type::tensor collocation_matrix(matrix_shape, std::move(buffer));
+    auto rv = results();
+    return pt::wrap_results(rv, std::move(collocation_matrix));
+};
+
+} // namespace scf::xc

src/scf/xc/density2grid.cpp

@@ -14,6 +14,7 @@
  * limitations under the License.
  */
 
+#include "libxc/libxc.hpp"
 #include "xc.hpp"
 #include <simde/integration_grids/collocation_matrix.hpp>
 
@@ -25,38 +26,6 @@ DensityCollocationMatrix
 -----------------
 )";
 
-struct Kernel {
-    using buffer_base = tensorwrapper::buffer::BufferBase;
-
-    template<typename FloatType>
-    auto run(const buffer_base& aos_on_grid, const buffer_base& X,
-             parallelzone::runtime::RuntimeView& rv) {
-        tensorwrapper::allocator::Eigen<FloatType> 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;
@@ -78,16 +47,11 @@ MODULE_RUN(Density2Grid) {
     auto& ao2grid_mod = submods.at("AOs on a grid");
     auto aos_on_grid  = ao2grid_mod.run_as<ao2grid_pt>(grid, aos);
 
-    simde::type::tensor X;
-    X("m,i") = rho("m,n") * aos_on_grid("n,i");
+    simde::type::tensor X, rho2;
+    rho2("m,n") = rho("m,n") * 2.0; // Assumes restricted orbitals
+    X("m,i")    = rho2("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 rho_on_grid = libxc::batched_dot(aos_on_grid, X);
 
     auto rv = results();
     return pt::wrap_results(rv, std::move(rho_on_grid));

src/scf/xc/gau2grid.cpp

@@ -24,6 +24,9 @@ const auto desc = R"(
 
 CollocationMatrix
 -----------------
+
+Warning! Gau2Grid assumes that the primitive normalization has been folded into
+the contraction coefficients!!!
 )";
 
 template<typename T>
@@ -73,6 +76,7 @@ MODULE_RUN(Gau2Grid) {
             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();
@@ -84,12 +88,14 @@ MODULE_RUN(Gau2Grid) {
             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 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);
+            gg_collocation(static_cast<int>(L), static_cast<int>(n_points),
+                           flattened_grid.data(), 3,
+                           static_cast<int>(n_primitives), coefficients.data(),
+                           exponents.data(), center.data(), order,
+                           shell_i_data);
+
             ao_i += n_aos;
         }
     }