Re-enables DFT

CMakeLists.txt

@@ -42,7 +42,7 @@ if("${BUILD_TAMM_SCF}")
     set(DEPENDENCIES simde gauxc tamm exachem chemcache)
     include(get_libint2)
 else()
-    set(DEPENDENCIES simde)
+    set(DEPENDENCIES simde gauxc)
 endif()
 
 foreach(dependency_i ${DEPENDENCIES})

src/scf/fock_operator/fock_operator.hpp

@@ -22,28 +22,45 @@ namespace scf::fock_operator {
 template<typename DensityType, typename ElectronType>
 DECLARE_MODULE(Restricted);
 
+template<typename DensityType, typename ElectronType>
+DECLARE_MODULE(RKohnSham);
+
 inline void load_modules(pluginplay::ModuleManager& mm) {
     using simde::type::decomposable_e_density;
     using simde::type::e_density;
     using simde::type::electron;
     using simde::type::many_electrons;
-    using e_many = Restricted<e_density, many_electrons>;
-    using e_e    = Restricted<e_density, electron>;
-    using d_many = Restricted<decomposable_e_density, many_electrons>;
-    using d_e    = Restricted<decomposable_e_density, electron>;
-
+    using e_many    = Restricted<e_density, many_electrons>;
+    using e_e       = Restricted<e_density, electron>;
+    using d_many    = Restricted<decomposable_e_density, many_electrons>;
+    using d_e       = Restricted<decomposable_e_density, electron>;
+    using e_ks_many = RKohnSham<e_density, many_electrons>;
+    using e_ks_e    = RKohnSham<e_density, electron>;
+    using d_ks_many = RKohnSham<decomposable_e_density, many_electrons>;
+    using d_ks_e    = RKohnSham<decomposable_e_density, electron>;
     mm.add_module<e_many>("AO Restricted Fock Op");
     mm.add_module<d_many>("Restricted Fock Op");
     mm.add_module<e_e>("AO Restricted One-Electron Fock Op");
     mm.add_module<d_e>("Restricted One-Electron Fock Op");
+    mm.add_module<e_ks_many>("AO Restricted Kohn-Sham Op");
+    mm.add_module<d_ks_many>("Restricted Kohn-Sham Op");
+    mm.add_module<e_ks_e>("AO Restricted One-Electron Kohn-Sham Op");
+    mm.add_module<d_ks_e>("Restricted One-Electron Kohn-Sham Op");
 }
 
-#define EXTERN_RESTRICTED(density)                                          \
+#define EXTERN_RESTRICTED(density)                                           \
     extern template struct Restricted<density, simde::type::many_electrons>; \
     extern template struct Restricted<density, simde::type::electron>
 
+#define EXTERN_RKOHN_SHAM(density)                                          \
+    extern template struct RKohnSham<density, simde::type::many_electrons>; \
+    extern template struct RKohnSham<density, simde::type::electron>
+
 EXTERN_RESTRICTED(simde::type::e_density);
 EXTERN_RESTRICTED(simde::type::decomposable_e_density);
+EXTERN_RKOHN_SHAM(simde::type::e_density);
+EXTERN_RKOHN_SHAM(simde::type::decomposable_e_density);
 
+#undef EXTERN_RKOHN_SHAM
 #undef EXTERN_RESTRICTED
 } // namespace scf::fock_operator

src/scf/fock_operator/fock_operator_common.hpp

