diff --git a/CMakeLists.txt b/CMakeLists.txt index f413b57..3bc1679 100644 --- a/CMakeLists.txt +++ b/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}) diff --git a/src/scf/fock_operator/fock_operator.hpp b/src/scf/fock_operator/fock_operator.hpp index 313f248..c2327f1 100644 --- a/src/scf/fock_operator/fock_operator.hpp +++ b/src/scf/fock_operator/fock_operator.hpp @@ -22,28 +22,45 @@ namespace scf::fock_operator { template DECLARE_MODULE(Restricted); +template +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; - using e_e = Restricted; - using d_many = Restricted; - using d_e = Restricted; - + using e_many = Restricted; + using e_e = Restricted; + using d_many = Restricted; + using d_e = Restricted; + using e_ks_many = RKohnSham; + using e_ks_e = RKohnSham; + using d_ks_many = RKohnSham; + using d_ks_e = RKohnSham; 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"); + mm.add_module("AO Restricted Kohn-Sham Op"); + mm.add_module("Restricted Kohn-Sham Op"); + mm.add_module("AO Restricted One-Electron Kohn-Sham Op"); + mm.add_module("Restricted One-Electron Kohn-Sham Op"); } -#define EXTERN_RESTRICTED(density) \ +#define EXTERN_RESTRICTED(density) \ extern template struct Restricted; \ extern template struct Restricted +#define EXTERN_RKOHN_SHAM(density) \ + extern template struct RKohnSham; \ + extern template struct RKohnSham + 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 \ No newline at end of file diff --git a/src/scf/fock_operator/fock_operator_common.hpp b/src/scf/fock_operator/fock_operator_common.hpp new file mode 100644 index 0000000..84a7fee --- /dev/null +++ b/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 + +namespace scf::fock_operator { + +using namespace chemist::qm_operator; + +template +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>(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; + auto v = std::make_unique(get_e_(V_en), rhs); + m_pF_->emplace_back(1.0, std::move(v)); + } + +protected: + 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_; +}; +} // 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 79da8d3..ea16ba9 100644 --- a/src/scf/fock_operator/restricted.cpp +++ b/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 -class RestrictedVisitor : public chemist::qm_operator::OperatorVisitor { +class RestrictedVisitor : public FockOperatorCommon { +private: + using base_type = FockOperatorCommon; + 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>(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; - auto v = std::make_unique(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; 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)); + const auto& electrons = this->get_e_(V_ee); + auto j = std::make_unique(electrons, *m_prho_); + auto k = std::make_unique(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 - 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_; }; diff --git a/src/scf/fock_operator/rkohn_sham.cpp b/src/scf/fock_operator/rkohn_sham.cpp new file mode 100644 index 0000000..0d1ced3 --- /dev/null +++ b/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 +class RKohnShamVisitor : public FockOperatorCommon { +private: + using base_type = FockOperatorCommon; + +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; + using xc_type = ExchangeCorrelation; + + const auto& electrons = this->get_e_(V_ee); + auto j = std::make_unique(electrons, *m_prho_); + auto xc = std::make_unique(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 +TEMPLATED_MODULE_CTOR(RKohnSham, DensityType, ElectronType) { + using pt = simde::FockOperator; + satisfies_property_type(); + description(desc); + + add_input("XC Potential"); +} + +template +TEMPLATED_MODULE_RUN(RKohnSham, DensityType, ElectronType) { + using pt = simde::FockOperator; + + const auto& [H, rho] = pt::unwrap_inputs(inputs); + auto H_elec = H.electronic_hamiltonian(); + auto func = inputs.at("XC Potential").value(); + simde::type::fock F; + RKohnShamVisitor visitor(func, F, rho); + H_elec.visit(visitor); + + auto rv = results(); + return pt::wrap_results(rv, F); +} + +#define RESTRICTED(density) \ + template struct RKohnSham; \ + template struct RKohnSham + +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/matrix_builder/electronic_energy.cpp b/src/scf/matrix_builder/electronic_energy.cpp index b456549..aa3be00 100644 --- a/src/scf/matrix_builder/electronic_energy.cpp +++ b/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 using pt = simde::eval_braket; +template +using xc_pt = simde::eval_braket; + template using det_pt = simde::eval_braket; +template +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_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>(); add_submodule>("determinant driver"); + add_submodule>("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 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>(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>(O_ij); + temp("") = o("") * ci; + } if(i == 0) e = temp; else { e("") = e("") + temp(""); } diff --git a/src/scf/matrix_builder/matrix_builder.hpp b/src/scf/matrix_builder/matrix_builder.hpp index b8da0aa..3139a8c 100644 --- a/src/scf/matrix_builder/matrix_builder.hpp +++ b/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 \ No newline at end of file diff --git a/src/scf/matrix_builder/scf_integrals_driver.cpp b/src/scf/matrix_builder/scf_integrals_driver.cpp index 1314f34..99c6163 100644 --- a/src/scf/matrix_builder/scf_integrals_driver.cpp +++ b/src/scf/matrix_builder/scf_integrals_driver.cpp @@ -34,15 +34,16 @@ using simde::type::tensor; using pt = simde::aos_op_base_aos; using f_e_pt = simde::aos_f_e_aos; using rho_e_pt = simde::aos_rho_e_aos; +using xc_e_pt = simde::aos_xc_e_aos; class AODispatcher : public chemist::qm_operator::OperatorVisitor { private: using base_type = chemist::qm_operator::OperatorVisitor; public: - using f_e_type = simde::type::fock; - using rho_e_type = simde::type::rho_e; - + using f_e_type = simde::type::fock; + using rho_e_type = simde::type::rho_e; + using xc_e_type = simde::type::xc_e_type; using submods_type = pluginplay::type::submodule_map; AODispatcher(const aos& bra, const aos& ket, submodule_map& submods, @@ -50,6 +51,7 @@ class AODispatcher : public chemist::qm_operator::OperatorVisitor { base_type(false), m_pbra_(&bra), m_pket_(&ket), + m_evaluated_(false), m_psubmods_(&submods), m_ptensor_(&t) {} @@ -67,6 +69,13 @@ class AODispatcher : public chemist::qm_operator::OperatorVisitor { m_evaluated_ = true; } + void run(const xc_e_type& xc_e) { + chemist::braket::BraKet input(*m_pbra_, xc_e, *m_pket_); + const auto key = "XC Potential"; + *m_ptensor_ = m_psubmods_->at(key).run_as(input); + m_evaluated_ = true; + } + bool evaluated() const noexcept { return m_evaluated_; } private: @@ -81,8 +90,8 @@ MODULE_CTOR(SCFIntegralsDriver) { description(desc); satisfies_property_type(); add_submodule("Fundamental matrices"); - add_submodule("Fock matrix"); add_submodule("Density matrix"); + add_submodule("XC Potential"); } MODULE_RUN(SCFIntegralsDriver) { diff --git a/src/scf/scf_mm.cpp b/src/scf/scf_mm.cpp index bb9bf52..472dc4e 100644 --- a/src/scf/scf_mm.cpp +++ b/src/scf/scf_mm.cpp @@ -21,6 +21,7 @@ #include "matrix_builder/matrix_builder.hpp" #include "scf_modules.hpp" #include "update/update.hpp" +#include "xc/xc.hpp" #include namespace scf { @@ -32,6 +33,7 @@ void load_modules(pluginplay::ModuleManager& mm) { guess::load_modules(mm); matrix_builder::load_modules(mm); update::load_modules(mm); + xc::load_modules(mm); mm.add_module("Coulomb's Law"); @@ -44,6 +46,7 @@ void load_modules(pluginplay::ModuleManager& mm) { guess::set_defaults(mm); matrix_builder::set_defaults(mm); update::set_defaults(mm); + xc::set_defaults(mm); } } // namespace scf diff --git a/src/scf/xc/gauxc/basis_conversion.cpp b/src/scf/xc/gauxc/basis_conversion.cpp new file mode 100644 index 0000000..b8228cf --- /dev/null +++ b/src/scf/xc/gauxc/basis_conversion.cpp @@ -0,0 +1,75 @@ +/* + * 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 "gauxc.hpp" + +namespace scf::xc::gauxc { +namespace { + +const auto desc = R"( +Chemist AOBasisSet to GauXC BasisSet +==================================== + +TODO: Write me!! +)"; + +} + +// Module to convert NWX Basis -> GauXC Basis +MODULE_CTOR(BasisConversion) { + satisfies_property_type(); + description(desc); +} +MODULE_RUN(BasisConversion) { + const auto& [nwx_basis] = gauxc_basis_conversion_t::unwrap_inputs(inputs); + + GauXC::BasisSet gauxc_basis; + auto nshells = nwx_basis.n_shells(); + gauxc_basis.reserve(nshells); + for(decltype(nshells) i_sh = 0; i_sh < nshells; ++i_sh) { + const auto shell = nwx_basis.shell(i_sh); + const auto nprims = shell.n_primitives(); + const auto first_prim = shell.primitive(0); + const auto last_prim = shell.primitive(nprims - 1); + const int l = shell.l(); + const bool pure = shell.pure() == chemist::ShellType::pure; + + if(nprims > GauXC::detail::shell_nprim_max) { + std::stringstream ss; + ss << "GauXC received a shell with nprim = " << nprims + << "with NPRIM_MAX = " << GauXC::detail::shell_nprim_max + << ". Please reconfigure GauXC."; + throw std::runtime_error(ss.str()); + } + + GauXC::Shell::prim_array alpha, coeff; + GauXC::Shell::cart_array center = { + shell.center().x(), shell.center().y(), shell.center().z()}; + + std::copy(&first_prim.exponent(), &last_prim.exponent() + 1, + alpha.begin()); + std::copy(&first_prim.coefficient(), &last_prim.coefficient() + 1, + coeff.begin()); + + gauxc_basis.emplace_back( + GauXC::PrimSize(nprims), GauXC::AngularMomentum(l), + GauXC::SphericalType(pure), alpha, coeff, center); + } + + auto rv = results(); + return gauxc_basis_conversion_t::wrap_results(rv, gauxc_basis); +} +} // namespace scf::xc::gauxc diff --git a/src/scf/xc/gauxc/gauxc.cpp b/src/scf/xc/gauxc/gauxc.cpp new file mode 100644 index 0000000..e39efce --- /dev/null +++ b/src/scf/xc/gauxc/gauxc.cpp @@ -0,0 +1,57 @@ +/* + * 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 "gauxc.hpp" + +namespace scf::xc::gauxc { + +void load_modules(pluginplay::ModuleManager& mm) { + // Conversion Utilities + mm.add_module("GauXC Basis Converter"); + mm.add_module("GauXC Molecule Converter"); + + // Data geneation utilities + mm.add_module("GauXC Quadrature Batches"); + + // XC Integration + mm.add_module("GauXC Driver"); + mm.add_module("GauXC XC Potential"); + mm.add_module("GauXC XC Energy"); + + // sn-LinK Integration + mm.add_module("snLinK"); +} + +void set_defaults(pluginplay::ModuleManager& mm) { + // Data geneation utilities + mm.change_submod("GauXC Quadrature Batches", "GauXC Basis Converter", + "GauXC Basis Converter"); + mm.change_submod("GauXC Quadrature Batches", "GauXC Molecule Converter", + "GauXC Molecule Converter"); + + // XC Integration + mm.change_submod("GauXC Driver", "Quadrature Batches", + "GauXC Quadrature Batches"); + + mm.change_submod("GauXC XC Potential", "XC Driver", "GauXC Driver"); + mm.change_submod("GauXC XC Energy", "XC Driver", "GauXC Driver"); + + // sn-LinK Integration + mm.change_submod("snLinK", "Quadrature Batches", + "GauXC Quadrature Batches"); +} + +} // namespace scf::xc::gauxc \ No newline at end of file diff --git a/src/scf/xc/gauxc/gauxc.hpp b/src/scf/xc/gauxc/gauxc.hpp new file mode 100644 index 0000000..433e8cd --- /dev/null +++ b/src/scf/xc/gauxc/gauxc.hpp @@ -0,0 +1,41 @@ +/* + * 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. + */ + +#pragma once +#include "gauxc_property_types.hpp" +#include + +namespace scf::xc::gauxc { + +DECLARE_MODULE(GauXCDriver); +DECLARE_MODULE(XCPotential); +DECLARE_MODULE(XCEnergy); + +DECLARE_MODULE(snLinK); + +// Module to convert Chemist Basis -> GauXC Basis +DECLARE_MODULE(BasisConversion); + +// Module to convert Chemist AO Basis -> GauXC Molecule +DECLARE_MODULE(MoleculeConversion); + +// Modules to generate XC Quadrature Batches +DECLARE_MODULE(QuadratureBatches); + +void set_defaults(pluginplay::ModuleManager& mm); +void load_modules(pluginplay::ModuleManager& mm); + +} // namespace scf::xc::gauxc \ No newline at end of file diff --git a/src/scf/xc/gauxc/gauxc_driver.cpp b/src/scf/xc/gauxc/gauxc_driver.cpp new file mode 100644 index 0000000..8b55677 --- /dev/null +++ b/src/scf/xc/gauxc/gauxc_driver.cpp @@ -0,0 +1,87 @@ +/* + * 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 "gauxc.hpp" +#include "utilities.hpp" +#include +#include +#include + +namespace scf::xc::gauxc { + +using chemist::qm_operator::xc_functional; + +#define GAUXC_NWX_XC_PAIR(FUNC) \ + if(in == xc_functional::FUNC) return ExchCXX::Functional::FUNC + +ExchCXX::Functional nwx_to_exchcxx(xc_functional in) { + GAUXC_NWX_XC_PAIR(SVWN3); + else GAUXC_NWX_XC_PAIR(SVWN5); + else GAUXC_NWX_XC_PAIR(BLYP); + else GAUXC_NWX_XC_PAIR(B3LYP); + else GAUXC_NWX_XC_PAIR(PBE); + else GAUXC_NWX_XC_PAIR(revPBE); + else GAUXC_NWX_XC_PAIR(PBE0); + else throw std::out_of_range("No such functional"); +} + +#undef GAUXC_NWX_XC_PAIR + +// XC Integration +MODULE_CTOR(GauXCDriver) { + satisfies_property_type(); + + add_submodule("Quadrature Batches") + .set_description("Generate XC Quadrature Batches"); + + /// TODO: Replace with information from the runtime + add_input("On GPU").set_default(false); +} + +MODULE_RUN(GauXCDriver) { + const auto& [nwx_func, aos, P] = XCDriver::unwrap_inputs(inputs); + + auto P_eigen = tw_to_eigen(P); + + // Generate Partitioned Molecular Quadrature + const auto& lb = + submods.at("Quadrature Batches").run_as(aos); + + // Create XC functional + auto func_spec = nwx_to_exchcxx(nwx_func); + GauXC::functional_type func(ExchCXX::Backend::builtin, func_spec, + ExchCXX::Spin::Unpolarized); + + // Create XCIntegrator instance + auto on_gpu = inputs.at("On GPU").value(); + auto ex_space = + (on_gpu) ? GauXC::ExecutionSpace::Device : GauXC::ExecutionSpace::Host; + GauXC::XCIntegratorFactory integrator_factory( + ex_space, "Replicated", "Default", "Default", "Default"); + auto integrator = integrator_factory.get_instance(func, lb); + + // Do the XC integration + auto [EXC_float, VXC_eigen] = integrator.eval_exc_vxc(P_eigen); + + auto VXC = eigen_to_tw(VXC_eigen, get_runtime()); + simde::type::tensor EXC(EXC_float); + + // Wrap results + auto rv = results(); + return XCDriver::wrap_results(rv, EXC, VXC); +} + +} // namespace scf::xc::gauxc diff --git a/src/scf/xc/gauxc/gauxc_property_types.hpp b/src/scf/xc/gauxc/gauxc_property_types.hpp new file mode 100644 index 0000000..0be43d5 --- /dev/null +++ b/src/scf/xc/gauxc/gauxc_property_types.hpp @@ -0,0 +1,70 @@ +/* + * 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. + */ + +#pragma once +#include +#include +#include +#include +#include + +namespace scf::xc::gauxc { + +using gauxc_basis_conversion_t = + simde::Convert, simde::type::ao_basis_set>; + +using basis_to_gauxc_molecule_conversion_t = + simde::Convert; + +DECLARE_PROPERTY_TYPE(XCDriver); + +PROPERTY_TYPE_INPUTS(XCDriver) { + using functional_type = chemist::qm_operator::xc_functional; + using basis_type = simde::type::ao_basis_set; + using tensor_type = simde::type::tensor; + auto rv = pluginplay::declare_input() + .add_field("Functional Name") + .add_field("AO Basis") + .add_field("density"); + return rv; +} + +PROPERTY_TYPE_RESULTS(XCDriver) { + using tensor_type = simde::type::tensor; + auto rv = pluginplay::declare_result() + .add_field("XC Energy") + .add_field("XC Potential"); + + return rv; +} + +DECLARE_PROPERTY_TYPE(XCQuadratureBatches); + +PROPERTY_TYPE_INPUTS(XCQuadratureBatches) { + using basis_type = simde::type::ao_basis_set; + + auto rv = + pluginplay::declare_input().add_field("AO Basis"); + return rv; +} + +PROPERTY_TYPE_RESULTS(XCQuadratureBatches) { + using xc_quad_batch_type = GauXC::LoadBalancer; + return pluginplay::declare_result().add_field( + "Quadrature Batches"); +} + +} // namespace scf::xc::gauxc diff --git a/src/scf/xc/gauxc/gauxc_quadrature.cpp b/src/scf/xc/gauxc/gauxc_quadrature.cpp new file mode 100644 index 0000000..0714efe --- /dev/null +++ b/src/scf/xc/gauxc/gauxc_quadrature.cpp @@ -0,0 +1,126 @@ +// /* +// * 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 "gauxc.hpp" + +#include +#include +#include + +namespace scf::xc::gauxc { + +using xc_grid_type = GauXC::AtomicGridSizeDefault; + +/** + * All default grids are product quadratures consisting of NR Mura-Knowles + (MK) + * radial quadratures and NA Lebedev (L) spherical quadratures. + * + * Fine(Grid) + * - 75 (MK) x 302 (L) grid for each atom + * - Target Accuracy 1e-6 + * UltraFine(Grid) + * - 99 (MK) x 590 (L) grid for each atom + * - Target Accuracy 1e-8 + * SuperFine(Grid) + * - 175 (MK) x 974 (L) grid for atoms with Z <= 2 + * - 250 (MK) x 974 (L) grid for atoms with Z >= 3 + * - Target Accuracy 1e-10 + */ +static utilities::CaseInsensitiveMap atomic_grid_types = { + {"Fine", xc_grid_type::FineGrid}, + {"UltraFine", xc_grid_type::UltraFineGrid}, + {"SuperFine", xc_grid_type::SuperFineGrid}, + {"FineGrid", xc_grid_type::FineGrid}, + {"UltraFineGrid", xc_grid_type::UltraFineGrid}, + {"SuperFineGrid", xc_grid_type::SuperFineGrid}}; + +MODULE_CTOR(QuadratureBatches) { + satisfies_property_type(); + + add_submodule("GauXC Basis Converter") + .set_description("Converts NWX Basis -> GauXC Basis"); + add_submodule( + "GauXC Molecule Converter") + .set_description("Converts NWX Basis -> GauXC Molecule"); + + add_input("Grid Type").set_default("UltraFine"); + + add_input("Pruning Scheme") + .set_default(GauXC::PruningScheme::Unpruned); + + add_input("Radial Quadrature Type") + .set_default(GauXC::RadialQuad::MuraKnowles); + + /// TODO: Replace with information from the runtime + add_input("On GPU").set_default(false); +} + +MODULE_RUN(QuadratureBatches) { + const auto& [obs] = XCQuadratureBatches::unwrap_inputs(inputs); + + // Unpack module-specific inputs + auto grid_spec = + atomic_grid_types.at(inputs.at("Grid Type").value()); + + auto pruning_scheme = + inputs.at("Pruning Scheme").value(); + + auto radial_quad = + inputs.at("Radial Quadrature Type").value(); + + auto on_gpu = inputs.at("On GPU").value(); + + auto ex_space = (on_gpu) ? + GauXC::ExecutionSpace::Device : + GauXC::ExecutionSpace::Host; // TODO handle CUDA/HIP + auto comm = get_runtime().mpi_comm(); + std::shared_ptr gauxc_rt = +#ifdef GAUXC_ENABLE_DEVICE + (on_gpu) ? std::make_shared(comm, 0.8) : +#endif + std::make_shared(comm); + + // Convert NWX Molecule -> GauXC Molecule + auto& mol_submod = submods.at("GauXC Molecule Converter"); + const auto& gauxc_mol = + mol_submod.run_as(obs); + + // Convert NWX Basis -> GauXC Basis + auto& basis_submod = submods.at("GauXC Basis Converter"); + const auto& gauxc_basis = + basis_submod.run_as(obs); + + // Set up integration grid + const auto bsz = GauXC::BatchSize(512); + auto mol_grid = GauXC::MolGridFactory::create_default_molgrid( + gauxc_mol, pruning_scheme, bsz, radial_quad, grid_spec); + + // Generate LoadBalancer + GauXC::LoadBalancerFactory lb_factory(ex_space, "Default"); + auto lb = + lb_factory.get_instance(*gauxc_rt, gauxc_mol, mol_grid, gauxc_basis); + + // Apply Molecular Partition Weights + GauXC::MolecularWeightsFactory mw_factory( + ex_space, "Default", GauXC::MolecularWeightsSettings{}); + mw_factory.get_instance().modify_weights(lb); + + auto rv = results(); + return XCQuadratureBatches::wrap_results(rv, lb); +} + +} // namespace scf::xc::gauxc diff --git a/src/scf/xc/gauxc/molecule_conversion.cpp b/src/scf/xc/gauxc/molecule_conversion.cpp new file mode 100644 index 0000000..3923c87 --- /dev/null +++ b/src/scf/xc/gauxc/molecule_conversion.cpp @@ -0,0 +1,55 @@ +/* + * 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 "gauxc.hpp" + +namespace scf::xc::gauxc { +namespace { + +const auto desc = R"( +Chemist Molecule to GauXC Molecule +---------------------------------- +TODO: Write me!!! +)"; + +} + +using pt = basis_to_gauxc_molecule_conversion_t; + +// Module to convert NWX AO Basis -> GauXC Molecule +MODULE_CTOR(MoleculeConversion) { + satisfies_property_type(); + description(desc); +} + +MODULE_RUN(MoleculeConversion) { + const auto& [nwx_basis] = pt::unwrap_inputs(inputs); + + GauXC::Molecule gauxc_mol; + gauxc_mol.reserve(nwx_basis.size()); + for(const auto& center : nwx_basis) { + auto Z_optional = center.atomic_number(); + auto Z = Z_optional.has_value() ? Z_optional.value() : 0; + auto x = center.center().coord(0); + auto y = center.center().coord(1); + auto z = center.center().coord(2); + gauxc_mol.emplace_back(GauXC::AtomicNumber(Z), x, y, z); + } + auto rv = results(); + return pt::wrap_results(rv, gauxc_mol); +} + +} // namespace scf::xc::gauxc diff --git a/src/scf/xc/gauxc/snlink.cpp b/src/scf/xc/gauxc/snlink.cpp new file mode 100644 index 0000000..0ca2863 --- /dev/null +++ b/src/scf/xc/gauxc/snlink.cpp @@ -0,0 +1,73 @@ +/* + * 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 "gauxc.hpp" +#include "utilities.hpp" +#include +#include +#include + +namespace scf::xc::gauxc { + +using k_type = simde::aos_k_e_aos; + +// K Integration +MODULE_CTOR(snLinK) { + satisfies_property_type(); + + add_submodule("Quadrature Batches") + .set_description("Generate XC Quadrature Batches"); + + /// TODO: Replace with information from the runtime + add_input("On GPU").set_default(false); +} + +MODULE_RUN(snLinK) { + const auto&& [braket] = k_type::unwrap_inputs(inputs); + const auto& bra = braket.bra(); + const auto& P = braket.op().rhs_particle(); + const auto& ket = braket.ket(); + + if(bra != ket || P.basis_set() != bra) + throw std::runtime_error("GauXC: bra must be equal to ket"); + + // Extract P into replicated Eigen matrix + auto P_eigen = tw_to_eigen(P.value()); + + // Generate Partitioned Molecular Quadrature + auto& lb_mod = submods.at("Quadrature Batches"); + const auto& lb = lb_mod.run_as(bra.ao_basis_set()); + + // Create XCIntegrator instance + auto on_gpu = inputs.at("On GPU").value(); + auto ex_space = (on_gpu) ? + GauXC::ExecutionSpace::Device : + GauXC::ExecutionSpace::Host; // TODO handle CUDA/HIP + GauXC::functional_type func; // Dummy functional + GauXC::XCIntegratorFactory integrator_factory( + ex_space, "Replicated", "Default", "Default", "Default"); + auto integrator = integrator_factory.get_instance(func, lb); + + // Do the K integration + GauXC::IntegratorSettingsSNLinK sn_link_settings; + auto K_eigen = integrator.eval_exx(P_eigen, sn_link_settings); + auto K = eigen_to_tw(K_eigen, get_runtime()); + + // Wrap results + auto rv = results(); + return k_type::wrap_results(rv, K); +} + +} // namespace scf::xc::gauxc diff --git a/src/scf/xc/gauxc/utilities.hpp b/src/scf/xc/gauxc/utilities.hpp new file mode 100644 index 0000000..f686c90 --- /dev/null +++ b/src/scf/xc/gauxc/utilities.hpp @@ -0,0 +1,57 @@ +/* + * 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. + */ + +#pragma once +#include +#include + +namespace scf::xc::gauxc { + +template +auto tw_to_eigen(const tensorwrapper::Tensor& t) { + auto n = t.logical_layout().shape().as_smooth().extent(0); + auto m = t.logical_layout().shape().as_smooth().extent(1); + Eigen::Matrix t_eigen(n, m); + + const auto& rt = t.buffer().allocator().runtime(); + tensorwrapper::allocator::Eigen alloc(rt); + + auto t_vector = alloc.rebind(t.buffer()); + auto begin = t_vector.get_immutable_data(); + auto end = begin + (n * m); + std::copy(begin, end, t_eigen.data()); + return t_eigen; +} + +template +auto eigen_to_tw( + const Eigen::Matrix& t_eigen, + const parallelzone::runtime::RuntimeView& rt) { + tensorwrapper::allocator::Eigen alloc(rt); + auto n = static_cast(t_eigen.rows()); + auto m = static_cast(t_eigen.cols()); + + tensorwrapper::shape::Smooth shape{n, m}; + tensorwrapper::layout::Physical layout(shape); + auto pbuffer = alloc.allocate(layout); + for(std::size_t i = 0; i < n; ++i) { + for(std::size_t j = 0; j < m; ++j) + pbuffer->set_elem({i, j}, t_eigen(i, j)); + } + return tensorwrapper::Tensor(std::move(shape), std::move(pbuffer)); +} + +} // namespace scf::xc::gauxc \ No newline at end of file diff --git a/src/scf/xc/gauxc/xc_energy.cpp b/src/scf/xc/gauxc/xc_energy.cpp new file mode 100644 index 0000000..792e4f4 --- /dev/null +++ b/src/scf/xc/gauxc/xc_energy.cpp @@ -0,0 +1,55 @@ +/* + * 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 "gauxc.hpp" + +namespace scf::xc::gauxc { + +using XC_e_t = simde::type::XC_e_type; + +template +using pt = simde::eval_braket; + +MODULE_CTOR(XCEnergy) { + using wf_type = simde::type::rscf_wf; + satisfies_property_type>(); + + add_submodule("XC Driver"); +} + +MODULE_RUN(XCEnergy) { + using wf_type = simde::type::rscf_wf; + const auto& [braket] = pt::unwrap_inputs(inputs); + + const auto& bra_wf = braket.bra(); + const auto& xc_op = braket.op(); + const auto& ket_wf = braket.ket(); + + if(bra_wf != ket_wf) + throw std::runtime_error("Expected the same basis set"); + + const auto func = xc_op.functional_name(); + const auto& P = xc_op.rhs_particle(); + const auto& aos = P.basis_set().ao_basis_set(); + + auto& driver = submods.at("XC Driver"); + const auto& [exc, vxc] = driver.run_as(func, aos, P.value()); + + auto rv = results(); + return pt::wrap_results(rv, exc); +} + +} // namespace scf::xc::gauxc \ No newline at end of file diff --git a/src/scf/xc/gauxc/xc_potential.cpp b/src/scf/xc/gauxc/xc_potential.cpp new file mode 100644 index 0000000..15199a1 --- /dev/null +++ b/src/scf/xc/gauxc/xc_potential.cpp @@ -0,0 +1,57 @@ +/* + * 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 "gauxc.hpp" + +namespace scf::xc::gauxc { +namespace { +const auto desc = R"( +Exchange-Correlation Potential via Driver +========================================= +)"; +} + +using pt = simde::aos_xc_e_aos; + +MODULE_CTOR(XCPotential) { + satisfies_property_type(); + description(desc); + + add_submodule("XC Driver"); +} + +MODULE_RUN(XCPotential) { + const auto& [braket] = pt::unwrap_inputs(inputs); + + const auto& bra_aos = braket.bra(); + const auto& xc_op = braket.op(); + const auto& ket_aos = braket.ket(); + + if(bra_aos != ket_aos) + throw std::runtime_error("Expected the same basis set!"); + + const auto func = xc_op.functional_name(); + const auto& P = xc_op.rhs_particle(); + const auto& aos = bra_aos.ao_basis_set(); + + auto& driver = submods.at("XC Driver"); + const auto& [exc, vxc] = driver.run_as(func, aos, P.value()); + + auto rv = results(); + return pt::wrap_results(rv, vxc); +} + +} // namespace scf::xc::gauxc \ No newline at end of file diff --git a/src/scf/xc/xc.cpp b/src/scf/xc/xc.cpp new file mode 100644 index 0000000..e74036e --- /dev/null +++ b/src/scf/xc/xc.cpp @@ -0,0 +1,24 @@ +/* + * 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 "gauxc/gauxc.hpp" +#include "xc.hpp" + +namespace scf::xc { +void load_modules(pluginplay::ModuleManager& mm) { gauxc::load_modules(mm); } + +void set_defaults(pluginplay::ModuleManager& mm) { gauxc::set_defaults(mm); } +} // namespace scf::xc diff --git a/src/scf/xc/xc.hpp b/src/scf/xc/xc.hpp new file mode 100644 index 0000000..74ade90 --- /dev/null +++ b/src/scf/xc/xc.hpp @@ -0,0 +1,25 @@ +/* + * 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. + */ + +#pragma once +#include + +namespace scf::xc { + +void set_defaults(pluginplay::ModuleManager& mm); +void load_modules(pluginplay::ModuleManager& mm); + +} // namespace scf::xc diff --git a/tests/cxx/integration_tests/driver/scf_driver.cpp b/tests/cxx/integration_tests/driver/scf_driver.cpp index 3d6bfd3..d7a562f 100644 --- a/tests/cxx/integration_tests/driver/scf_driver.cpp +++ b/tests/cxx/integration_tests/driver/scf_driver.cpp @@ -17,22 +17,41 @@ #include "../integration_tests.hpp" using pt = simde::AOEnergy; +using tensorwrapper::operations::approximately_equal; TEMPLATE_LIST_TEST_CASE("SCFDriver", "", test_scf::float_types) { using float_type = TestType; auto mm = test_scf::load_modules(); - auto h2 = test_scf::make_h2(); - auto aos = test_scf::h2_aos().ao_basis_set(); tensorwrapper::allocator::Eigen alloc(mm.get_runtime()); tensorwrapper::shape::Smooth shape_corr{}; auto pcorr = alloc.allocate(tensorwrapper::layout::Physical(shape_corr)); - using tensorwrapper::operations::approximately_equal; - pcorr->set_elem({}, -1.1167592336); - tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); - const auto e = mm.template run_as("SCF Driver", aos, h2); - REQUIRE(approximately_equal(corr, e, 1E-6)); + SECTION("H2") { + auto h2 = test_scf::make_h2(); + auto aos = test_scf::h2_aos().ao_basis_set(); + + SECTION("SCF") { + pcorr->set_elem({}, -1.1167592336); + simde::type::tensor corr(shape_corr, std::move(pcorr)); + const auto e = mm.template run_as("SCF Driver", aos, h2); + REQUIRE(approximately_equal(corr, e, 1E-6)); + } + SECTION("DFT") { + auto func = chemist::qm_operator::xc_functional::PBE; + const auto RKS_op = "Restricted Kohn-Sham Op"; + const auto rks_op = "Restricted One-Electron Kohn-Sham Op"; + mm.change_input(RKS_op, "XC Potential", func); + mm.change_input(rks_op, "XC Potential", func); + mm.change_submod("Loop", "One-electron Fock operator", rks_op); + mm.change_submod("Loop", "Fock operator", RKS_op); + mm.change_submod("Core guess", "Build Fock Operator", rks_op); + const auto e = mm.template run_as("SCF Driver", aos, h2); + pcorr->set_elem({}, -1.15207); + simde::type::tensor corr(shape_corr, std::move(pcorr)); + REQUIRE(approximately_equal(corr, e, 1E-5)); + } + } SECTION("H2 Dimer") { simde::type::nucleus h0("H", 1ul, 1836.15, 0.0, 0.0, 0.0); @@ -45,7 +64,8 @@ TEMPLATE_LIST_TEST_CASE("SCFDriver", "", test_scf::float_types) { simde::type::chemical_system h2_dimer_sys(h2_dimer_mol); const auto e = mm.template run_as("SCF Driver", ao_bs, h2_dimer_sys); - alloc.rebind(corr.buffer()).set_elem({}, -2.2260535919670001); + pcorr->set_elem({}, -2.2260535919670001); + simde::type::tensor corr(shape_corr, std::move(pcorr)); REQUIRE(approximately_equal(corr, e, 1E-6)); } } \ No newline at end of file diff --git a/tests/cxx/integration_tests/driver/scf_loop.cpp b/tests/cxx/integration_tests/driver/scf_loop.cpp index e7d2473..f786021 100644 --- a/tests/cxx/integration_tests/driver/scf_loop.cpp +++ b/tests/cxx/integration_tests/driver/scf_loop.cpp @@ -39,6 +39,7 @@ TEMPLATE_LIST_TEST_CASE("SCFLoop", "", test_scf::float_types) { SECTION("H2") { wf_type psi0(index_set{0}, test_scf::h2_cmos()); + auto H = test_scf::h2_hamiltonian(); chemist::braket::BraKet H_00(psi0, H, psi0); diff --git a/tests/cxx/integration_tests/integration_tests.hpp b/tests/cxx/integration_tests/integration_tests.hpp index 1ef8cec..a71084c 100644 --- a/tests/cxx/integration_tests/integration_tests.hpp +++ b/tests/cxx/integration_tests/integration_tests.hpp @@ -39,9 +39,6 @@ pluginplay::ModuleManager load_modules() { mm.change_submod(ao_driver_key, "Fundamental matrices", "AO integral driver"); - mm.change_submod("Fock Matrix Builder", "Two center evaluator", - "AO integral driver"); - mm.change_submod("Diagonalization Fock update", "Overlap matrix builder", "Overlap"); diff --git a/tests/cxx/integration_tests/matrix_builder/electronic_energy.cpp b/tests/cxx/integration_tests/matrix_builder/electronic_energy.cpp index 8f869a3..364044b 100644 --- a/tests/cxx/integration_tests/matrix_builder/electronic_energy.cpp +++ b/tests/cxx/integration_tests/matrix_builder/electronic_energy.cpp @@ -44,14 +44,28 @@ TEMPLATE_LIST_TEST_CASE("ElectronicEnergy", "", test_scf::float_types) { auto rho = test_scf::h2_density(); simde::type::J_e_type J_e(es, rho); - simde::type::K_e_type K_e(es, rho); - simde::type::electronic_hamiltonian H_e(T_e * 2.0 + V_en * 2.0 + J_e * 2.0 - - K_e); - chemist::braket::BraKet braket(psi, H_e, psi); + SECTION("RHF") { + simde::type::K_e_type K_e(es, rho); + simde::type::electronic_hamiltonian H_e(T_e * 2.0 + V_en * 2.0 + + J_e * 2.0 - K_e); + chemist::braket::BraKet braket(psi, H_e, psi); - const auto& E_elec = mod.template run_as>(braket); + const auto& E_elec = mod.template run_as>(braket); - pcorr->set_elem({}, -1.90066758625308307); - tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); - REQUIRE(approximately_equal(corr, E_elec, 1E-6)); + pcorr->set_elem({}, -1.90066758625308307); + tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); + REQUIRE(approximately_equal(corr, E_elec, 1E-6)); + } + SECTION("RKS") { + if constexpr(std::is_same_v) { + auto func = chemist::qm_operator::xc_functional::PBE; + simde::type::XC_e_type XC_e(func, es, rho); + simde::type::electronic_hamiltonian H_e(T_e * 2.0 + V_en * 2.0 + + J_e * 2.0 + XC_e); + chemist::braket::BraKet braket(psi, H_e, psi); + const auto& E_elec = mod.template run_as>(braket); + tensorwrapper::Tensor corr(-1.90692); + REQUIRE(approximately_equal(corr, E_elec, 1E-5)); + } + } } \ No newline at end of file diff --git a/tests/cxx/integration_tests/matrix_builder/scf_integrals_driver.cpp b/tests/cxx/integration_tests/matrix_builder/scf_integrals_driver.cpp index d3a0b46..7bbc767 100644 --- a/tests/cxx/integration_tests/matrix_builder/scf_integrals_driver.cpp +++ b/tests/cxx/integration_tests/matrix_builder/scf_integrals_driver.cpp @@ -96,4 +96,16 @@ TEMPLATE_LIST_TEST_CASE("SCFIntegralsDriver", "", test_scf::float_types) { // const auto& F_corr = fmod.run_as(braket); // compare_matrices(F, F_corr); // } + + SECTION("Calling XC Potential") { + auto func = chemist::qm_operator::xc_functional::PBE0; + auto rho = test_scf::h2_density(); + const auto& aos = rho.basis_set(); + simde::type::xc_e_type xc_op(func, e, rho); + chemist::braket::BraKet xc_ij(aos, xc_op, aos); + erased_type copy_braket(xc_ij); + auto vxc = mod.template run_as(copy_braket); + simde::type::tensor corr{{-0.357302, -0.23347}, {-0.23347, -0.357302}}; + REQUIRE(approximately_equal(vxc, corr, 1E-5)); + } } \ No newline at end of file diff --git a/tests/cxx/test_scf.hpp b/tests/cxx/test_scf.hpp index 5507f63..508c5d4 100644 --- a/tests/cxx/test_scf.hpp +++ b/tests/cxx/test_scf.hpp @@ -122,7 +122,7 @@ inline auto he_basis(NucleiType& heliums) { center_type coords(hi.x(), hi.y(), hi.z()); contracted_gaussian_type he_cg(he_coefs.begin(), he_coefs.end(), he_exps.begin(), he_exps.end(), coords); - atomic_basis_type he_basis("STO-3G", 1, coords); + atomic_basis_type he_basis("STO-3G", 2, coords); he_basis.add_shell(cartesian, l0, he_cg); rv.add_center(he_basis); } @@ -232,6 +232,19 @@ inline auto h2_density() { return density_type(std::move(t), h2_mos()); } +template +inline auto he_density() { + using density_type = simde::type::decomposable_e_density; + using tensor_type = typename density_type::value_type; + using allocator_type = tensorwrapper::allocator::Eigen; + allocator_type alloc(parallelzone::runtime::RuntimeView{}); + tensorwrapper::shape::Smooth shape{1, 1}; + tensorwrapper::layout::Physical l(shape); + auto pbuffer = alloc.construct(l, 1.000); + tensor_type t(shape, std::move(pbuffer)); + return density_type(std::move(t), he_mos()); +} + /// The Fock matrix consistent with h2_hamiltonian and h2_density template inline auto h2_fock() { diff --git a/tests/cxx/unit_tests/xc/gauxc/basis_conversion.cpp b/tests/cxx/unit_tests/xc/gauxc/basis_conversion.cpp new file mode 100644 index 0000000..05dfc28 --- /dev/null +++ b/tests/cxx/unit_tests/xc/gauxc/basis_conversion.cpp @@ -0,0 +1,49 @@ +/* + * 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 "../test_scf.hpp" +#include +#include + +using pt = scf::xc::gauxc::gauxc_basis_conversion_t; + +TEST_CASE("BasisConversion") { + pluginplay::ModuleManager mm; + scf::load_modules(mm); + auto& mod = mm.at("GauXC Basis Converter"); + + using nuclei_type = simde::type::nuclei; + + SECTION("he") { + auto mol = test_scf::make_he(); + auto bs = test_scf::he_basis(mol); + + auto rv = mod.run_as(bs); + REQUIRE(rv.nbf() == 1); + REQUIRE(rv.nshells() == 1); + REQUIRE(rv.max_l() == 0); + } + + SECTION("h2") { + auto mol = test_scf::make_h2(); + auto bs = test_scf::h_basis(mol); + + auto rv = mod.run_as(bs); + REQUIRE(rv.nbf() == 2); + REQUIRE(rv.nshells() == 2); + REQUIRE(rv.max_l() == 0); + } +} \ No newline at end of file diff --git a/tests/cxx/unit_tests/xc/gauxc/gauxc_driver.cpp b/tests/cxx/unit_tests/xc/gauxc/gauxc_driver.cpp new file mode 100644 index 0000000..5c18d6a --- /dev/null +++ b/tests/cxx/unit_tests/xc/gauxc/gauxc_driver.cpp @@ -0,0 +1,62 @@ +/* + * 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 "../../../test_scf.hpp" +#include +#include +#include + +using namespace scf; + +using pt = scf::xc::gauxc::XCDriver; +using tensorwrapper::operations::approximately_equal; + +TEST_CASE("GauXCDriver") { + pluginplay::ModuleManager mm; + scf::load_modules(mm); + auto& mod = mm.at("GauXC Driver"); + + auto func = chemist::qm_operator::xc_functional::PBE0; + + SECTION("h2") { + auto rho = test_scf::h2_density(); + const auto& aos = rho.basis_set().ao_basis_set(); + auto [exc, vxc] = mod.run_as(func, aos, rho.value()); + simde::type::tensor exc_corr(-0.587164); + REQUIRE(approximately_equal(exc_corr, exc, 1E-5)); + + simde::type::tensor vxc_corr{{-0.357302, -0.23347}, + {-0.23347, -0.357302}}; + REQUIRE(approximately_equal(vxc, vxc_corr, 1E-5)); + } + + SECTION("he") { + auto rho = test_scf::he_density(); + const auto& aos = rho.basis_set().ao_basis_set(); + + auto [exc, vxc] = mod.run_as(func, aos, rho.value()); + simde::type::tensor exc_corr(-0.819986); + REQUIRE(approximately_equal(exc_corr, exc, 1E-5)); + + tensorwrapper::allocator::Eigen alloc(mm.get_runtime()); + tensorwrapper::shape::Smooth shape_corr{1, 1}; + auto pcorr = + alloc.allocate(tensorwrapper::layout::Physical(shape_corr)); + pcorr->set_elem({0, 0}, -0.526535); + simde::type::tensor vxc_corr(shape_corr, std::move(pcorr)); + REQUIRE(approximately_equal(vxc, vxc_corr, 1E-5)); + } +} diff --git a/tests/cxx/unit_tests/xc/gauxc/molecule_conversion.cpp b/tests/cxx/unit_tests/xc/gauxc/molecule_conversion.cpp new file mode 100644 index 0000000..9d5c723 --- /dev/null +++ b/tests/cxx/unit_tests/xc/gauxc/molecule_conversion.cpp @@ -0,0 +1,47 @@ +/* + * 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 "../test_scf.hpp" +#include +#include + +using pt = scf::xc::gauxc::basis_to_gauxc_molecule_conversion_t; + +TEST_CASE("MoleculeConversion") { + pluginplay::ModuleManager mm; + scf::load_modules(mm); + auto& mod = mm.at("GauXC Molecule Converter"); + + using nuclei_type = simde::type::nuclei; + + SECTION("he") { + auto mol = test_scf::make_he(); + auto bs = test_scf::he_basis(mol); + + auto rv = mod.run_as(bs); + REQUIRE(rv.natoms() == 1); + REQUIRE(rv.maxZ() == 2); + } + + SECTION("h2") { + auto mol = test_scf::make_h2(); + auto bs = test_scf::h_basis(mol); + + auto rv = mod.run_as(bs); + REQUIRE(rv.natoms() == 2); + REQUIRE(rv.maxZ() == 1); + } +} \ No newline at end of file diff --git a/tests/cxx/unit_tests/xc/gauxc/snlink.cpp b/tests/cxx/unit_tests/xc/gauxc/snlink.cpp new file mode 100644 index 0000000..c3b59fd --- /dev/null +++ b/tests/cxx/unit_tests/xc/gauxc/snlink.cpp @@ -0,0 +1,39 @@ +/* + * 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 "../../../test_scf.hpp" + +using k_pt = simde::aos_k_e_aos; +using tensorwrapper::operations::approximately_equal; + +TEST_CASE("snLinK") { + pluginplay::ModuleManager mm; + scf::load_modules(mm); + auto& mod = mm.at("snLinK"); + + SECTION("h2") { + auto rho = test_scf::h2_density(); + const auto& aos = rho.basis_set(); + simde::type::electron e; + simde::type::k_e_type k(e, rho); + chemist::braket::BraKet k_ij(aos, k, aos); + const auto& K = mod.run_as(k_ij); + + simde::type::tensor corr_k( + {{0.627264, 0.561828}, {0.561828, 0.627264}}); + REQUIRE(approximately_equal(K, corr_k, 1E-5)); + } +} \ No newline at end of file diff --git a/tests/cxx/unit_tests/xc/gauxc/xc_energy.cpp b/tests/cxx/unit_tests/xc/gauxc/xc_energy.cpp new file mode 100644 index 0000000..a6a3378 --- /dev/null +++ b/tests/cxx/unit_tests/xc/gauxc/xc_energy.cpp @@ -0,0 +1,62 @@ +/* + * 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 "../../../test_scf.hpp" +#include +#include + +using namespace scf; + +using XC_e_t = simde::type::XC_e_type; + +template +using pt = simde::eval_braket; + +using tensorwrapper::operations::approximately_equal; + +TEST_CASE("XCEnergy") { + using wf_type = simde::type::rscf_wf; + using index_set = typename wf_type::orbital_index_set_type; + using float_type = double; + + pluginplay::ModuleManager mm; + scf::load_modules(mm); + auto& mod = mm.at("GauXC XC Energy"); + + auto func = chemist::qm_operator::xc_functional::PBE0; + + SECTION("h2") { + wf_type psi(index_set{0}, test_scf::h2_cmos()); + auto rho = test_scf::h2_density(); + simde::type::many_electrons es{2}; + XC_e_t xc_op(func, es, rho); + chemist::braket::BraKet xc_ij(psi, xc_op, psi); + auto exc = mod.run_as>(xc_ij); + simde::type::tensor corr(-0.587164); + REQUIRE(approximately_equal(exc, corr, 1E-5)); + } + + SECTION("he") { + wf_type psi(index_set{0}, test_scf::he_cmos()); + auto rho = test_scf::he_density(); + simde::type::many_electrons es{2}; + XC_e_t xc_op(func, es, rho); + chemist::braket::BraKet xc_ij(psi, xc_op, psi); + auto exc = mod.run_as>(xc_ij); + simde::type::tensor corr(-0.819986); + REQUIRE(approximately_equal(exc, corr, 1E-5)); + } +} \ No newline at end of file diff --git a/tests/cxx/unit_tests/xc/gauxc/xc_potential.cpp b/tests/cxx/unit_tests/xc/gauxc/xc_potential.cpp new file mode 100644 index 0000000..6aa0bde --- /dev/null +++ b/tests/cxx/unit_tests/xc/gauxc/xc_potential.cpp @@ -0,0 +1,60 @@ +/* + * 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 "../../../test_scf.hpp" +#include +#include + +using namespace scf; + +using pt = simde::aos_xc_e_aos; +using tensorwrapper::operations::approximately_equal; + +TEST_CASE("XCPotential") { + pluginplay::ModuleManager mm; + scf::load_modules(mm); + auto& mod = mm.at("GauXC XC Potential"); + + simde::type::electron e; + auto func = chemist::qm_operator::xc_functional::PBE0; + + SECTION("h2") { + auto rho = test_scf::h2_density(); + const auto& aos = rho.basis_set(); + simde::type::xc_e_type xc_op(func, e, rho); + chemist::braket::BraKet xc_ij(aos, xc_op, aos); + auto vxc = mod.run_as(xc_ij); + + simde::type::tensor corr{{-0.357302, -0.23347}, {-0.23347, -0.357302}}; + REQUIRE(approximately_equal(vxc, corr, 1E-5)); + } + + SECTION("he") { + auto rho = test_scf::he_density(); + const auto& aos = rho.basis_set(); + simde::type::xc_e_type xc_op(func, e, rho); + chemist::braket::BraKet xc_ij(aos, xc_op, aos); + auto vxc = mod.run_as(xc_ij); + + tensorwrapper::allocator::Eigen alloc(mm.get_runtime()); + tensorwrapper::shape::Smooth shape_corr{1, 1}; + auto pcorr = + alloc.allocate(tensorwrapper::layout::Physical(shape_corr)); + pcorr->set_elem({0, 0}, -0.526535); + simde::type::tensor corr(shape_corr, std::move(pcorr)); + REQUIRE(approximately_equal(vxc, corr, 1E-5)); + } +} \ No newline at end of file