Core Guess

src/scf/fock_operator/fock_operator.hpp

@@ -19,17 +19,31 @@
 
 namespace scf::fock_operator {
 
-template<typename DensityType>
+template<typename DensityType, typename ElectronType>
 DECLARE_MODULE(Restricted);
 
 inline void load_modules(pluginplay::ModuleManager& mm) {
     using simde::type::decomposable_e_density;
     using simde::type::e_density;
-    mm.add_module<Restricted<e_density>>("AO Restricted Fock Op");
-    mm.add_module<Restricted<decomposable_e_density>>("Restricted Fock Op");
+    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>;
+
+    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");
 }
 
-extern template class Restricted<simde::type::e_density>;
-extern template class Restricted<simde::type::decomposable_e_density>;
+#define EXTERN_RESTRICTED(density)                                          \
+    extern template class Restricted<density, simde::type::many_electrons>; \
+    extern template class Restricted<density, simde::type::electron>
+
+EXTERN_RESTRICTED(simde::type::e_density);
+EXTERN_RESTRICTED(simde::type::decomposable_e_density);
 
+#undef EXTERN_RESTRICTED
 } // namespace scf::fock_operator

src/scf/fock_operator/restricted.cpp

@@ -18,13 +18,11 @@
 
 namespace scf::fock_operator {
 
-template<typename DensityType>
+using namespace chemist::qm_operator;
+
+template<typename DensityType, typename ElectronType>
 class RestrictedVisitor : public chemist::qm_operator::OperatorVisitor {
 public:
-    using many_electrons = simde::type::many_electrons;
-    using J_type = chemist::qm_operator::Coulomb<many_electrons, DensityType>;
-    using K_type = chemist::qm_operator::Exchange<many_electrons, DensityType>;
-
     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;
@@ -33,17 +31,39 @@ class RestrictedVisitor : public chemist::qm_operator::OperatorVisitor {
     RestrictedVisitor(simde::type::fock& F, const DensityType& rho) :
       m_pF_(&F), m_prho_(&rho) {}
 
-    void run(const T_e_term& T_e) { m_pF_->emplace_back(1.0, T_e.clone()); }
+    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) { m_pF_->emplace_back(1.0, V_en.clone()); }
+    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));
+    }
 
     void run(const V_ee_term& V_ee) {
-        auto es = V_ee.at<0>();
-        m_pF_->emplace_back(2.0, std::make_unique<J_type>(es, *m_prho_));
-        m_pF_->emplace_back(-1.0, std::make_unique<K_type>(es, *m_prho_));
+        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));
     }
 
 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_;
 };
@@ -58,32 +78,40 @@ electronic Hamiltonian to itself with the exception of the electron-electron
 repulsion. The electron-electron repulsion is mapped to 2J-K where J is the
 Coulomb operator for the electrons interacting with a density and K is the 
 exchange operator for the same density.
+
+N.b. Empty densities are treated as zero densities and this module will skip
+adding 2J-K if provided an empty density.
 )";
 