@@ -0,0 +1,56 @@
+/*
+ * 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 <simde/simde.hpp>
+
+namespace scf::fock_operator {
+
+using namespace chemist::qm_operator;
+
+template<typename ElectronType>
+class FockOperatorCommon : public chemist::qm_operator::OperatorVisitor {
+public:
+    using T_e_term  = simde::type::T_e_type;
+    using V_en_term = simde::type::V_en_type;
+    using V_nn_term = simde::type::V_nn_type;
+
+    explicit FockOperatorCommon(simde::type::fock& F) : m_pF_(&F) {}
+
+    void run(const T_e_term& T_e) {
+        auto t = std::make_unique<Kinetic<ElectronType>>(get_e_(T_e));
+        m_pF_->emplace_back(1.0, std::move(t));
+    }
+
+    void run(const V_en_term& V_en) {
+        auto rhs     = V_en.rhs_particle().as_nuclei();
+        using v_type = Coulomb<ElectronType, simde::type::nuclei>;
+        auto v       = std::make_unique<v_type>(get_e_(V_en), rhs);
+        m_pF_->emplace_back(1.0, std::move(v));
+    }
+
+protected:
+    template<typename T>
+    auto get_e_(T&& op) {
+        if constexpr(std::is_same_v<ElectronType, simde::type::electron>) {
+            return simde::type::electron{};
+        } else {
+            return op.template at<0>();
+        }
+    }
+
+    simde::type::fock* m_pF_;
+};
+} // namespace scf::fock_operator

src/scf/fock_operator/restricted.cpp

@@ -15,56 +15,34 @@
  */
 
 #include "fock_operator.hpp"
