Density Fitting Modules

src/integrals/ao_integrals/ao_integrals.hpp

@@ -22,17 +22,32 @@ namespace integrals::ao_integrals {
 DECLARE_MODULE(AOIntegralsDriver);
 DECLARE_MODULE(JFourCenter);
 DECLARE_MODULE(KFourCenter);
+DECLARE_MODULE(JDensityFitted);
+DECLARE_MODULE(KDensityFitted);
+DECLARE_MODULE(DFIntegral);
+DECLARE_MODULE(CoulombMetric);
 
 inline void set_defaults(pluginplay::ModuleManager& mm) {
-    const auto ao_driver = "AO integral driver";
-    mm.change_submod(ao_driver, "Coulomb matrix", "Four center J builder");
-    mm.change_submod(ao_driver, "Exchange matrix", "Four center K builder");
+    mm.change_submod("AO integral driver", "Coulomb matrix",
+                     "Four center J builder");
+    mm.change_submod("AO integral driver", "Exchange matrix",
+                     "Four center K builder");
+    mm.change_submod("Density Fitted J builder", "DF ERI",
+                     "Density Fitting Integral");
+    mm.change_submod("Density Fitted K builder", "DF ERI",
+                     "Density Fitting Integral");
+    mm.change_submod("Density Fitting Integral", "Coulomb Metric",
+                     "Coulomb Metric");
 }
 
 inline void load_modules(pluginplay::ModuleManager& mm) {
     mm.add_module<AOIntegralsDriver>("AO integral driver");
     mm.add_module<JFourCenter>("Four center J builder");
     mm.add_module<KFourCenter>("Four center K builder");
+    mm.add_module<JDensityFitted>("Density Fitted J builder");
+    mm.add_module<KDensityFitted>("Density Fitted K builder");
+    mm.add_module<DFIntegral>("Density Fitting Integral");
+    mm.add_module<CoulombMetric>("Coulomb Metric");
     set_defaults(mm);
 }
 

src/integrals/ao_integrals/coulomb_metric.cpp

@@ -0,0 +1,74 @@
+/*
+ * 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 "ao_integrals.hpp"
+#include <Eigen/Eigen>
+
+namespace integrals::ao_integrals {
+
+using pt = simde::ERI2;
+
+namespace {
+
+auto desc = R"(
+Inverse Coulomb Metric
+---------------------
+)";
+
+}
+
+MODULE_CTOR(CoulombMetric) {
+    description(desc);
+    satisfies_property_type<pt>();
+    add_submodule<pt>("Two-center ERI");
+}
+
+MODULE_RUN(CoulombMetric) {
+    const auto& [braket] = pt::unwrap_inputs(inputs);
+    auto& eri2_mod       = submods.at("Two-center ERI");
+
+    // Get ERI2
+    const auto& I = eri2_mod.run_as<pt>(braket);
+
+    // Cholesky Decomp
+    tensorwrapper::allocator::Eigen<double> allocator(get_runtime());
+    const auto& eigen_I = allocator.rebind(I.buffer());
+    const auto* pI      = eigen_I.data();
+    const auto& shape_I = eigen_I.layout().shape().as_smooth();
+    auto rows           = shape_I.extent(0);
+    auto cols           = shape_I.extent(1);
+
+    using emat_t = Eigen::MatrixXd;
+    Eigen::Map<const emat_t> I_map(pI, rows, cols);
+    Eigen::LLT<emat_t> lltOfI(I_map);
+    emat_t U    = lltOfI.matrixU();
+    emat_t Linv = U.inverse().transpose();
+
+    // Wrap result
+    tensorwrapper::shape::Smooth matrix_shape{rows, cols};
+    tensorwrapper::layout::Physical matrix_layout(matrix_shape);
+    auto pM_buffer = allocator.allocate(matrix_layout);
+
+    for(auto i = 0; i < rows; ++i) {
+        for(auto j = 0; j < cols; ++j) { pM_buffer->at(i, j) = Linv(i, j); }
+    }
+    simde::type::tensor M(matrix_shape, std::move(pM_buffer));
+
+    auto rv = results();
+    return pt::wrap_results(rv, M);
+}
+
+} // namespace integrals::ao_integrals

src/integrals/ao_integrals/df_integral.cpp

@@ -0,0 +1,59 @@
+/*
+ * 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 "ao_integrals.hpp"
+
+namespace integrals::ao_integrals {
+
+using pt    = simde::ERI3;
+using pt_2c = simde::ERI2;
+
+namespace {
+
+auto desc = R"(
+Three-index ERI with Coulomb metric transformation
+---------------------
+)";
+
+}
+
+MODULE_CTOR(DFIntegral) {
+    description(desc);
+    satisfies_property_type<pt>();
+    add_submodule<pt>("Three-center ERI");
+    add_submodule<pt_2c>("Coulomb Metric");
+}
+
+MODULE_RUN(DFIntegral) {
+    const auto& [braket] = pt::unwrap_inputs(inputs);
+    auto bra             = braket.bra();
+    auto ket             = braket.ket();
+    auto& op             = braket.op();
+    auto& eri2_mod       = submods.at("Coulomb Metric");
+    auto& eri3_mod       = submods.at("Three-center ERI");
+
+    chemist::braket::BraKet aux_v_aux(bra, op, bra);
+    const auto& M = eri2_mod.run_as<pt_2c>(aux_v_aux);
+    const auto& I = eri3_mod.run_as<pt>(braket);
+
+    simde::type::tensor L;
+    L("i,k,l") = M("i,j") * I("j,k,l");
+
+    auto rv = results();
+    return pt::wrap_results(rv, std::move(L));
+}
+
+} // namespace integrals::ao_integrals

src/integrals/ao_integrals/j_density_fitted.cpp

@@ -0,0 +1,68 @@
+/*
+ * Copyright 2024 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 "ao_integrals.hpp"
+
+namespace integrals::ao_integrals {
+
+using pt    = simde::aos_j_e_aos;
+using pt_3c = simde::ERI3;
+using aos_t = simde::type::aos;
+
+namespace {
+
+auto desc = R"(
+Density Fitted J Builder
+---------------------
+)";
+
+}
+MODULE_CTOR(JDensityFitted) {
+    description(desc);
+    satisfies_property_type<pt>();
+    add_submodule<pt_3c>("DF ERI");
+    add_input<aos_t>("Auxiliary Basis Set");
+}
+
+MODULE_RUN(JDensityFitted) {
+    const auto&& [braket] = pt::unwrap_inputs(inputs);
+
+    auto bra          = braket.bra();
+    auto ket          = braket.ket();
+    const auto& j_hat = braket.op();
+    const auto& rho   = j_hat.rhs_particle();
+    const auto& p     = rho.value();
+    auto rho_aos      = rho.basis_set();
+    auto aux          = inputs.at("Auxiliary Basis Set").value<aos_t>();
+    auto& eri_mod     = submods.at("DF ERI");
+
+    simde::type::v_ee_type v_ee;
+    simde::type::aos_squared ij_pair(bra, ket);
+    simde::type::aos_squared kl_pair(rho_aos, rho_aos);
+    chemist::braket::BraKet aux_v_ij(aux, v_ee, ij_pair);
+    chemist::braket::BraKet aux_v_kl(aux, v_ee, kl_pair);
+    const auto& I_akl = eri_mod.run_as<pt_3c>(std::move(aux_v_kl));
+    const auto& I_aij = eri_mod.run_as<pt_3c>(std::move(aux_v_ij));
+
+    simde::type::tensor j;
+    j("a")   = p("k,l") * I_akl("a,k,l");
+    j("i,j") = j("a") * I_aij("a,i,j");
+
+    auto rv = results();
+    return pt::wrap_results(rv, std::move(j));
+}
+
+} // namespace integrals::ao_integrals

src/integrals/ao_integrals/k_density_fitted.cpp

@@ -0,0 +1,69 @@
+/*
+ * Copyright 2024 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 "ao_integrals.hpp"
+
+namespace integrals::ao_integrals {
+
+using pt    = simde::aos_k_e_aos;
+using pt_3c = simde::ERI3;
+using aos_t = simde::type::aos;
+
+namespace {
+
+auto desc = R"(
+Density Fitted K Builder
+---------------------
+)";
+
+}
+
+MODULE_CTOR(KDensityFitted) {
+    description(desc);
+    satisfies_property_type<pt>();
+    add_submodule<pt_3c>("DF ERI");
+    add_input<aos_t>("Auxiliary Basis Set");
+}
+
+MODULE_RUN(KDensityFitted) {
+    const auto&& [braket] = pt::unwrap_inputs(inputs);
+
+    auto bra          = braket.bra();
+    auto ket          = braket.ket();
+    const auto& j_hat = braket.op();
+    const auto& rho   = j_hat.rhs_particle();
+    const auto& p     = rho.value();
+    auto rho_aos      = rho.basis_set();
+    auto aux          = inputs.at("Auxiliary Basis Set").value<aos_t>();
+    auto& eri_mod     = submods.at("DF ERI");
+
+    simde::type::v_ee_type v_ee;
+    simde::type::aos_squared ik_pair(bra, rho_aos);
+    simde::type::aos_squared jl_pair(ket, rho_aos);
+    chemist::braket::BraKet aux_v_ik(aux, v_ee, ik_pair);
+    chemist::braket::BraKet aux_v_jl(aux, v_ee, jl_pair);
+    const auto& I_aik = eri_mod.run_as<pt_3c>(std::move(aux_v_ik));
+    const auto& I_ajl = eri_mod.run_as<pt_3c>(std::move(aux_v_jl));
+
+    simde::type::tensor k;
+    k("a,i,l") = p("k,l") * I_aik("a,i,k");
+    k("i,j")   = k("a,i,l") * I_ajl("a,j,l");
+
+    auto rv = results();
+    return pt::wrap_results(rv, std::move(k));
+}
+
+} // namespace integrals::ao_integrals

src/integrals/integrals_mm.cpp

@@ -27,11 +27,13 @@ namespace integrals {
  *  @throw none No throw guarantee
  */
 void set_defaults(pluginplay::ModuleManager& mm) {
-    const auto ao_driver = "AO integral driver";
-    mm.change_submod(ao_driver, "Kinetic", "Kinetic");
-    mm.change_submod(ao_driver, "Electron-Nuclear attraction", "Nuclear");
+    mm.change_submod("AO integral driver", "Kinetic", "Kinetic");
+    mm.change_submod("AO integral driver", "Electron-Nuclear attraction",
+                     "Nuclear");
     mm.change_submod("Four center J builder", "Four-center ERI", "ERI4");
     mm.change_submod("Four center K builder", "Four-center ERI", "ERI4");
+    mm.change_submod("Density Fitting Integral", "Three-center ERI", "ERI3");
+    mm.change_submod("Coulomb Metric", "Two-center ERI", "ERI2");
 }
 
 void load_modules(pluginplay::ModuleManager& mm) {

tests/cxx/unit/integrals/ao_integrals/coulomb_metric.cpp

@@ -0,0 +1,47 @@
+/*
+ * Copyright 2022 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 "../testing.hpp"
+
+TEST_CASE("Coulomb Metric") {
+    using test_pt = simde::ERI2;
+
+    pluginplay::ModuleManager mm;
+    integrals::load_modules(mm);
+    REQUIRE(mm.count("Coulomb Metric"));
+
+    // Get basis set
+    auto mol  = test::h2_molecule();
+    auto aobs = test::h2_sto3g_basis_set();
+
+    // Make AOS object
+    simde::type::aos aos(aobs);
+
+    // Make Operator
+    simde::type::v_ee_type op{};
+
+    // Make BraKet Input
+    chemist::braket::BraKet braket(aos, op, aos);
+
+    // Call module
+    auto T = mm.at("Coulomb Metric").run_as<test_pt>(braket);
+
+    auto t = test::eigen_tensor<2>(T.buffer());
+    REQUIRE(t(0, 0) == Catch::Approx(0.15824726).margin(1E-6));
+    REQUIRE(t(0, 1) == Catch::Approx(0.0).margin(1E-6));
+    REQUIRE(t(1, 0) == Catch::Approx(-0.23097095).margin(1E-6));
+    REQUIRE(t(1, 1) == Catch::Approx(0.27998174).margin(1E-6));
+}