diff --git a/src/scf/driver/driver.hpp b/src/scf/driver/driver.hpp new file mode 100644 index 0000000..869ade9 --- /dev/null +++ b/src/scf/driver/driver.hpp @@ -0,0 +1,43 @@ +/* + * 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::driver { + +DECLARE_MODULE(SCFDriver); +DECLARE_MODULE(SCFLoop); + +inline void load_modules(pluginplay::ModuleManager& mm) { + mm.add_module("SCF Driver"); + mm.add_module("Loop"); +} + +inline void set_defaults(pluginplay::ModuleManager& mm) { + mm.change_submod("Loop", "Electronic energy", "Electronic energy"); + mm.change_submod("Loop", "Density matrix", "Density matrix builder"); + mm.change_submod("Loop", "Guess update", "Diagonalization Fock update"); + mm.change_submod("Loop", "One-electron Fock operator", + "Restricted One-Electron Fock op"); + mm.change_submod("Loop", "Fock operator", "Restricted Fock Op"); + mm.change_submod("Loop", "Charge-charge", "Coulomb's Law"); + + mm.change_submod("SCF Driver", "Guess", "Core guess"); + mm.change_submod("SCF Driver", "Optimizer", "Loop"); +} + +} // namespace scf::driver \ No newline at end of file diff --git a/src/scf/driver/scf_driver.cpp b/src/scf/driver/scf_driver.cpp new file mode 100644 index 0000000..116ee06 --- /dev/null +++ b/src/scf/driver/scf_driver.cpp @@ -0,0 +1,75 @@ +/* + * 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 "driver.hpp" + +namespace scf::driver { +namespace { + +const auto desc = R"( +)"; + +} + +using simde::type::hamiltonian; + +using pt = simde::AOEnergy; + +using ham_pt = simde::Convert; + +template +using guess_pt = simde::InitialGuess; + +template +using braket_t = chemist::braket::BraKet; + +template +using egy_pt = simde::EvaluateBraKet>; + +template +using opt_pt = simde::Optimize, WfType>; + +MODULE_CTOR(SCFDriver) { + using wf_type = simde::type::rscf_wf; + + description(desc); + satisfies_property_type(); + add_submodule("Hamiltonian"); + add_submodule>("Guess"); + add_submodule>("Optimizer"); +} + +MODULE_RUN(SCFDriver) { + using wf_type = simde::type::rscf_wf; + + const auto& [ao_params, sys] = pt::unwrap_inputs(inputs); + simde::type::aos aos(ao_params); + + auto& ham_mod = submods.at("Hamiltonian"); + const auto& H = ham_mod.run_as(sys); + + auto& guess_mod = submods.at("Guess"); + const auto& Psi0 = guess_mod.run_as>(H, aos); + + auto& opt_mod = submods.at("Optimizer"); + braket_t H_00(Psi0, H, Psi0); + const auto&& [e, Psi] = opt_mod.run_as>(H_00, Psi0); + + auto rv = results(); + return pt::wrap_results(rv, e); +} + +} // namespace scf::driver \ No newline at end of file diff --git a/src/scf/driver/scf_loop.cpp b/src/scf/driver/scf_loop.cpp new file mode 100644 index 0000000..c005c96 --- /dev/null +++ b/src/scf/driver/scf_loop.cpp @@ -0,0 +1,161 @@ +/* + * 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 "driver.hpp" + +namespace scf::driver { +namespace { + +const auto desc = R"( +)"; + +} + +using simde::type::electronic_hamiltonian; +using simde::type::hamiltonian; +using simde::type::op_base_type; + +template +using egy_pt = simde::eval_braket; + +template +using elec_egy_pt = simde::eval_braket; + +template +using pt = simde::Optimize, WfType>; + +template +using update_pt = simde::UpdateGuess; + +using density_t = simde::type::decomposable_e_density; + +using fock_pt = simde::FockOperator; + +using density_pt = simde::aos_rho_e_aos; + +using v_nn_pt = simde::charge_charge_interaction; + +struct GrabNuclear : chemist::qm_operator::OperatorVisitor { + using V_nn_type = simde::type::V_nn_type; + + GrabNuclear() : chemist::qm_operator::OperatorVisitor(false) {} + + void run(const V_nn_type& V_nn) { m_pv = &V_nn; } + + const V_nn_type* m_pv; +}; + +MODULE_CTOR(SCFLoop) { + using wf_type = simde::type::rscf_wf; + description(desc); + satisfies_property_type>(); + + add_submodule>("Electronic energy"); + add_submodule("Density matrix"); + add_submodule>("Guess update"); + add_submodule("One-electron Fock operator"); + add_submodule("Fock operator"); + add_submodule("Charge-charge"); +} + +MODULE_RUN(SCFLoop) { + using wf_type = simde::type::rscf_wf; + using density_op_type = simde::type::rho_e; + const auto&& [braket, psi0] = pt::unwrap_inputs(inputs); + // TODO: Assert bra == ket == psi0 + const auto& H = braket.op(); + const auto& H_elec = H.electronic_hamiltonian(); + const auto& H_core = H_elec.core_hamiltonian(); + const auto& aos = psi0.orbitals().from_space(); + + auto& egy_mod = submods.at("Electronic energy"); + auto& density_mod = submods.at("Density matrix"); + auto& update_mod = submods.at("Guess update"); + auto& fock_mod = submods.at("One-electron Fock operator"); + auto& Fock_mod = submods.at("Fock operator"); + auto& V_nn_mod = submods.at("Charge-charge"); + + // Step 1: Nuclear-nuclear repulsion + GrabNuclear visitor; + H.visit(visitor); + // TODO: Clean up charges class to make this easier... + const auto& V_nn = *visitor.m_pv; + const auto n_lhs = V_nn.lhs_particle().as_nuclei(); + const auto qs_lhs_view = n_lhs.charges(); + const auto n_rhs = V_nn.rhs_particle().as_nuclei(); + const auto qs_rhs_view = n_rhs.charges(); + simde::type::charges qs_lhs; + simde::type::charges qs_rhs; + for(const auto q_i : qs_lhs_view) { + qs_lhs.push_back(q_i.as_point_charge()); + } + for(const auto q_i : qs_rhs_view) { + qs_rhs.push_back(q_i.as_point_charge()); + } + auto e_nuclear = V_nn_mod.run_as(qs_lhs, qs_rhs); + + wf_type psi_old = psi0; + double e_old = 0; + const unsigned int max_iter = 3; + unsigned int iter = 0; + + while(iter < max_iter) { + // Step 2: Build old density + density_op_type rho_hat(psi_old.orbitals(), psi_old.occupations()); + chemist::braket::BraKet P_mn(aos, rho_hat, aos); + const auto& P = density_mod.run_as(P_mn); + density_t rho_old(P, psi_old.orbitals()); + + // Step 3: Old density is used to create the new Fock operator + // TODO: Make easier to go from many-electron to one-electron + // TODO: template fock_pt on Hamiltonian type and only pass H_elec + const auto& f_new = fock_mod.run_as(H, rho_old); + const auto& F_new = Fock_mod.run_as(H, rho_old); + + // Step 4: New Fock operator is used to compute the new wavefunction + auto new_psi = update_mod.run_as>(f_new, psi_old); + const auto& new_cmos = new_psi.orbitals(); + const auto& new_evals = new_cmos.diagonalized_matrix(); + const auto& new_c = new_cmos.transform(); + + // Step 5: New electronic energy + // Step 5a: New Fock operator to new electronic Hamiltonian + // TODO: Should just be H_core + F; + electronic_hamiltonian H_new; + for(std::size_t i = 0; i < H_core.size(); ++i) + H_new.emplace_back(H_core.coefficient(i), + H_core.get_operator(i).clone()); + for(std::size_t i = 0; i < F_new.size(); ++i) + H_new.emplace_back(F_new.coefficient(i), + F_new.get_operator(i).clone()); + + // Step 5b: New electronic hamiltonian to new electronic energy + chemist::braket::BraKet H_00(new_psi, H_new, new_psi); + auto e_new = egy_mod.run_as>(H_00); + + // Step 6: Converged? + // TODO: gradient and energy differences + + // Step 7: Not converged so reset + e_old = e_new; + psi_old = new_psi; + ++iter; + } + auto rv = results(); + return pt::wrap_results(rv, e_old + e_nuclear, psi_old); +} + +} // namespace scf::driver \ No newline at end of file diff --git a/src/scf/eigen_solver/eigen_generalized.cpp b/src/scf/eigen_solver/eigen_generalized.cpp index ec47bed..941ed6a 100644 --- a/src/scf/eigen_solver/eigen_generalized.cpp +++ b/src/scf/eigen_solver/eigen_generalized.cpp @@ -53,7 +53,7 @@ MODULE_RUN(EigenGeneralized) { Eigen::Map B_map(pB, rows, cols); // Compute - Eigen::GeneralizedEigenSolver ges(A_map, B_map); + Eigen::GeneralizedSelfAdjointEigenSolver ges(A_map, B_map); auto eigen_values = ges.eigenvalues(); auto eigen_vectors = ges.eigenvectors(); @@ -69,9 +69,9 @@ MODULE_RUN(EigenGeneralized) { typename vector_buffer::data_type vector_tensor(rows); typename matrix_buffer::data_type matrix_tensor(rows, cols); for(auto i = 0; i < rows; ++i) { - vector_tensor(i) = eigen_values(i).real(); + vector_tensor(i) = eigen_values(i); for(auto j = 0; j < cols; ++j) { - matrix_tensor(i, j) = eigen_vectors(i, j).real(); + matrix_tensor(i, j) = eigen_vectors(i, j); } } diff --git a/src/scf/matrix_builder/ao_integrals_driver.cpp b/src/scf/matrix_builder/ao_integrals_driver.cpp index b589980..bc7b296 100644 --- a/src/scf/matrix_builder/ao_integrals_driver.cpp +++ b/src/scf/matrix_builder/ao_integrals_driver.cpp @@ -31,20 +31,22 @@ 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; +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; +using rho_e_pt = simde::aos_rho_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 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 rho_e_type = simde::type::rho_e; using submods_type = pluginplay::type::submodule_map; @@ -81,6 +83,12 @@ class AODispatcher : public chemist::qm_operator::OperatorVisitor { // *m_ptensor_ = m_psubmods_->at(key).run_as(input); // } + void run(const rho_e_type& rho_e) { + chemist::braket::BraKet input(*m_pbra_, rho_e, *m_pket_); + const auto key = "Density matrix"; + *m_ptensor_ = m_psubmods_->at(key).run_as(input); + } + private: const aos* m_pbra_; const aos* m_pket_; @@ -96,6 +104,7 @@ MODULE_CTOR(AOIntegralsDriver) { add_submodule("Coulomb matrix"); add_submodule("Exchange matrix"); // add_submodule("Fock matrix"); + add_submodule("Density matrix"); } MODULE_RUN(AOIntegralsDriver) { diff --git a/src/scf/matrix_builder/density_matrix.cpp b/src/scf/matrix_builder/density_matrix.cpp new file mode 100644 index 0000000..027af3f --- /dev/null +++ b/src/scf/matrix_builder/density_matrix.cpp @@ -0,0 +1,91 @@ +/* + * 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_rho_e_aos; + +MODULE_CTOR(DensityMatrix) { + description(desc); + satisfies_property_type(); + add_input("cutoff") + .set_description( + "The cutoff for considering a state as part of the ensemble.") + .set_default(1E-16); +} + +MODULE_RUN(DensityMatrix) { + const auto&& [braket] = pt::unwrap_inputs(inputs); + const auto& bra_aos = braket.bra(); + const auto& op = braket.op(); + const auto& ket_aos = braket.ket(); + const auto& cutoff = inputs.at("cutoff").value(); + + const auto& mos = op.orbitals(); + const auto& c = mos.transform(); + const auto& weights = op.weights(); + auto n_aos = c.logical_layout().shape().as_smooth().extent(0); + + // TODO: Compare the aos and make sure they're all the same + + // Step 1: Figure out which orbitals we need to grab + using size_type = std::size_t; + std::vector participants; + for(size_type i = 0; i < weights.size(); ++i) + if(std::fabs(weights[i]) >= cutoff) participants.push_back(i); + + // For now we require that participants == [0, participants.size()) + for(size_type i = 0; i < participants.size(); ++i) + if(participants[i] != i) + throw std::runtime_error( + "Please shuffle your orbitals so that the ensemble is a " + "contiguous slice of the orbitals."); + + // Step 2: Grab the orbitals in the ensemble + using allocator_type = tensorwrapper::allocator::Eigen; + using buffer_type = typename allocator_type::eigen_buffer_type; + using eigen_tensor_type = typename buffer_type::data_type; + + const auto& c_buffer = allocator_type::rebind(c.buffer()); + const auto& c_eigen = c_buffer.value(); + + using slice_t = Eigen::array; + slice_t offsets{0, 0}; + slice_t extents{n_aos, Eigen::Index(participants.size())}; + eigen_tensor_type slice = c_eigen.slice(offsets, extents); + + // Step 3: CC_dagger + using index_pair_t = Eigen::IndexPair; + Eigen::array modes{index_pair_t(1, 1)}; + eigen_tensor_type p_eigen = slice.contract(slice, modes); + tensorwrapper::shape::Smooth p_shape(n_aos, n_aos); + tensorwrapper::layout::Physical l(p_shape); + buffer_type p_buffer(p_eigen, l); + simde::type::tensor p(p_shape, p_buffer); + + auto rv = results(); + return pt::wrap_results(rv, p); +} + +} // namespace scf::matrix_builder \ No newline at end of file diff --git a/src/scf/matrix_builder/determinant_driver.cpp b/src/scf/matrix_builder/determinant_driver.cpp new file mode 100644 index 0000000..45ac0d6 --- /dev/null +++ b/src/scf/matrix_builder/determinant_driver.cpp @@ -0,0 +1,143 @@ +/* + * 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 ao_pt = simde::aos_op_base_aos; +using fock_pt = simde::aos_f_e_aos; + +template +using pt = simde::eval_braket; + +template +using density_op_t = simde::type::rho_e; + +using ao_type = simde::type::aos; +using op_type = simde::type::op_base_type; +using braket_base = simde::type::braket; + +template +class DeterminantDispatcher : 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 submodule_map = pluginplay::type::submodule_map; + + DeterminantDispatcher(const WFType& bra, const WFType& ket, + submodule_map& submods) : + m_pbra_(&bra), m_pket_(&ket), m_psubmods_(&submods) {} + + void run(const T_e_type& T_e) { + simde::type::electron e; + simde::type::t_e_type t_e(e); + run_(t_e); + } + + void run(const V_en_type& V_en) { + simde::type::electron e; + simde::type::v_en_type v_en(e, V_en.rhs_particle().as_nuclei()); + run_(v_en); + } + + void run(const J_e_type& J_e) { + simde::type::electron e; + simde::type::j_e_type j_e(e, J_e.rhs_particle()); + run_(j_e); + } + + void run(const K_e_type& K_e) { + simde::type::electron e; + simde::type::k_e_type k_e(e, K_e.rhs_particle()); + run_(k_e); + } + + template + void run_(const OpType& o) { + const auto& aos = m_pbra_->orbitals().from_space(); + braket_base o_mn(aos, o, aos); + auto& mod = m_psubmods_->at("Two center evaluator"); + const auto& t = mod.run_as(o_mn); + m_pt = t; + } + + simde::type::tensor m_pt; + +private: + const WFType* m_pbra_; + const WFType* m_pket_; + submodule_map* m_psubmods_; +}; + +} // namespace + +MODULE_CTOR(DeterminantDriver) { + using wf_type = simde::type::rscf_wf; + description(desc); + satisfies_property_type>(); + add_submodule("Two center evaluator"); + // For now Fock matrix is special + add_submodule("Fock matrix"); +} + +MODULE_RUN(DeterminantDriver) { + using wf_type = simde::type::rscf_wf; + const auto&& [braket] = pt::unwrap_inputs(inputs); + const auto& bra = braket.bra(); + const auto& op = braket.op(); + const auto& ket = braket.ket(); + + // Step 1: Build density matrix + const auto& mos = bra.orbitals(); + const auto& aos = mos.from_space(); + const auto& occ = bra.occupations(); + + density_op_t rho_hat(mos, occ); + braket_base p_mn(aos, rho_hat, aos); + + auto& mod = submods.at("Two center evaluator"); + const auto& rho = mod.run_as(p_mn); + + // Step 2: Build the other matrix + DeterminantDispatcher visitor(bra, ket, submods); + op.visit(visitor); + const auto& t = visitor.m_pt; + + // Step 3: Contract + using allocator_type = tensorwrapper::allocator::Eigen; + using scalar_type = tensorwrapper::buffer::Eigen::data_type; + + const auto& t_eigen = allocator_type::rebind(t.buffer()).value(); + const auto& p_eigen = allocator_type::rebind(rho.buffer()).value(); + + using index_pair_t = Eigen::IndexPair; + Eigen::array modes{index_pair_t(0, 0), index_pair_t(1, 1)}; + scalar_type x = p_eigen.contract(t_eigen, modes); + + auto rv = results(); + return pt::wrap_results(rv, x()); +} + +} // namespace scf::matrix_builder diff --git a/src/scf/matrix_builder/electronic_energy.cpp b/src/scf/matrix_builder/electronic_energy.cpp new file mode 100644 index 0000000..66b0fc6 --- /dev/null +++ b/src/scf/matrix_builder/electronic_energy.cpp @@ -0,0 +1,67 @@ +/* + * 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 simde::type::electronic_hamiltonian; + +template +using pt = simde::eval_braket; + +template +using det_pt = simde::eval_braket; + +MODULE_CTOR(ElectronicEnergy) { + using wf_type = simde::type::rscf_wf; + description(desc); + satisfies_property_type>(); + add_submodule>("determinant driver"); +} + +MODULE_RUN(ElectronicEnergy) { + using wf_type = simde::type::rscf_wf; + const auto&& [H_ij] = pt::unwrap_inputs(inputs); + const auto& bra = H_ij.bra(); + const auto& H = H_ij.op(); + const auto& ket = H_ij.ket(); + + auto& mod = submods.at("determinant driver"); + double e = 0.0; + + 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); + e += ci * o; + } + + auto rv = results(); + return pt::wrap_results(rv, e); +} + +} // 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 index b1e6345..e4467df 100644 --- a/src/scf/matrix_builder/matrix_builder.hpp +++ b/src/scf/matrix_builder/matrix_builder.hpp @@ -20,12 +20,18 @@ namespace scf::matrix_builder { DECLARE_MODULE(AOIntegralsDriver); +DECLARE_MODULE(DensityMatrix); +DECLARE_MODULE(DeterminantDriver); +DECLARE_MODULE(ElectronicEnergy); 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("Density matrix builder"); + mm.add_module("Determinant driver"); + mm.add_module("Electronic energy"); mm.add_module("Fock matrix builder"); mm.add_module("Four center J builder"); mm.add_module("Four center K builder"); @@ -37,8 +43,15 @@ inline void set_defaults(pluginplay::ModuleManager& mm) { 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(ao_driver, "Density matrix", "Density matrix builder"); 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); } } // 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 34019be..bb9bf52 100644 --- a/src/scf/scf_mm.cpp +++ b/src/scf/scf_mm.cpp @@ -14,6 +14,7 @@ * limitations under the License. */ +#include "driver/driver.hpp" #include "eigen_solver/eigen_solver.hpp" #include "fock_operator/fock_operator.hpp" #include "guess/guess.hpp" @@ -25,6 +26,7 @@ namespace scf { void load_modules(pluginplay::ModuleManager& mm) { + driver::load_modules(mm); eigen_solver::load_modules(mm); fock_operator::load_modules(mm); guess::load_modules(mm); @@ -32,10 +34,13 @@ void load_modules(pluginplay::ModuleManager& mm) { update::load_modules(mm); mm.add_module("Coulomb's Law"); + + // mm.add_module("Driver"); #ifdef BUILD_TAMM_SCF mm.add_module("SCF Energy via TAMM"); #endif + driver::set_defaults(mm); guess::set_defaults(mm); matrix_builder::set_defaults(mm); update::set_defaults(mm); diff --git a/src/scf/scf_modules.hpp b/src/scf/scf_modules.hpp index e808a79..5508461 100644 --- a/src/scf/scf_modules.hpp +++ b/src/scf/scf_modules.hpp @@ -23,6 +23,8 @@ namespace scf { DECLARE_MODULE(TAMMEnergy); #endif +// DECLARE_MODULE(SCFDriver); +// DECLARE_MODULE(SCFLoop); DECLARE_MODULE(CoulombsLaw); } // namespace scf diff --git a/tests/cxx/integration_tests/coulombs_law.cpp b/tests/cxx/integration_tests/coulombs_law.cpp new file mode 100644 index 0000000..bf7525c --- /dev/null +++ b/tests/cxx/integration_tests/coulombs_law.cpp @@ -0,0 +1,34 @@ +/* + * 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::charge_charge_interaction; +using Catch::Matchers::WithinAbs; + +TEST_CASE("CoulombsLaw") { + auto mm = test_scf::load_modules(); + auto& mod = mm.at("Coulomb's Law"); + + auto h2_nuclei = test_scf::make_h2(); + // TODO: Conversions are missing in Chemist. Use those when they're in place + simde::type::charges qs; + for(const auto& nucleus : h2_nuclei) + qs.push_back(nucleus.as_point_charge()); + + auto e_nuclear = mod.run_as(qs, qs); + REQUIRE_THAT(e_nuclear, WithinAbs(0.71510297482837526, 1E-6)); +} \ No newline at end of file diff --git a/tests/cxx/integration_tests/driver/scf_driver.cpp b/tests/cxx/integration_tests/driver/scf_driver.cpp new file mode 100644 index 0000000..417f187 --- /dev/null +++ b/tests/cxx/integration_tests/driver/scf_driver.cpp @@ -0,0 +1,30 @@ +/* + * 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 Catch::Matchers::WithinAbs; + +using pt = simde::AOEnergy; + +TEST_CASE("SCFDriver") { + auto mm = test_scf::load_modules(); + auto h2 = test_scf::make_h2(); + auto aos = test_scf::h2_aos().ao_basis_set(); + + const auto e = mm.run_as("SCF Driver", aos, h2); + REQUIRE_THAT(e, WithinAbs(-1.1167592336, 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 new file mode 100644 index 0000000..5c9cbfd --- /dev/null +++ b/tests/cxx/integration_tests/driver/scf_loop.cpp @@ -0,0 +1,40 @@ +/* + * 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 Catch::Matchers::WithinAbs; + +template +using egy_pt = simde::eval_braket; + +template +using pt = simde::Optimize, WFType>; + +TEST_CASE("SCFLoop") { + using wf_type = simde::type::rscf_wf; + auto mm = test_scf::load_modules(); + auto& mod = mm.at("Loop"); + + using index_set = typename wf_type::orbital_index_set_type; + wf_type psi0(index_set{0}, test_scf::h2_cmos()); + + auto H = test_scf::h2_hamiltonian(); + + chemist::braket::BraKet H_00(psi0, H, psi0); + const auto& [e, psi] = mod.run_as>(H_00, psi0); + REQUIRE_THAT(e, WithinAbs(-1.1167592336, 1E-6)); +} \ No newline at end of file diff --git a/tests/cxx/integration_tests/integration_tests.hpp b/tests/cxx/integration_tests/integration_tests.hpp index eb9be3e..ead0faa 100644 --- a/tests/cxx/integration_tests/integration_tests.hpp +++ b/tests/cxx/integration_tests/integration_tests.hpp @@ -17,6 +17,7 @@ #pragma once #include "../test_scf.hpp" #include +#include #include #include @@ -27,6 +28,10 @@ inline auto load_modules() { pluginplay::ModuleManager mm; scf::load_modules(mm); integrals::load_modules(mm); + nux::load_modules(mm); + + mm.change_submod("SCF Driver", "Hamiltonian", + "Born-Oppenheimer approximation"); mm.change_submod("Four center J builder", "Four-center ERI", "ERI4"); mm.change_submod("Four center K builder", "Four-center ERI", "ERI4"); diff --git a/tests/cxx/integration_tests/matrix_builder/ao_integrals_driver.cpp b/tests/cxx/integration_tests/matrix_builder/ao_integrals_driver.cpp index e90765e..7cd63ce 100644 --- a/tests/cxx/integration_tests/matrix_builder/ao_integrals_driver.cpp +++ b/tests/cxx/integration_tests/matrix_builder/ao_integrals_driver.cpp @@ -91,6 +91,20 @@ TEST_CASE("AOIntegralsDriver") { compare_matrices(K, K_corr); } + SECTION("Calling density matrix") { + auto& pmod = mm.at("Density matrix builder"); + auto aos = test_scf::h2_aos(); + auto cmos = test_scf::h2_cmos(); + std::vector occs{1, 0}; + simde::type::rho_e rho_hat(cmos, occs); + chemist::braket::BraKet braket(aos, rho_hat, aos); + erased_type copy_braket(braket); + const auto& P = mod.run_as(copy_braket); + using op_pt = simde::aos_rho_e_aos; + const auto& P_corr = pmod.run_as(braket); + compare_matrices(P, P_corr); + } + // Re-enable when PluginPlay doesn't choke on loops in modules // SECTION("Calling Fock Matrix") { // auto& fmod = mm.at("Fock matrix builder"); diff --git a/tests/cxx/integration_tests/matrix_builder/determinant_driver.cpp b/tests/cxx/integration_tests/matrix_builder/determinant_driver.cpp new file mode 100644 index 0000000..3655ed8 --- /dev/null +++ b/tests/cxx/integration_tests/matrix_builder/determinant_driver.cpp @@ -0,0 +1,73 @@ +/* + * 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 Catch::Matchers::WithinAbs; + +template +using pt = simde::eval_braket; + +using simde::type::tensor; + +template +using erased_type = + chemist::braket::BraKet; + +TEST_CASE("DeterminantDriver") { + auto mm = test_scf::load_modules(); + auto mod = mm.at("Determinant driver"); + + using wf_type = simde::type::rscf_wf; + using index_set = typename wf_type::orbital_index_set_type; + + wf_type psi(index_set{0}, test_scf::h2_cmos()); + simde::type::many_electrons es{2}; + + SECTION("Calling Kinetic") { + simde::type::T_e_type T_e(es); + chemist::braket::BraKet braket(psi, T_e, psi); + erased_type copy_braket(braket); + const auto& T = mod.run_as>(copy_braket); + REQUIRE_THAT(T, WithinAbs(0.637692, 1E-6)); + } + + SECTION("Calling Electron-Nuclear Attraction") { + auto h2_nuclei = test_scf::make_h2(); + simde::type::V_en_type V_en(es, h2_nuclei); + chemist::braket::BraKet braket(psi, V_en, psi); + erased_type copy_braket(braket); + const auto& V = mod.run_as>(copy_braket); + REQUIRE_THAT(V, WithinAbs(-1.96830777255516853, 1E-6)); + } + + SECTION("Calling J") { + auto rho = test_scf::h2_density(); + simde::type::J_e_type J_e(es, rho); + chemist::braket::BraKet braket(psi, J_e, psi); + erased_type copy_braket(braket); + const auto& J = mod.run_as>(copy_braket); + REQUIRE_THAT(J, WithinAbs(0.76056339681664897, 1E-6)); + } + + SECTION("Calling K") { + auto rho = test_scf::h2_density(); + simde::type::K_e_type K_e(es, rho); + chemist::braket::BraKet braket(psi, K_e, psi); + erased_type copy_braket(braket); + const auto& K = mod.run_as>(copy_braket); + REQUIRE_THAT(K, WithinAbs(0.76056339681664897, 1E-6)); + } +} \ No newline at end of file diff --git a/tests/cxx/integration_tests/matrix_builder/electronic_energy.cpp b/tests/cxx/integration_tests/matrix_builder/electronic_energy.cpp new file mode 100644 index 0000000..28406f2 --- /dev/null +++ b/tests/cxx/integration_tests/matrix_builder/electronic_energy.cpp @@ -0,0 +1,48 @@ +/* + * 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 Catch::Matchers::WithinAbs; + +template +using pt = + simde::eval_braket; + +TEST_CASE("ElectronicEnergy") { + auto mm = test_scf::load_modules(); + auto mod = mm.at("Electronic energy"); + + using wf_type = simde::type::rscf_wf; + using index_set = typename wf_type::orbital_index_set_type; + + wf_type psi(index_set{0}, test_scf::h2_cmos()); + simde::type::many_electrons es{2}; + + simde::type::T_e_type T_e(es); + + auto h2_nuclei = test_scf::make_h2(); + simde::type::V_en_type V_en(es, h2_nuclei); + + 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); + + const auto& E_elec = mod.run_as>(braket); + REQUIRE_THAT(E_elec, WithinAbs(-1.90066758625308307, 1E-6)); +} \ No newline at end of file diff --git a/tests/cxx/integration_tests/update/diagonalization.cpp b/tests/cxx/integration_tests/update/diagonalization.cpp index 24cf377..4b46ec0 100644 --- a/tests/cxx/integration_tests/update/diagonalization.cpp +++ b/tests/cxx/integration_tests/update/diagonalization.cpp @@ -40,12 +40,15 @@ TEST_CASE("Diagaonalization") { 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); + + // Check orbital energies const auto& evals = psi.orbitals().diagonalized_matrix(); - using allocator_type = tensorwrapper::allocator::Eigen; - const auto& eval_buffer = allocator_type::rebind(evals.buffer()); + using vector_allocator = tensorwrapper::allocator::Eigen; + const auto& eval_buffer = vector_allocator::rebind(evals.buffer()); const auto tol = 1E-6; using Catch::Matchers::WithinAbs; diff --git a/tests/cxx/test_scf.hpp b/tests/cxx/test_scf.hpp index 4e5a4d6..62f8c31 100644 --- a/tests/cxx/test_scf.hpp +++ b/tests/cxx/test_scf.hpp @@ -18,6 +18,7 @@ #include #include #include +#include #include namespace test_scf { @@ -97,6 +98,13 @@ inline auto h2_mos() { return mos_type(h2_aos(), std::move(c)); } +inline auto h2_cmos() { + using cmos_type = simde::type::cmos; + using tensor_type = typename cmos_type::transform_type; + tensor_type e({-1.25330893, -0.47506974}); + return cmos_type(e, h2_aos(), h2_mos().transform()); +} + inline auto h2_density() { using density_type = simde::type::decomposable_e_density; typename density_type::value_type rho( diff --git a/tests/cxx/unit_tests/matrix_builder/density_matrix.cpp b/tests/cxx/unit_tests/matrix_builder/density_matrix.cpp new file mode 100644 index 0000000..bad620d --- /dev/null +++ b/tests/cxx/unit_tests/matrix_builder/density_matrix.cpp @@ -0,0 +1,40 @@ +/* + * 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" + +using pt = simde::aos_rho_e_aos; + +TEST_CASE("Density Matrix Builder") { + pluginplay::ModuleManager mm; + scf::load_modules(mm); + auto& mod = mm.at("Density matrix builder"); + auto aos = test_scf::h2_aos(); + auto cmos = test_scf::h2_cmos(); + std::vector occs{1, 0}; + simde::type::rho_e rho_hat(cmos, occs); + + chemist::braket::BraKet p_mn(aos, rho_hat, aos); + const auto& P = mod.run_as(p_mn); + using allocator_type = tensorwrapper::allocator::Eigen; + const auto& P_eigen = allocator_type::rebind(P.buffer()).value(); + + using Catch::Matchers::WithinAbs; + REQUIRE_THAT(P_eigen(0, 0), WithinAbs(0.31980835, 1E-6)); + REQUIRE_THAT(P_eigen(0, 1), WithinAbs(0.31980835, 1E-6)); + REQUIRE_THAT(P_eigen(1, 0), WithinAbs(0.31980835, 1E-6)); + REQUIRE_THAT(P_eigen(1, 1), WithinAbs(0.31980835, 1E-6)); +} \ No newline at end of file