diff --git a/src/scf/fock_operator/fock_operator.hpp b/src/scf/fock_operator/fock_operator.hpp index 48841b7..5276e29 100644 --- a/src/scf/fock_operator/fock_operator.hpp +++ b/src/scf/fock_operator/fock_operator.hpp @@ -19,17 +19,31 @@ namespace scf::fock_operator { -template +template DECLARE_MODULE(Restricted); inline void load_modules(pluginplay::ModuleManager& mm) { using simde::type::decomposable_e_density; using simde::type::e_density; - mm.add_module>("AO Restricted Fock Op"); - mm.add_module>("Restricted Fock Op"); + using simde::type::electron; + using simde::type::many_electrons; + using e_many = Restricted; + using e_e = Restricted; + using d_many = Restricted; + using d_e = Restricted; + + mm.add_module("AO Restricted Fock Op"); + mm.add_module("Restricted Fock Op"); + mm.add_module("AO Restricted One-Electron Fock Op"); + mm.add_module("Restricted One-Electron Fock Op"); } -extern template class Restricted; -extern template class Restricted; +#define EXTERN_RESTRICTED(density) \ + extern template class Restricted; \ + extern template class Restricted + +EXTERN_RESTRICTED(simde::type::e_density); +EXTERN_RESTRICTED(simde::type::decomposable_e_density); +#undef EXTERN_RESTRICTED } // namespace scf::fock_operator \ No newline at end of file diff --git a/src/scf/fock_operator/restricted.cpp b/src/scf/fock_operator/restricted.cpp index 0bb90e2..a82b94d 100644 --- a/src/scf/fock_operator/restricted.cpp +++ b/src/scf/fock_operator/restricted.cpp @@ -18,13 +18,11 @@ namespace scf::fock_operator { -template +using namespace chemist::qm_operator; + +template class RestrictedVisitor : public chemist::qm_operator::OperatorVisitor { public: - using many_electrons = simde::type::many_electrons; - using J_type = chemist::qm_operator::Coulomb; - using K_type = chemist::qm_operator::Exchange; - 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>(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; + auto v = std::make_unique(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(es, *m_prho_)); - m_pF_->emplace_back(-1.0, std::make_unique(es, *m_prho_)); + if(*m_prho_ == DensityType{}) return; + using j_type = Coulomb; + using k_type = Exchange; + + auto j = std::make_unique(get_e_(V_ee), *m_prho_); + auto k = std::make_unique(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 + auto get_e_(T&& op) { + if constexpr(std::is_same_v) { + 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 -TEMPLATED_MODULE_CTOR(Restricted, DensityType) { +template +TEMPLATED_MODULE_CTOR(Restricted, DensityType, ElectronType) { using pt = simde::FockOperator; satisfies_property_type(); description(desc); } -template -TEMPLATED_MODULE_RUN(Restricted, DensityType) { +template +TEMPLATED_MODULE_RUN(Restricted, DensityType, ElectronType) { using pt = simde::FockOperator; - using simde::type::many_electrons; const auto& [H, rho] = pt::unwrap_inputs(inputs); auto H_elec = H.electronic_hamiltonian(); simde::type::fock F; - RestrictedVisitor visitor(F, rho); + RestrictedVisitor visitor(F, rho); H_elec.visit(visitor); auto rv = results(); return pt::wrap_results(rv, F); } -template class Restricted; -template class Restricted; +#define RESTRICTED(density) \ + template class Restricted; \ + template class Restricted + +RESTRICTED(simde::type::e_density); +RESTRICTED(simde::type::decomposable_e_density); + +#undef RESTRICTED } // namespace scf::fock_operator \ No newline at end of file diff --git a/src/scf/guess/core.cpp b/src/scf/guess/core.cpp new file mode 100644 index 0000000..ace3108 --- /dev/null +++ b/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; +using fock_op_pt = simde::FockOperator; +using update_pt = simde::UpdateGuess; + +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(); + add_submodule("Build Fock operator"); + add_submodule("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(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(f, zero_guess); + + auto rv = results(); + return pt::wrap_results(rv, Psi0); +} + +} // namespace scf::guess \ No newline at end of file diff --git a/src/scf/guess/guess.hpp b/src/scf/guess/guess.hpp new file mode 100644 index 0000000..dec266e --- /dev/null +++ b/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 + +namespace scf::guess { + +DECLARE_MODULE(Core); + +inline void load_modules(pluginplay::ModuleManager& mm) { + mm.add_module("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 \ No newline at end of file diff --git a/src/scf/matrix_builder/ao_integrals_driver.cpp b/src/scf/matrix_builder/ao_integrals_driver.cpp new file mode 100644 index 0000000..b589980 --- /dev/null +++ b/src/scf/matrix_builder/ao_integrals_driver.cpp @@ -0,0 +1,115 @@ +/* + * 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 "matrix_builder.hpp" + +namespace scf::matrix_builder { + +namespace { + +const auto desc = R"( +AO Integrals Driver +------------------- +)"; + +} + +using pluginplay::type::submodule_map; +using simde::type::aos; +using simde::type::tensor; + +using pt = simde::aos_op_base_aos; +using t_e_pt = simde::aos_t_e_aos; +using v_en_pt = simde::aos_v_en_aos; +using j_e_pt = simde::aos_j_e_aos; +using k_e_pt = simde::aos_k_e_aos; +using f_e_pt = simde::aos_f_e_aos; + +class AODispatcher : public chemist::qm_operator::OperatorVisitor { +public: + using t_e_type = simde::type::t_e_type; + using v_en_type = simde::type::v_en_type; + using j_e_type = simde::type::j_e_type; + using k_e_type = simde::type::k_e_type; + using f_e_type = simde::type::fock; + + using submods_type = pluginplay::type::submodule_map; + + AODispatcher(const aos& bra, const aos& ket, submodule_map& submods, + tensor& t) : + m_pbra_(&bra), m_pket_(&ket), m_psubmods_(&submods), m_ptensor_(&t) {} + + void run(const t_e_type& t_e) { + chemist::braket::BraKet input(*m_pbra_, t_e, *m_pket_); + *m_ptensor_ = m_psubmods_->at("Kinetic").run_as(input); + } + + void run(const v_en_type& v_en) { + chemist::braket::BraKet input(*m_pbra_, v_en, *m_pket_); + const auto key = "Electron-Nuclear attraction"; + *m_ptensor_ = m_psubmods_->at(key).run_as(input); + } + + void run(const j_e_type& j_e) { + chemist::braket::BraKet input(*m_pbra_, j_e, *m_pket_); + const auto key = "Coulomb matrix"; + *m_ptensor_ = m_psubmods_->at(key).run_as(input); + } + + void run(const k_e_type& k_e) { + chemist::braket::BraKet input(*m_pbra_, k_e, *m_pket_); + const auto key = "Exchange matrix"; + *m_ptensor_ = m_psubmods_->at(key).run_as(input); + } + + // void run(const f_e_type& f_e) { + // chemist::braket::BraKet input(*m_pbra_, f_e, *m_pket_); + // const auto key = "Fock matrix"; + // *m_ptensor_ = m_psubmods_->at(key).run_as(input); + // } + +private: + const aos* m_pbra_; + const aos* m_pket_; + submodule_map* m_psubmods_; + tensor* m_ptensor_; +}; + +MODULE_CTOR(AOIntegralsDriver) { + description(desc); + satisfies_property_type(); + add_submodule("Kinetic"); + add_submodule("Electron-Nuclear attraction"); + add_submodule("Coulomb matrix"); + add_submodule("Exchange matrix"); + // add_submodule("Fock matrix"); +} + +MODULE_RUN(AOIntegralsDriver) { + const auto&& [braket] = pt::unwrap_inputs(inputs); + const auto& bra = braket.bra(); + const auto& op = braket.op(); + const auto& ket = braket.ket(); + + tensor t; + AODispatcher visitor(bra, ket, submods, t); + op.visit(visitor); + + auto rv = results(); + return pt::wrap_results(rv, std::move(t)); +} + +} // namespace scf::matrix_builder \ No newline at end of file diff --git a/src/scf/matrix_builder/fock.cpp b/src/scf/matrix_builder/fock.cpp new file mode 100644 index 0000000..fc28897 --- /dev/null +++ b/src/scf/matrix_builder/fock.cpp @@ -0,0 +1,71 @@ +/* + * 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 "matrix_builder.hpp" + +namespace scf::matrix_builder { +namespace { + +const auto desc = R"( +)"; + +} + +using pt = simde::aos_f_e_aos; +using ao_pt = simde::aos_op_base_aos; + +MODULE_CTOR(Fock) { + description(desc); + satisfies_property_type(); + add_submodule("Two center evaluator"); +} + +MODULE_RUN(Fock) { + const auto&& [braket] = pt::unwrap_inputs(inputs); + const auto& bra_aos = braket.bra(); + const auto& f = braket.op(); + const auto& ket_aos = braket.ket(); + auto& ao_dispatcher = submods.at("Two center evaluator"); + + simde::type::tensor F; + auto n_terms = f.size(); + if(n_terms > 0) { + // Initialize F to zero-th term + const auto c0 = f.coefficient(0); + const auto& op_0 = f.get_operator(0); + chemist::braket::BraKet term0(bra_aos, op_0, ket_aos); + F = ao_dispatcher.run_as(term0); + + using allocator_type = tensorwrapper::allocator::Eigen; + auto& f_buffer = allocator_type::rebind(F.buffer()); + f_buffer.value() = f_buffer.value() * c0; + + for(decltype(n_terms) i = 1; i < n_terms; ++i) { + const auto ci = f.coefficient(i); + const auto& op_i = f.get_operator(i); + chemist::braket::BraKet termi(bra_aos, op_i, ket_aos); + auto matrix = ao_dispatcher.run_as(termi); + + auto& matrix_buffer = allocator_type::rebind(matrix.buffer()); + f_buffer.value() += ci * matrix_buffer.value(); + } + } + + auto rv = results(); + return pt::wrap_results(rv, F); +} + +} // namespace scf::matrix_builder \ No newline at end of file diff --git a/src/scf/matrix_builder/j_four_center.cpp b/src/scf/matrix_builder/j_four_center.cpp new file mode 100644 index 0000000..909e246 --- /dev/null +++ b/src/scf/matrix_builder/j_four_center.cpp @@ -0,0 +1,83 @@ +/* + * 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 "matrix_builder.hpp" + +namespace scf::matrix_builder { + +using pt = simde::aos_j_e_aos; +using pt_4c = simde::ERI4; + +namespace { + +auto desc = R"( +Four-Center J Builder +--------------------- +)"; + +} +MODULE_CTOR(JFourCenter) { + description(desc); + satisfies_property_type(); + add_submodule("Four-center ERI"); +} + +MODULE_RUN(JFourCenter) { + const auto&& [braket] = pt::unwrap_inputs(inputs); + // TODO: avoid copying AOs + simde::type::aos bra_e0 = braket.bra(); + const auto& j_hat = braket.op(); + simde::type::aos ket_e0 = braket.ket(); + const auto& rho = j_hat.rhs_particle(); + simde::type::aos aos_e1 = rho.basis_set(); + const auto& p = rho.value(); + auto& eri_mod = submods.at("Four-center ERI"); + + // auto aos2_v_aos2 = (bra_e0 * ket_e0 | v_ee | aos_e1 * aos_e1); + simde::type::v_ee_type v_ee; + simde::type::aos_squared e0_pair(bra_e0, ket_e0); + simde::type::aos_squared e1_pair(aos_e1, aos_e1); + chemist::braket::BraKet aos2_v_aos2(e0_pair, v_ee, e1_pair); + const auto& I = eri_mod.run_as(std::move(aos2_v_aos2)); + + // This goes away when j("m,n") = p("l,s")*I("m,n,s,l") works + // { + using eri_alloc = tensorwrapper::allocator::Eigen; + using rho_alloc = tensorwrapper::allocator::Eigen; + using index_pair = Eigen::IndexPair; + using array_t = Eigen::array; + + const auto& I_eigen = eri_alloc::rebind(I.buffer()).value(); + const auto& p_buffer = rho_alloc::rebind(p.buffer()); + const auto& p_eigen = p_buffer.value(); + + array_t contract_modes{index_pair{0, 3}, index_pair(1, 2)}; + using eigen_buffer_type = typename rho_alloc::eigen_buffer_type; + using eigen_tensor_type = typename eigen_buffer_type::data_type; + + eigen_tensor_type j_eigen = p_eigen.contract(I_eigen, contract_modes); + auto j_buffer = std::make_unique(std::move(j_eigen), + p_buffer.layout()); + + using logical_type = tensorwrapper::layout::Logical; + auto playout = p.logical_layout().clone_as(); + simde::type::tensor j(std::move(playout), std::move(j_buffer)); + //} + auto rv = results(); + return pt::wrap_results(rv, std::move(j)); +} + +} // namespace scf::matrix_builder \ No newline at end of file diff --git a/src/scf/matrix_builder/k_four_center.cpp b/src/scf/matrix_builder/k_four_center.cpp new file mode 100644 index 0000000..19f9b9f --- /dev/null +++ b/src/scf/matrix_builder/k_four_center.cpp @@ -0,0 +1,84 @@ +/* + * 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 "matrix_builder.hpp" + +namespace scf::matrix_builder { + +using pt = simde::aos_k_e_aos; +using pt_4c = simde::ERI4; + +namespace { + +auto desc = R"( +Four-Center K Builder +--------------------- +)"; + +} + +MODULE_CTOR(KFourCenter) { + description(desc); + satisfies_property_type(); + add_submodule("Four-center ERI"); +} + +MODULE_RUN(KFourCenter) { + const auto&& [braket] = pt::unwrap_inputs(inputs); + // TODO: avoid copying AOs + simde::type::aos bra_e0 = braket.bra(); + const auto& k_hat = braket.op(); + simde::type::aos ket_e0 = braket.ket(); + const auto& rho = k_hat.rhs_particle(); + simde::type::aos aos_e1 = rho.basis_set(); + const auto& p = rho.value(); + auto& eri_mod = submods.at("Four-center ERI"); + + // auto aos2_v_aos2 = (bra_e0 * aos_e1 | v_ee | aos_e1 * ket_e0); + simde::type::v_ee_type v_ee; + simde::type::aos_squared e0_pair(bra_e0, aos_e1); + simde::type::aos_squared e1_pair(aos_e1, ket_e0); + chemist::braket::BraKet aos2_v_aos2(e0_pair, v_ee, e1_pair); + const auto& I = eri_mod.run_as(std::move(aos2_v_aos2)); + + // This goes away when k("m,n") = p("l,s")*I("m,l,s,n") works + // { + using eri_alloc = tensorwrapper::allocator::Eigen; + using rho_alloc = tensorwrapper::allocator::Eigen; + using index_pair = Eigen::IndexPair; + using array_t = Eigen::array; + + const auto& I_eigen = eri_alloc::rebind(I.buffer()).value(); + const auto& p_buffer = rho_alloc::rebind(p.buffer()); + const auto& p_eigen = p_buffer.value(); + + array_t contract_modes{index_pair{0, 1}, index_pair(1, 2)}; + using eigen_buffer_type = typename rho_alloc::eigen_buffer_type; + using eigen_tensor_type = typename eigen_buffer_type::data_type; + + eigen_tensor_type k_eigen = p_eigen.contract(I_eigen, contract_modes); + auto k_buffer = std::make_unique(std::move(k_eigen), + p_buffer.layout()); + + using logical_type = tensorwrapper::layout::Logical; + auto playout = p.logical_layout().clone_as(); + simde::type::tensor k(std::move(playout), std::move(k_buffer)); + //} + auto rv = results(); + return pt::wrap_results(rv, std::move(k)); +} + +} // namespace scf::matrix_builder \ No newline at end of file diff --git a/src/scf/matrix_builder/matrix_builder.hpp b/src/scf/matrix_builder/matrix_builder.hpp new file mode 100644 index 0000000..b1e6345 --- /dev/null +++ b/src/scf/matrix_builder/matrix_builder.hpp @@ -0,0 +1,44 @@ +/* + * 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 + +namespace scf::matrix_builder { + +DECLARE_MODULE(AOIntegralsDriver); +DECLARE_MODULE(Fock); +DECLARE_MODULE(JFourCenter); +DECLARE_MODULE(KFourCenter); + +inline void load_modules(pluginplay::ModuleManager& mm) { + mm.add_module("AO integral driver"); + mm.add_module("Fock matrix builder"); + mm.add_module("Four center J builder"); + mm.add_module("Four center K builder"); +} + +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"); + // TODO: Re-enable when PluginPlay doesn't choke on loops in modules + // mm.change_submod(ao_driver, "Fock matrix", "Fock Matrix Builder"); + + mm.change_submod("Fock matrix builder", "Two center evaluator", ao_driver); +} + +} // namespace scf::matrix_builder \ No newline at end of file diff --git a/src/scf/scf_mm.cpp b/src/scf/scf_mm.cpp index 1cc4260..34019be 100644 --- a/src/scf/scf_mm.cpp +++ b/src/scf/scf_mm.cpp @@ -16,7 +16,10 @@ #include "eigen_solver/eigen_solver.hpp" #include "fock_operator/fock_operator.hpp" +#include "guess/guess.hpp" +#include "matrix_builder/matrix_builder.hpp" #include "scf_modules.hpp" +#include "update/update.hpp" #include namespace scf { @@ -24,10 +27,18 @@ namespace scf { void load_modules(pluginplay::ModuleManager& mm) { eigen_solver::load_modules(mm); fock_operator::load_modules(mm); + guess::load_modules(mm); + matrix_builder::load_modules(mm); + update::load_modules(mm); + mm.add_module("Coulomb's Law"); #ifdef BUILD_TAMM_SCF mm.add_module("SCF Energy via TAMM"); #endif + + guess::set_defaults(mm); + matrix_builder::set_defaults(mm); + update::set_defaults(mm); } } // namespace scf diff --git a/src/scf/update/diagonalization.cpp b/src/scf/update/diagonalization.cpp new file mode 100644 index 0000000..9beeaf3 --- /dev/null +++ b/src/scf/update/diagonalization.cpp @@ -0,0 +1,65 @@ +/* + * 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 "update.hpp" + +namespace scf::update { + +using rscf_wf = simde::type::rscf_wf; +using fock_matrix_pt = simde::aos_f_e_aos; +using pt = simde::UpdateGuess; +using diagonalizer_pt = simde::GeneralizedEigenSolve; +using s_pt = simde::aos_s_e_aos; + +const auto desc = R"( +)"; + +MODULE_CTOR(Diagonalization) { + description(desc); + satisfies_property_type(); + add_submodule("Fock matrix builder"); + add_submodule("Diagonalizer"); + add_submodule("Overlap matrix builder"); +} + +MODULE_RUN(Diagonalization) { + const auto&& [f, old_guess] = pt::unwrap_inputs(inputs); + const auto& aos = old_guess.orbitals().from_space(); + + // Comput F + chemist::braket::BraKet f_mn(aos, f, aos); + auto& fock_matrix_mod = submods.at("Fock matrix builder"); + const auto& f_matrix = fock_matrix_mod.run_as(f_mn); + + // Compute S + chemist::braket::BraKet s_mn(aos, simde::type::s_e_type{}, aos); + auto& s_mod = submods.at("Overlap matrix builder"); + const auto& s_matrix = s_mod.run_as(s_mn); + + // Diagonalize + auto& diagonalizer_mod = submods.at("Diagonalizer"); + const auto&& [evalues, evectors] = + diagonalizer_mod.run_as(f_matrix, s_matrix); + + // Create new guess + simde::type::cmos cmos(evalues, aos, evectors); + rscf_wf new_guess(old_guess.orbital_indices(), cmos); + + auto rv = results(); + return pt::wrap_results(rv, new_guess); +} + +} // namespace scf::update \ No newline at end of file diff --git a/src/scf/update/update.hpp b/src/scf/update/update.hpp new file mode 100644 index 0000000..20dceff --- /dev/null +++ b/src/scf/update/update.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 + +namespace scf::update { + +DECLARE_MODULE(Diagonalization); + +inline void load_modules(pluginplay::ModuleManager& mm) { + mm.add_module("Diagonalization Fock update"); +} + +inline void set_defaults(pluginplay::ModuleManager& mm) { + mm.change_submod("Diagonalization Fock update", "Diagonalizer", + "Generalized eigensolve via Eigen"); + mm.change_submod("Diagonalization Fock update", "Fock matrix builder", + "Fock matrix builder"); +} + +} // namespace scf::update \ No newline at end of file diff --git a/tests/cxx/integration_tests/guess/core.cpp b/tests/cxx/integration_tests/guess/core.cpp new file mode 100644 index 0000000..591c657 --- /dev/null +++ b/tests/cxx/integration_tests/guess/core.cpp @@ -0,0 +1,42 @@ +/* + * 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 "../integration_tests.hpp" + +using rscf_wf = simde::type::rscf_wf; +using pt = simde::InitialGuess; +using simde::type::tensor; + +TEST_CASE("Core") { + auto mm = test_scf::load_modules(); + auto aos = test_scf::h2_aos(); + auto H = test_scf::h2_hamiltonian(); + + auto mod = mm.at("Core guess"); + auto psi = mod.run_as(H, aos); + + typename rscf_wf::orbital_index_set_type occs{0}; + REQUIRE(psi.orbital_indices() == occs); + REQUIRE(psi.orbitals().from_space() == aos); + const auto& evals = psi.orbitals().diagonalized_matrix(); + using allocator_type = tensorwrapper::allocator::Eigen; + const auto& eval_buffer = allocator_type::rebind(evals.buffer()); + + const auto tol = 1E-6; + using Catch::Matchers::WithinAbs; + REQUIRE_THAT(eval_buffer.value()(0), WithinAbs(-1.25330893, tol)); + REQUIRE_THAT(eval_buffer.value()(1), WithinAbs(-0.47506974, tol)); +} \ No newline at end of file diff --git a/tests/cxx/integration_tests/integration_tests.hpp b/tests/cxx/integration_tests/integration_tests.hpp new file mode 100644 index 0000000..eb9be3e --- /dev/null +++ b/tests/cxx/integration_tests/integration_tests.hpp @@ -0,0 +1,44 @@ +/* + * 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 "../test_scf.hpp" +#include +#include +#include + +namespace test_scf { + +/// Factors out setting submodules for SCF plugin from other plugins +inline auto load_modules() { + pluginplay::ModuleManager mm; + scf::load_modules(mm); + integrals::load_modules(mm); + + mm.change_submod("Four center J builder", "Four-center ERI", "ERI4"); + mm.change_submod("Four center K builder", "Four-center ERI", "ERI4"); + + const auto ao_driver_key = "AO integral driver"; + mm.change_submod(ao_driver_key, "Kinetic", "Kinetic"); + mm.change_submod(ao_driver_key, "Electron-Nuclear attraction", "Nuclear"); + + mm.change_submod("Diagonalization Fock update", "Overlap matrix builder", + "Overlap"); + + return mm; +} + +} // namespace test_scf \ No newline at end of file diff --git a/tests/cxx/integration_tests/matrix_builder/ao_integrals_driver.cpp b/tests/cxx/integration_tests/matrix_builder/ao_integrals_driver.cpp new file mode 100644 index 0000000..e90765e --- /dev/null +++ b/tests/cxx/integration_tests/matrix_builder/ao_integrals_driver.cpp @@ -0,0 +1,104 @@ +/* + * 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 "../integration_tests.hpp" + +using pt = simde::aos_op_base_aos; +using simde::type::tensor; + +namespace { + +void compare_matrices(const tensor& A, const tensor& A_corr) { + using Catch::Matchers::WithinAbs; + using alloc_type = tensorwrapper::allocator::Eigen; + const auto& A_buffer = alloc_type::rebind(A.buffer()); + const auto& A_corr_buffer = alloc_type::rebind(A_corr.buffer()); + const auto& A_eigen = A_buffer.value(); + const auto& A_corr_eigen = A_corr_buffer.value(); + + const auto tol = 1E-6; + + REQUIRE_THAT(A_eigen(0, 0), WithinAbs(A_corr_eigen(0, 0), 1E-6)); + REQUIRE_THAT(A_eigen(0, 1), WithinAbs(A_corr_eigen(0, 1), 1E-6)); + REQUIRE_THAT(A_eigen(1, 0), WithinAbs(A_corr_eigen(1, 0), 1E-6)); + REQUIRE_THAT(A_eigen(1, 1), WithinAbs(A_corr_eigen(1, 1), 1E-6)); +} + +} // namespace + +using erased_type = + chemist::braket::BraKet; + +TEST_CASE("AOIntegralsDriver") { + auto mm = test_scf::load_modules(); + auto aos = test_scf::h2_aos(); + auto mod = mm.at("AO integral driver"); + simde::type::electron e; + auto rho = test_scf::h2_density(); + + SECTION("Calling Kinetic") { + auto& tmod = mm.at("Kinetic"); + simde::type::t_e_type t_e(e); + chemist::braket::BraKet braket(aos, t_e, aos); + erased_type copy_braket(braket); + const auto& T = mod.run_as(copy_braket); + const auto& T_corr = tmod.run_as(braket); + compare_matrices(T, T_corr); + } + + SECTION("Calling Electron-Nuclear Attraction") { + auto& tmod = mm.at("Nuclear"); + auto h2_nuclei = test_scf::make_h2(); + simde::type::v_en_type v_en(e, h2_nuclei); + chemist::braket::BraKet braket(aos, v_en, aos); + erased_type copy_braket(braket); + const auto& V = mod.run_as(copy_braket); + const auto& V_corr = tmod.run_as(braket); + compare_matrices(V, V_corr); + } + + SECTION("Calling J Matrix") { + auto& jmod = mm.at("Four center J builder"); + simde::type::j_e_type j_e(e, rho); + chemist::braket::BraKet braket(aos, j_e, aos); + erased_type copy_braket(braket); + const auto& J = mod.run_as(copy_braket); + const auto& J_corr = jmod.run_as(braket); + compare_matrices(J, J_corr); + } + + SECTION("Calling K Matrix") { + auto& kmod = mm.at("Four center K builder"); + simde::type::k_e_type k_e(e, rho); + chemist::braket::BraKet braket(aos, k_e, aos); + erased_type copy_braket(braket); + const auto& K = mod.run_as(copy_braket); + const auto& K_corr = kmod.run_as(braket); + compare_matrices(K, K_corr); + } + + // Re-enable when PluginPlay doesn't choke on loops in modules + // SECTION("Calling Fock Matrix") { + // auto& fmod = mm.at("Fock matrix builder"); + // auto f_e = test_scf::h2_fock(); + // chemist::braket::BraKet braket(aos, f_e, aos); + // erased_type copy_braket(braket); + // const auto& F = mod.run_as(copy_braket); + // const auto& F_corr = fmod.run_as(braket); + // compare_matrices(F, F_corr); + // } +} \ No newline at end of file diff --git a/tests/cxx/integration_tests/matrix_builder/fock.cpp b/tests/cxx/integration_tests/matrix_builder/fock.cpp new file mode 100644 index 0000000..8283428 --- /dev/null +++ b/tests/cxx/integration_tests/matrix_builder/fock.cpp @@ -0,0 +1,61 @@ +/* + * 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 "../integration_tests.hpp" + +using pt = simde::aos_f_e_aos; + +using simde::type::t_e_type; +using simde::type::v_en_type; + +TEST_CASE("Fock Matrix Builder") { + auto mm = test_scf::load_modules(); + auto& mod = mm.at("Fock Matrix Builder"); + auto aos = test_scf::h2_aos(); + + SECTION("No J or K") { + auto h2 = test_scf::make_h2(); + simde::type::electron e; + + simde::type::fock f_e; + f_e.emplace_back(1.0, std::make_unique(e)); + f_e.emplace_back(1.0, std::make_unique(e, h2)); + const auto& F = mod.run_as(chemist::braket::BraKet(aos, f_e, aos)); + + using alloc_type = tensorwrapper::allocator::Eigen; + const auto& F_eigen = alloc_type::rebind(F.buffer()); + using Catch::Matchers::WithinAbs; + REQUIRE_THAT(F_eigen.value()(0, 0), WithinAbs(-1.120958, 1E-6)); + REQUIRE_THAT(F_eigen.value()(0, 1), WithinAbs(-0.959374, 1E-6)); + REQUIRE_THAT(F_eigen.value()(1, 0), WithinAbs(-0.959374, 1E-6)); + REQUIRE_THAT(F_eigen.value()(1, 1), WithinAbs(-1.120958, 1E-6)); + } + + SECTION("With J and K") { + auto f_e = test_scf::h2_fock(); + + const auto& F = mod.run_as(chemist::braket::BraKet(aos, f_e, aos)); + + using alloc_type = tensorwrapper::allocator::Eigen; + const auto& F_eigen = alloc_type::rebind(F.buffer()); + + using Catch::Matchers::WithinAbs; + REQUIRE_THAT(F_eigen.value()(0, 0), WithinAbs(-0.319459, 1E-6)); + REQUIRE_THAT(F_eigen.value()(0, 1), WithinAbs(-0.571781, 1E-6)); + REQUIRE_THAT(F_eigen.value()(1, 0), WithinAbs(-0.571781, 1E-6)); + REQUIRE_THAT(F_eigen.value()(1, 1), WithinAbs(-0.319459, 1E-6)); + } +} \ No newline at end of file diff --git a/tests/cxx/integration_tests/matrix_builder/j_four_center.cpp b/tests/cxx/integration_tests/matrix_builder/j_four_center.cpp new file mode 100644 index 0000000..ea7edf3 --- /dev/null +++ b/tests/cxx/integration_tests/matrix_builder/j_four_center.cpp @@ -0,0 +1,36 @@ +/* + * 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 "../integration_tests.hpp" + +using pt = simde::aos_j_e_aos; + +TEST_CASE("JFourCenter") { + auto mm = test_scf::load_modules(); + auto& mod = mm.at("Four center J builder"); + auto aos = test_scf::h2_aos(); + + simde::type::j_e_type j_e(simde::type::electron{}, test_scf::h2_density()); + const auto& J = mod.run_as(chemist::braket::BraKet(aos, j_e, aos)); + + using alloc_type = tensorwrapper::allocator::Eigen; + const auto& J_eigen = alloc_type::rebind(J.buffer()); + using Catch::Matchers::WithinAbs; + REQUIRE_THAT(J_eigen.value()(0, 0), WithinAbs(0.71438149, 1E-6)); + REQUIRE_THAT(J_eigen.value()(0, 1), WithinAbs(0.47471072, 1E-6)); + REQUIRE_THAT(J_eigen.value()(1, 0), WithinAbs(0.47471072, 1E-6)); + REQUIRE_THAT(J_eigen.value()(1, 1), WithinAbs(0.71438149, 1E-6)); +} \ No newline at end of file diff --git a/tests/cxx/integration_tests/matrix_builder/k_four_center.cpp b/tests/cxx/integration_tests/matrix_builder/k_four_center.cpp new file mode 100644 index 0000000..73f3580 --- /dev/null +++ b/tests/cxx/integration_tests/matrix_builder/k_four_center.cpp @@ -0,0 +1,36 @@ +/* + * 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 "../integration_tests.hpp" + +using pt = simde::aos_k_e_aos; + +TEST_CASE("KFourCenter") { + auto mm = test_scf::load_modules(); + auto& mod = mm.at("Four center K builder"); + auto aos = test_scf::h2_aos(); + + simde::type::k_e_type k_e(simde::type::electron{}, test_scf::h2_density()); + const auto& K = mod.run_as(chemist::braket::BraKet(aos, k_e, aos)); + + using alloc_type = tensorwrapper::allocator::Eigen; + const auto& K_eigen = alloc_type::rebind(K.buffer()); + using Catch::Matchers::WithinAbs; + REQUIRE_THAT(K_eigen.value()(0, 0), WithinAbs(0.627264, 1E-6)); + REQUIRE_THAT(K_eigen.value()(0, 1), WithinAbs(0.561828, 1E-6)); + REQUIRE_THAT(K_eigen.value()(1, 0), WithinAbs(0.561828, 1E-6)); + REQUIRE_THAT(K_eigen.value()(1, 1), WithinAbs(0.627264, 1E-6)); +} \ No newline at end of file diff --git a/tests/cxx/integration_tests/test_scf.cpp b/tests/cxx/integration_tests/test_scf.cpp index 00c7b02..2d9494c 100644 --- a/tests/cxx/integration_tests/test_scf.cpp +++ b/tests/cxx/integration_tests/test_scf.cpp @@ -21,54 +21,58 @@ #include #include -TEST_CASE("SCF") { +#ifdef BUILD_TAMM_SCF - // Populate modules - pluginplay::ModuleManager mm; - chemcache::load_modules(mm); - scf::load_modules(mm); +TEST_CASE("SCF") { + // Populate modules + pluginplay::ModuleManager mm; + chemcache::load_modules(mm); + scf::load_modules(mm); - // Create ChemicalSystem - std::string mol_name = "water"; - auto mol = mm.at("NWX Molecules").run_as(mol_name); - simde::type::chemical_system cs(mol); + // Create ChemicalSystem + std::string mol_name = "water"; + auto mol = + mm.at("NWX Molecules").run_as(mol_name); + simde::type::chemical_system cs(mol); - // Create BasisSet - std::string basis_name = + // Create BasisSet + std::string basis_name = "sto-3g"; // This is the only supported basis in ChemCache - auto aos = mm.at(basis_name).run_as(mol); + auto aos = mm.at(basis_name).run_as(mol); - // Run module - mm.change_input("SCF Energy via TAMM", "molecule_name", mol_name); - auto E = mm.at("SCF Energy via TAMM").run_as(aos, cs); - std::cout << "SCF Energy = " << E << " Hartree" << std::endl; + // Run module + mm.change_input("SCF Energy via TAMM", "molecule_name", mol_name); + auto E = mm.at("SCF Energy via TAMM").run_as(aos, cs); + std::cout << "SCF Energy = " << E << " Hartree" << std::endl; - REQUIRE(E == Catch::Approx(-74.3670617803483).margin(1.0e-6)); + REQUIRE(E == Catch::Approx(-74.3670617803483).margin(1.0e-6)); } TEST_CASE("DFT") { + // Populate modules + pluginplay::ModuleManager mm; + chemcache::load_modules(mm); + scf::load_modules(mm); - // Populate modules - pluginplay::ModuleManager mm; - chemcache::load_modules(mm); - scf::load_modules(mm); - - // Create ChemicalSystem - std::string mol_name = "water"; - auto mol = mm.at("NWX Molecules").run_as(mol_name); - simde::type::chemical_system cs(mol); + // Create ChemicalSystem + std::string mol_name = "water"; + auto mol = + mm.at("NWX Molecules").run_as(mol_name); + simde::type::chemical_system cs(mol); - // Create BasisSet - std::string basis_name = + // Create BasisSet + std::string basis_name = "sto-3g"; // This is the only supported basis in ChemCache - auto aos = mm.at(basis_name).run_as(mol); + auto aos = mm.at(basis_name).run_as(mol); - // Run module - std::vector xc_type = {"pbe0"}; - mm.change_input("SCF Energy via TAMM", "xc_type", xc_type); - mm.change_input("SCF Energy via TAMM", "molecule_name", mol_name); - auto E = mm.at("SCF Energy via TAMM").run_as(aos, cs); - std::cout << "SCF Energy = " << E << " Hartree" << std::endl; + // Run module + std::vector xc_type = {"pbe0"}; + mm.change_input("SCF Energy via TAMM", "xc_type", xc_type); + mm.change_input("SCF Energy via TAMM", "molecule_name", mol_name); + auto E = mm.at("SCF Energy via TAMM").run_as(aos, cs); + std::cout << "SCF Energy = " << E << " Hartree" << std::endl; + + REQUIRE(E == Catch::Approx(-74.81168986385825).margin(1.0e-6)); +} - REQUIRE(E == Catch::Approx(-74.81168986385825).margin(1.0e-6)); -} \ No newline at end of file +#endif \ No newline at end of file diff --git a/tests/cxx/integration_tests/update/diagonalization.cpp b/tests/cxx/integration_tests/update/diagonalization.cpp new file mode 100644 index 0000000..24cf377 --- /dev/null +++ b/tests/cxx/integration_tests/update/diagonalization.cpp @@ -0,0 +1,55 @@ +/* + * 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 "../integration_tests.hpp" + +using wf_t = simde::type::rscf_wf; +using pt = simde::UpdateGuess; +using orbital_index_set = typename wf_t::orbital_index_set_type; + +using simde::type::t_e_type; +using simde::type::v_en_type; + +TEST_CASE("Diagaonalization") { + auto mm = test_scf::load_modules(); + auto& mod = mm.at("Diagonalization Fock update"); + + auto aos = test_scf::h2_aos(); + auto h2 = test_scf::make_h2(); + simde::type::electron e; + orbital_index_set occs{0}; + + simde::type::fock f_e; + f_e.emplace_back(1.0, std::make_unique(e)); + f_e.emplace_back(1.0, std::make_unique(e, h2)); + + SECTION("No old guess") { + simde::type::tensor empty; + simde::type::cmos cmos(empty, aos, empty); + simde::type::rscf_wf core_guess(occs, cmos); + const auto& psi = mod.run_as(f_e, core_guess); + REQUIRE(psi.orbital_indices() == occs); + REQUIRE(psi.orbitals().from_space() == aos); + const auto& evals = psi.orbitals().diagonalized_matrix(); + using allocator_type = tensorwrapper::allocator::Eigen; + const auto& eval_buffer = allocator_type::rebind(evals.buffer()); + + const auto tol = 1E-6; + using Catch::Matchers::WithinAbs; + REQUIRE_THAT(eval_buffer.value()(0), WithinAbs(-1.25330893, tol)); + REQUIRE_THAT(eval_buffer.value()(1), WithinAbs(-0.47506974, tol)); + } +} \ No newline at end of file diff --git a/tests/cxx/test_scf.hpp b/tests/cxx/test_scf.hpp index 40f6391..4e5a4d6 100644 --- a/tests/cxx/test_scf.hpp +++ b/tests/cxx/test_scf.hpp @@ -100,20 +100,31 @@ inline auto h2_mos() { inline auto h2_density() { using density_type = simde::type::decomposable_e_density; typename density_type::value_type rho( - {{0.92501791, -0.54009707}, {-0.28540122, 1.7505162}}); + {{0.31980835, 0.31980835}, {0.31980835, 0.31980835}}); return density_type(rho, h2_mos()); } /// The Fock matrix consistent with h2_hamiltonian and h2_density +template inline auto h2_fock() { - simde::type::many_electrons es(2); + ElectronType es; + if constexpr(std::is_same_v) { + es = simde::type::many_electrons(2); + } + auto h2 = make_h2(); auto rho = h2_density(); simde::type::fock F; - F.emplace_back(1.0, std::make_unique(es)); - F.emplace_back(1.0, std::make_unique(es, h2)); - F.emplace_back(2.0, std::make_unique(es, rho)); - F.emplace_back(-1.0, std::make_unique(es, rho)); + using namespace chemist::qm_operator; + using t_type = Kinetic; + using v_type = Coulomb; + using j_type = Coulomb; + using k_type = Exchange; + + F.emplace_back(1.0, std::make_unique(es)); + F.emplace_back(1.0, std::make_unique(es, h2)); + F.emplace_back(2.0, std::make_unique(es, rho)); + F.emplace_back(-1.0, std::make_unique(es, rho)); return F; } diff --git a/tests/cxx/unit_tests/fock_operator/restricted.cpp b/tests/cxx/unit_tests/fock_operator/restricted.cpp index fa71d66..26f7a17 100644 --- a/tests/cxx/unit_tests/fock_operator/restricted.cpp +++ b/tests/cxx/unit_tests/fock_operator/restricted.cpp @@ -24,15 +24,55 @@ TEST_CASE("Restricted") { using density_type = simde::type::decomposable_e_density; using pt = simde::FockOperator; + auto H = test_scf::h2_hamiltonian(); + auto h2 = test_scf::make_h2(); + auto rho = test_scf::h2_density(); - auto& mod = mm.at("Restricted Fock Op"); + SECTION("Many-electron version") { + auto& mod = mm.at("Restricted Fock Op"); - SECTION("H2 Molecule") { - auto H = test_scf::h2_hamiltonian(); - auto rho = test_scf::h2_density(); + SECTION("No density") { + density_type rho; + auto F = mod.run_as(H, rho); + simde::type::many_electrons es(2); - auto F = mod.run_as(H, rho); - auto F_corr = test_scf::h2_fock(); - REQUIRE(F == F_corr); + simde::type::fock F_corr; + using simde::type::T_e_type; + using simde::type::V_en_type; + F_corr.emplace_back(1.0, std::make_unique(es)); + F_corr.emplace_back(1.0, std::make_unique(es, h2)); + REQUIRE(F == F_corr); + } + + SECTION("Density") { + auto rho = test_scf::h2_density(); + + auto F = mod.run_as(H, rho); + auto F_corr = test_scf::h2_fock(); + REQUIRE(F == F_corr); + } + } + SECTION("One-electron version") { + auto& mod = mm.at("Restricted One-Electron Fock Op"); + + SECTION(" No density") { + density_type rho; + auto f = mod.run_as(H, rho); + simde::type::fock f_corr; + simde::type::electron e; + using simde::type::t_e_type; + using simde::type::v_en_type; + f_corr.emplace_back(1.0, std::make_unique(e)); + f_corr.emplace_back(1.0, std::make_unique(e, h2)); + REQUIRE(f == f_corr); + } + + SECTION("Density") { + auto rho = test_scf::h2_density(); + + auto f = mod.run_as(H, rho); + auto f_corr = test_scf::h2_fock(); + REQUIRE(f == f_corr); + } } } \ No newline at end of file