-
+#include "fock_operator_common.hpp"
 namespace scf::fock_operator {
 
-using namespace chemist::qm_operator;
-
 template<typename DensityType, typename ElectronType>
-class RestrictedVisitor : public chemist::qm_operator::OperatorVisitor {
+class RestrictedVisitor : public FockOperatorCommon<ElectronType> {
+private:
+    using base_type = FockOperatorCommon<ElectronType>;
+
 public:
-    using T_e_term  = simde::type::T_e_type;
-    using V_en_term = simde::type::V_en_type;
     using V_ee_term = simde::type::V_ee_type;
-    using V_nn_term = simde::type::V_nn_type;
 
     RestrictedVisitor(simde::type::fock& F, const DensityType& rho) :
-      m_pF_(&F), m_prho_(&rho) {}
-
-    void run(const T_e_term& T_e) {
-        auto t = std::make_unique<Kinetic<ElectronType>>(get_e_(T_e));
-        m_pF_->emplace_back(1.0, std::move(t));
-    }
-
-    void run(const V_en_term& V_en) {
-        auto rhs     = V_en.rhs_particle().as_nuclei();
-        using v_type = Coulomb<ElectronType, simde::type::nuclei>;
-        auto v       = std::make_unique<v_type>(get_e_(V_en), rhs);
-        m_pF_->emplace_back(1.0, std::move(v));
-    }
+      base_type(F), m_prho_(&rho) {}
 
-    void run(const V_ee_term& V_ee) {
+    using base_type::run;
+    void run(const V_ee_term& V_ee) override {
         if(*m_prho_ == DensityType{}) return;
         using j_type = Coulomb<ElectronType, DensityType>;
         using k_type = Exchange<ElectronType, DensityType>;
 
-        auto j = std::make_unique<j_type>(get_e_(V_ee), *m_prho_);
-        auto k = std::make_unique<k_type>(get_e_(V_ee), *m_prho_);
-        m_pF_->emplace_back(2.0, std::move(j));
-        m_pF_->emplace_back(-1.0, std::move(k));
+        const auto& electrons = this->get_e_(V_ee);
+        auto j                = std::make_unique<j_type>(electrons, *m_prho_);
+        auto k                = std::make_unique<k_type>(electrons, *m_prho_);
+        this->m_pF_->emplace_back(2.0, std::move(j));
+        this->m_pF_->emplace_back(-1.0, std::move(k));
     }
 
 private:
-    template<typename T>
-    auto get_e_(T&& op) {
-        if constexpr(std::is_same_v<ElectronType, simde::type::electron>) {
-            return simde::type::electron{};
-        } else {
-            return op.template at<0>();
-        }
-    }
-
-    simde::type::fock* m_pF_;
     const DensityType* m_prho_;
 };
 

src/scf/fock_operator/rkohn_sham.cpp

@@ -0,0 +1,101 @@
+/*
+ * 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 "fock_operator.hpp"
+#include "fock_operator_common.hpp"
+
+namespace scf::fock_operator {
+using xc_functional = chemist::qm_operator::xc_functional;
+
+template<typename DensityType, typename ElectronType>
+class RKohnShamVisitor : public FockOperatorCommon<ElectronType> {
+private:
+    using base_type = FockOperatorCommon<ElectronType>;
+
+public:
+    using V_ee_term = simde::type::V_ee_type;
+
+    RKohnShamVisitor(xc_functional func, simde::type::fock& F,
+                     const DensityType& rho) :
+      base_type(F), m_func_(func), m_prho_(&rho) {}
+
+    using base_type::run;
+    void run(const V_ee_term& V_ee) override {
+        if(*m_prho_ == DensityType{}) return;
+        using j_type  = Coulomb<ElectronType, DensityType>;
+        using xc_type = ExchangeCorrelation<ElectronType, DensityType>;
+
+        const auto& electrons = this->get_e_(V_ee);
+        auto j                = std::make_unique<j_type>(electrons, *m_prho_);
+        auto xc = std::make_unique<xc_type>(m_func_, electrons, *m_prho_);
+        this->m_pF_->emplace_back(2.0, std::move(j));
+        this->m_pF_->emplace_back(1.0, std::move(xc));
+    }
+
+private:
+    xc_functional m_func_;
+    const DensityType* m_prho_;
+};
+
+const auto desc = R"(
+Restricted Kohn-Sham Operator Builder
+-------------------------------------
+
+This module builds the Kohn-Sham operator given a restricted density. Under this 
+assumption, the Kohn-Sham operator is created by mapping each term in the 
+electronic Hamiltonian to itself with the exception of the electron-electron
+repulsion. The electron-electron repulsion is mapped to 2J + X where J is the
+Coulomb operator for the electrons interacting with a density and X is the 
+exchange-correlation potential for the same density.
+
+N.b. Empty densities are treated as zero densities and this module will skip
+adding 2J + X if provided an empty density.
+)";
+
+template<typename DensityType, typename ElectronType>
+TEMPLATED_MODULE_CTOR(RKohnSham, DensityType, ElectronType) {
+    using pt = simde::FockOperator<DensityType>;
+    satisfies_property_type<pt>();
+    description(desc);
+
+    add_input<xc_functional>("XC Potential");
+}
+
+template<typename DensityType, typename ElectronType>
+TEMPLATED_MODULE_RUN(RKohnSham, DensityType, ElectronType) {
+    using pt = simde::FockOperator<DensityType>;
+
+    const auto& [H, rho] = pt::unwrap_inputs(inputs);
+    auto H_elec          = H.electronic_hamiltonian();
+    auto func            = inputs.at("XC Potential").value<xc_functional>();
+    simde::type::fock F;
+    RKohnShamVisitor<DensityType, ElectronType> visitor(func, F, rho);
+    H_elec.visit(visitor);
+
+    auto rv = results();
+    return pt::wrap_results(rv, F);
+}
+
+#define RESTRICTED(density)                                          \
+    template struct RKohnSham<density, simde::type::many_electrons>; \
+    template struct RKohnSham<density, simde::type::electron>
+
+RESTRICTED(simde::type::e_density);
+RESTRICTED(simde::type::decomposable_e_density);
+
+#undef RESTRICTED
+
+} // namespace scf::fock_operator

src/scf/matrix_builder/electronic_energy.cpp

@@ -25,19 +25,51 @@ const auto desc = R"(
 
 }
 
+using pluginplay::type::submodule_map;
 using simde::type::electronic_hamiltonian;
 
+using XC_e_t = simde::type::XC_e_type;
+
 template<typename WFType>
 using pt = simde::eval_braket<WFType, electronic_hamiltonian, WFType>;
 
+template<typename WFType>
+using xc_pt = simde::eval_braket<WFType, XC_e_t, WFType>;
+
 template<typename WFType>
 using det_pt = simde::eval_braket<WFType, simde::type::op_base_type, WFType>;
 
+template<typename WFType>
+class EnergyDispatcher : public chemist::qm_operator::OperatorVisitor {
+public:
+    EnergyDispatcher(const WFType* pbra, const WFType* pket,
+                     submodule_map* psubmods) :
+      chemist::qm_operator::OperatorVisitor(false),
+      m_pbra_(pbra),
+      m_pket_(pket),
+      m_psubmods_(psubmods) {}
+
+    void run(const XC_e_t& XC_e) {
+        chemist::braket::BraKet XC_ij(*m_pbra_, XC_e, *m_pket_);
+        m_o   = m_psubmods_->at("XC Energy").run_as<xc_pt<WFType>>(XC_ij);
+        m_ran = true;
+    }
+
+    bool m_ran = false;
+    simde::type::tensor m_o;
+
+private:
+    const WFType* m_pbra_;
+    const WFType* m_pket_;
+    submodule_map* m_psubmods_;
+};
+
 MODULE_CTOR(ElectronicEnergy) {
     using wf_type = simde::type::rscf_wf;
     description(desc);
     satisfies_property_type<pt<wf_type>>();
     add_submodule<det_pt<wf_type>>("determinant driver");
+    add_submodule<xc_pt<wf_type>>("XC Energy");
 }
 
 MODULE_RUN(ElectronicEnergy) {
@@ -47,18 +79,24 @@ MODULE_RUN(ElectronicEnergy) {
     const auto& H       = H_ij.op();
     const auto& ket     = H_ij.ket();
 
-    auto& mod = submods.at("determinant driver");
     simde::type::tensor e;
+    auto& det_mod = submods.at("determinant driver");
+    EnergyDispatcher<wf_type> visitor(&bra, &ket, &submods);
 
     auto n_ops = H.size();
     for(decltype(n_ops) i = 0; i < n_ops; ++i) {
         const auto& ci  = H.coefficient(i);
         const auto& O_i = H.get_operator(i);
 
-        chemist::braket::BraKet O_ij(bra, O_i, ket);
-        const auto o = mod.run_as<det_pt<wf_type>>(O_ij);
+        O_i.visit(visitor);
         simde::type::tensor temp;
-        temp("") = o("") * ci;
+        if(visitor.m_ran) {
+            temp("") = visitor.m_o("") * ci;
+        } else {
+            chemist::braket::BraKet O_ij(bra, O_i, ket);
+            const auto& o = det_mod.run_as<det_pt<wf_type>>(O_ij);
+            temp("")      = o("") * ci;
+        }
         if(i == 0)
             e = temp;
         else { e("") = e("") + temp(""); }

src/scf/matrix_builder/matrix_builder.hpp

@@ -35,14 +35,18 @@ inline void load_modules(pluginplay::ModuleManager& mm) {
 
 inline void set_defaults(pluginplay::ModuleManager& mm) {
     const auto ao_driver = "SCF integral driver";
-    mm.change_submod(ao_driver, "Fock matrix", "Fock Matrix Builder");
+    // mm.change_submod(ao_driver, "Fock matrix", "Fock Matrix Builder");
     mm.change_submod(ao_driver, "Density matrix", "Density matrix builder");
+    mm.change_submod(ao_driver, "XC Potential", "GauXC XC Potential");
+
+    mm.change_submod("Fock Matrix Builder", "Two center evaluator", ao_driver);
 
     const auto det_driver = "Determinant driver";
     mm.change_submod(det_driver, "Two center evaluator", ao_driver);
     mm.change_submod(det_driver, "Fock matrix", "Fock matrix builder");
 
     mm.change_submod("Electronic energy", "determinant driver", det_driver);
+    mm.change_submod("Electronic energy", "XC Energy", "GauXC XC Energy");
 }
 
 } // namespace scf::matrix_builder