-template<typename DensityType>
-TEMPLATED_MODULE_CTOR(Restricted, DensityType) {
+template<typename DensityType, typename ElectronType>
+TEMPLATED_MODULE_CTOR(Restricted, DensityType, ElectronType) {
     using pt = simde::FockOperator<DensityType>;
     satisfies_property_type<pt>();
     description(desc);
 }
 
-template<typename DensityType>
-TEMPLATED_MODULE_RUN(Restricted, DensityType) {
+template<typename DensityType, typename ElectronType>
+TEMPLATED_MODULE_RUN(Restricted, DensityType, ElectronType) {
     using pt = simde::FockOperator<DensityType>;
-    using simde::type::many_electrons;
 
     const auto& [H, rho] = pt::unwrap_inputs(inputs);
     auto H_elec          = H.electronic_hamiltonian();
 
     simde::type::fock F;
-    RestrictedVisitor<DensityType> visitor(F, rho);
+    RestrictedVisitor<DensityType, ElectronType> visitor(F, rho);
     H_elec.visit(visitor);
 
     auto rv = results();
     return pt::wrap_results(rv, F);
 }
 
-template class Restricted<simde::type::e_density>;
-template class Restricted<simde::type::decomposable_e_density>;
+#define RESTRICTED(density)                                          \
+    template class Restricted<density, simde::type::many_electrons>; \
+    template class Restricted<density, simde::type::electron>
+
+RESTRICTED(simde::type::e_density);
+RESTRICTED(simde::type::decomposable_e_density);
+
+#undef RESTRICTED
 
 } // namespace scf::fock_operator

src/scf/guess/core.cpp

@@ -0,0 +1,98 @@
+/*
+ * 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 "guess.hpp"
+
+namespace scf::guess {
+namespace {
+const auto desc = R"(
+Core Guess
+----------
+
+TODO: Write me!!!
+)";
+}
+
+using rscf_wf    = simde::type::rscf_wf;
+using density_t  = simde::type::decomposable_e_density;
+using pt         = simde::InitialGuess<rscf_wf>;
+using fock_op_pt = simde::FockOperator<density_t>;
+using update_pt  = simde::UpdateGuess<rscf_wf>;
+
+using simde::type::tensor;
+
+// TODO: move to chemist?
+struct NElectronCounter : public chemist::qm_operator::OperatorVisitor {
+    NElectronCounter() : chemist::qm_operator::OperatorVisitor(false) {}
+
+    void run(const simde::type::T_e_type& T_e) { set_n(T_e.particle().size()); }
+
+    void run(const simde::type::V_en_type& V_en) {
+        set_n(V_en.lhs_particle().size());
+    }
+
+    void run(const simde::type::V_ee_type& V_ee) {
+        set_n(V_ee.lhs_particle().size());
+        set_n(V_ee.rhs_particle().size());
+    }
+
+    void set_n(unsigned int n) {
+        if(n_electrons == 0)
+            n_electrons = n;
+        else if(n_electrons != n) {
+            throw std::runtime_error("Deduced a different number of electrons");
+        }
+    }
+
+    unsigned int n_electrons = 0;
+};
+
+MODULE_CTOR(Core) {
+    description(desc);
+    satisfies_property_type<pt>();
+    add_submodule<fock_op_pt>("Build Fock operator");
+    add_submodule<update_pt>("Guess updater");
+}
+
+MODULE_RUN(Core) {
+    const auto&& [H, aos] = pt::unwrap_inputs(inputs);
+
+    // Step 1: Build Fock Operator with zero density
+    density_t rho;
+    auto& fock_op_mod = submods.at("Build Fock operator");
+    const auto& f     = fock_op_mod.run_as<fock_op_pt>(H, rho);
+
+    // Step 2: Get number of electrons and occupations
+    simde::type::cmos cmos(tensor{}, aos, tensor{});
+    NElectronCounter visitor;
+    H.visit(visitor);
+    auto n_electrons = visitor.n_electrons;
+    if(n_electrons % 2 != 0)
+        throw std::runtime_error("Assumed even number of electrons");
+
+    typename rscf_wf::orbital_index_set_type occs;
+    using value_type = typename rscf_wf::orbital_index_set_type::value_type;
+    for(value_type i = 0; i < n_electrons / 2; ++i) occs.insert(i);
+
+    rscf_wf zero_guess(occs, cmos);
+    auto& update_mod = submods.at("Guess updater");
+    const auto& Psi0 = update_mod.run_as<update_pt>(f, zero_guess);
+
+    auto rv = results();
+    return pt::wrap_results(rv, Psi0);
+}
+
+} // namespace scf::guess

src/scf/guess/guess.hpp

@@ -0,0 +1,35 @@
+/*
+ * 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.
+ */
+
+#pragma once
+#include <simde/simde.hpp>
+
+namespace scf::guess {
+
+DECLARE_MODULE(Core);
+
+inline void load_modules(pluginplay::ModuleManager& mm) {
+    mm.add_module<Core>("Core guess");
+}
+
+inline void set_defaults(pluginplay::ModuleManager& mm) {
+    mm.change_submod("Core guess", "Build Fock operator",
+                     "Restricted One-Electron Fock Op");
+    mm.change_submod("Core guess", "Guess updater",
+                     "Diagonalization Fock update");
+}
+
+} // namespace scf::guess