diff --git a/.gitignore b/.gitignore index 6478ba2..cc414ab 100644 --- a/.gitignore +++ b/.gitignore @@ -20,6 +20,7 @@ # These are directories used by IDEs for storing settings .idea/ .vscode/ +.cache/ # These are common Python virtual enviornment directory names venv/ diff --git a/include/scf/driver/convergence.hpp b/include/scf/driver/convergence.hpp new file mode 100644 index 0000000..f130e2e --- /dev/null +++ b/include/scf/driver/convergence.hpp @@ -0,0 +1,54 @@ +/* + * 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 { + +template +DECLARE_TEMPLATED_PROPERTY_TYPE(ConvergenceProp, KernelType); + +template +TEMPLATED_PROPERTY_TYPE_INPUTS(ConvergenceProp, KernelType) { + auto rv = pluginplay::declare_input() + .add_field("New Energy") + .template add_field("Old Energy") + .template add_field>("New Rho") + .template add_field>("Old Rho") + .template add_field("P_new") + .template add_field("S") + .template add_field("Fock Operator") + .template add_field("Wave Function") + .template add_field("K") + .template add_field("Energy Tolerance") + .template add_field("Density Tolerance") + .template add_field("Gradient Tolerance"); + return rv; +} + +template +TEMPLATED_PROPERTY_TYPE_RESULTS(ConvergenceProp, KernelType) { + auto rv = pluginplay::declare_result().add_field("Convergence Status"); + rv.at("Convergence Status") + .set_description("The molecule corresponding to the input string"); + return rv; +} + +} // namespace simde + diff --git a/src/scf/driver/convergence.cpp b/src/scf/driver/convergence.cpp new file mode 100644 index 0000000..1276805 --- /dev/null +++ b/src/scf/driver/convergence.cpp @@ -0,0 +1,87 @@ +/* + * 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" +#include + + +namespace scf::driver { + +template +TEMPLATED_MODULE_CTOR(ConvergenceMod, KernelType) { + satisfies_property_type>(); + + add_input("energy tolerance").set_default(1.0E-6); + add_input("density tolerance").set_default(1.0E-6); + add_input("gradient tolerance").set_default(1.0E-6); + add_submodule("Fock matrix builder"); +} + +template +TEMPLATED_MODULE_RUN(ConvergenceMod, KernelType) { + const auto&& [e_new, e_old, rho_new, rho_old, P_new, S, f_new, aos, k, e_tol, dp_tol, g_tol] = ConvergenceProp::unwrap_inputs(inputs); + + bool converged = false; + + auto& F_mod = submods.at("Fock matrix builder"); + + simde::type::tensor de; + de("") = e_new("") - e_old(""); + + // Change in the density + simde::type::tensor dp; + dp("m,n") = rho_new.value()("m,n") - rho_old.value()("m,n"); + auto dp_norm = tensorwrapper::operations::infinity_norm(dp); + + // Orbital gradient: FPS-SPF + // TODO: module satisfying BraKet(aos, Commutator(F,P), aos) + chemist::braket::BraKet F_mn(aos, f_new, aos); + const auto& F_matrix = F_mod.run_as(F_mn); + simde::type::tensor FPS; + FPS("m,l") = F_matrix("m,n") * P_new("n,l"); + FPS("m,l") = FPS("m,n") * S("n,l"); + + simde::type::tensor SPF; + SPF("m,l") = P_new("m,n") * F_matrix("n,l"); + SPF("m,l") = S("m,n") * SPF("n,l"); + + simde::type::tensor grad; + simde::type::tensor grad_norm; + grad("m,n") = FPS("m,n") - SPF("m,n"); + grad_norm("") = grad("m,n") * grad("n,m"); + + using tensorwrapper::utilities::floating_point_dispatch; + auto e_conv = floating_point_dispatch(k, de.buffer(), e_tol); + auto g_conv = floating_point_dispatch(k, grad_norm.buffer(), g_tol); + auto dp_conv = floating_point_dispatch(k, dp_norm.buffer(), dp_tol); + + // logger.log(" dE = " + de.to_string()); + // logger.log(" dP = " + dp_norm.to_string()); + // logger.log(" dG = " + grad_norm.to_string()); + + if(e_conv && g_conv && dp_conv) converged = true; + + auto rv = results(); + return scf::ConvergenceProp::wrap_inputs(rv, converged); +} +} + // // Step 6: Not converged so reset + // e_old = e_new; + // psi_old = psi_new; + // rho_old = rho_new; + // if(converged) break; + // ++iter; + // } diff --git a/src/scf/driver/driver.hpp b/src/scf/driver/driver.hpp index e5039e6..13b72d3 100644 --- a/src/scf/driver/driver.hpp +++ b/src/scf/driver/driver.hpp @@ -16,12 +16,16 @@ #pragma once #include +#include namespace scf::driver { DECLARE_MODULE(SCFDriver); DECLARE_MODULE(SCFLoop); +template +DECLARE_MODULE(ConvergenceMod); + inline void load_modules(pluginplay::ModuleManager& mm) { mm.add_module("SCF Driver"); mm.add_module("Loop"); @@ -41,4 +45,4 @@ inline void set_defaults(pluginplay::ModuleManager& mm) { mm.change_submod("SCF Driver", "Optimizer", "Loop"); } -} // namespace scf::driver \ No newline at end of file +} // namespace scf::driver diff --git a/src/scf/driver/scf_loop.cpp b/src/scf/driver/scf_loop.cpp index 39d6475..7eca8cf 100644 --- a/src/scf/driver/scf_loop.cpp +++ b/src/scf/driver/scf_loop.cpp @@ -92,6 +92,7 @@ MODULE_CTOR(SCFLoop) { add_submodule("Fock matrix builder"); add_submodule("Charge-charge"); add_submodule("Overlap matrix builder"); + add_submodule>("Converged Mod"); } MODULE_RUN(SCFLoop) { @@ -110,6 +111,7 @@ MODULE_RUN(SCFLoop) { auto& fock_mod = submods.at("One-electron Fock operator"); auto& Fock_mod = submods.at("Fock operator"); auto& V_nn_mod = submods.at("Charge-charge"); + auto& conv_mod = submods.at("Converged Mod"); // TODO: should be split off into orbital gradient module auto& F_mod = submods.at("Fock matrix builder"); @@ -191,53 +193,54 @@ MODULE_RUN(SCFLoop) { logger.log(e_msg); bool converged = false; - // Step 5: Converged? + // // Step 5: Converged? if(iter > 0) { // Change in the energy - simde::type::tensor de; - de("") = e_new("") - e_old(""); - - // Change in the density - simde::type::tensor dp; - dp("m,n") = rho_new.value()("m,n") - rho_old.value()("m,n"); - auto dp_norm = tensorwrapper::operations::infinity_norm(dp); - - // Orbital gradient: FPS-SPF - // TODO: module satisfying BraKet(aos, Commutator(F,P), aos) - chemist::braket::BraKet F_mn(aos, f_new, aos); - const auto& F_matrix = F_mod.run_as(F_mn); - simde::type::tensor FPS; - FPS("m,l") = F_matrix("m,n") * P_new("n,l"); - FPS("m,l") = FPS("m,n") * S("n,l"); - - simde::type::tensor SPF; - SPF("m,l") = P_new("m,n") * F_matrix("n,l"); - SPF("m,l") = S("m,n") * SPF("n,l"); - - simde::type::tensor grad; - simde::type::tensor grad_norm; - grad("m,n") = FPS("m,n") - SPF("m,n"); - grad_norm("") = grad("m,n") * grad("n,m"); + // simde::type::tensor de; + // de("") = e_new("") - e_old(""); + + // // Change in the density + // simde::type::tensor dp; + // dp("m,n") = rho_new.value()("m,n") - rho_old.value()("m,n"); + // auto dp_norm = tensorwrapper::operations::infinity_norm(dp); + + // // Orbital gradient: FPS-SPF + // // TODO: module satisfying BraKet(aos, Commutator(F,P), aos) + // chemist::braket::BraKet F_mn(aos, f_new, aos); + // const auto& F_matrix = F_mod.run_as(F_mn); + // simde::type::tensor FPS; + // FPS("m,l") = F_matrix("m,n") * P_new("n,l"); + // FPS("m,l") = FPS("m,n") * S("n,l"); + + // simde::type::tensor SPF; + // SPF("m,l") = P_new("m,n") * F_matrix("n,l"); + // SPF("m,l") = S("m,n") * SPF("n,l"); + + // simde::type::tensor grad; + // simde::type::tensor grad_norm; + // grad("m,n") = FPS("m,n") - SPF("m,n"); + // grad_norm("") = grad("m,n") * grad("n,m"); Kernel k(get_runtime()); - using tensorwrapper::utilities::floating_point_dispatch; - auto e_conv = floating_point_dispatch(k, de.buffer(), e_tol); - auto g_conv = floating_point_dispatch(k, grad_norm.buffer(), g_tol); - auto dp_conv = floating_point_dispatch(k, dp_norm.buffer(), dp_tol); + // using tensorwrapper::utilities::floating_point_dispatch; + // auto e_conv = floating_point_dispatch(k, de.buffer(), e_tol); + // auto g_conv = floating_point_dispatch(k, grad_norm.buffer(), g_tol); + // auto dp_conv = floating_point_dispatch(k, dp_norm.buffer(), dp_tol); - logger.log(" dE = " + de.to_string()); - logger.log(" dP = " + dp_norm.to_string()); - logger.log(" dG = " + grad_norm.to_string()); + // logger.log(" dE = " + de.to_string()); + // logger.log(" dP = " + dp_norm.to_string()); + // logger.log(" dG = " + grad_norm.to_string()); - if(e_conv && g_conv && dp_conv) converged = true; + // // if(e_conv && g_conv && dp_conv) converged = true; + auto converged = conv_mod.run_as>(e_new, e_old, rho_new, rho_old, P_new, S, f_new, aos, k, e_tol, dp_tol, g_tol); } + if (converged) break; // Step 6: Not converged so reset e_old = e_new; psi_old = psi_new; rho_old = rho_new; - if(converged) break; ++iter; } if(iter == max_iter) throw std::runtime_error("SCF failed to converge"); @@ -261,4 +264,4 @@ MODULE_RUN(SCFLoop) { return pt::wrap_results(rv, e_total, psi_old); } -} // namespace scf::driver \ No newline at end of file +} // namespace scf::driver