diff --git a/src/scf/driver/driver.hpp b/src/scf/driver/driver.hpp index 3c585e5..c4272d6 100644 --- a/src/scf/driver/driver.hpp +++ b/src/scf/driver/driver.hpp @@ -30,7 +30,9 @@ inline void load_modules(pluginplay::ModuleManager& mm) { 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", "Diagonalizer", + "Generalized eigensolve via Eigen"); + mm.change_submod("Loop", "Fock matrix builder", "Fock matrix builder"); mm.change_submod("Loop", "One-electron Fock operator", "Restricted One-Electron Fock op"); mm.change_submod("Loop", "Fock operator", "Restricted Fock Op"); diff --git a/src/scf/driver/scf_loop.cpp b/src/scf/driver/scf_loop.cpp index 1e90dcb..6b334c5 100644 --- a/src/scf/driver/scf_loop.cpp +++ b/src/scf/driver/scf_loop.cpp @@ -20,8 +20,12 @@ namespace scf::driver { namespace { +// Comparison between convergence metrics based on floating point type struct Kernel { - explicit Kernel(parallelzone::runtime::RuntimeView rv) : m_rv(rv) {} + using rv_t = parallelzone::runtime::RuntimeView; + + explicit Kernel(rv_t rv) : m_rv(rv) {} + template auto run(const tensorwrapper::buffer::BufferBase& a, double tol) { tensorwrapper::allocator::Eigen allocator(m_rv); @@ -29,7 +33,18 @@ struct Kernel { return tensorwrapper::types::fabs(eigen_a.get_data(0)) < FloatType(tol); } - parallelzone::runtime::RuntimeView m_rv; + rv_t m_rv; +}; + +// Pull out nuclear-nuclear interaction term, if there is one. +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 = nullptr; }; const auto desc = R"( @@ -37,19 +52,20 @@ const auto desc = R"( } // namespace +using tensor_t = simde::type::tensor; +using cmos_t = simde::type::cmos; +using diis_t = tensorwrapper::diis::DIIS; +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; +using fock_matrix_pt = simde::aos_f_e_aos; +using diagonalizer_pt = simde::GeneralizedEigenSolve; +using s_pt = simde::aos_s_e_aos; using simde::type::electronic_hamiltonian; -using simde::type::hamiltonian; -using simde::type::op_base_type; - -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; -using fock_matrix_pt = simde::aos_f_e_aos; -using s_pt = simde::aos_s_e_aos; template -using egy_pt = simde::eval_braket; +using egy_pt = simde::eval_braket; template using elec_egy_pt = simde::eval_braket; @@ -60,15 +76,7 @@ using pt = simde::Optimize, WfType>; template using update_pt = simde::UpdateGuess; -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 = nullptr; -}; +using tensorwrapper::utilities::floating_point_dispatch; MODULE_CTOR(SCFLoop) { using wf_type = simde::type::rscf_wf; @@ -81,47 +89,59 @@ MODULE_CTOR(SCFLoop) { add_input("density tolerance").set_default(1.0E-6); add_input("gradient tolerance").set_default(1.0E-6); + const std::size_t diis_sample_default = 5; + add_input("DIIS").set_default(true); + add_input("DIIS max samples").set_default(diis_sample_default); + add_submodule>("Electronic energy"); add_submodule("Density matrix"); - add_submodule>("Guess update"); + add_submodule("Overlap matrix builder"); + add_submodule("Fock matrix builder"); + add_submodule("Diagonalizer"); add_submodule("One-electron Fock operator"); add_submodule("Fock operator"); - add_submodule("Fock matrix builder"); add_submodule("Charge-charge"); - add_submodule("Overlap matrix builder"); } MODULE_RUN(SCFLoop) { - using wf_type = simde::type::rscf_wf; - using density_op_type = simde::type::rho_e; + 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(); const auto max_iter = inputs.at("max iterations").value(); const auto e_tol = inputs.at("energy tolerance").value(); const auto dp_tol = inputs.at("density tolerance").value(); const auto g_tol = inputs.at("gradient tolerance").value(); - 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"); - // TODO: should be split off into orbital gradient module - auto& F_mod = submods.at("Fock matrix builder"); - auto& S_mod = submods.at("Overlap matrix builder"); - - // Step 1: Nuclear-nuclear repulsion + auto& egy_mod = submods.at("Electronic energy"); + auto& density_mod = submods.at("Density matrix"); + auto& diagonalizer_mod = submods.at("Diagonalizer"); + 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& F_mod = submods.at("Fock matrix builder"); + auto& S_mod = submods.at("Overlap matrix builder"); + + // Unwrap Braket and trial wavefunction + // 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(); + + // DIIS settings + const auto diis_on = inputs.at("DIIS").value(); + const auto diis_max_samples = + inputs.at("DIIS max samples").value(); + diis_t diis(diis_max_samples); + + // Nuclear-nuclear repulsion GrabNuclear visitor; H.visit(visitor); bool has_nn = (visitor.m_pv != nullptr); // TODO: Clean up charges class to make this easier... - simde::type::tensor e_nuclear(0.0); + tensor_t e_nuclear(0.0); if(has_nn) { const auto& V_nn = *visitor.m_pv; const auto n_lhs = V_nn.lhs_particle().as_nuclei(); @@ -143,101 +163,113 @@ MODULE_RUN(SCFLoop) { chemist::braket::BraKet s_mn(aos, simde::type::s_e_type{}, aos); const auto& S = S_mod.run_as(s_mn); - wf_type psi_old = psi0; - simde::type::tensor e_old; - - 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()); - - auto& logger = get_runtime().logger(); + // For convergence checking + wf_type psi_old; + density_t rho_old; + tensor_t e_old; + tensor_t F_old; + // Initialize loop unsigned int iter = 0; + auto& logger = get_runtime().logger(); while(iter < max_iter) { - // Step 2: Old density is used to create the corresponding 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_old = fock_mod.run_as(H, rho_old); - - // Step 3: Fock operator used to compute new wavefunction/density - const auto& psi_new = - update_mod.run_as>(f_old, psi_old); + // Step 1: Generate trial wavefunction + wf_type psi; + if(iter > 0) { + // Diagonalize the Fock matrix + const auto&& [evalues, evectors] = + diagonalizer_mod.run_as(F_old, S); + + // Construct new trial wavefunction + cmos_t cmos(evalues, aos, evectors); + psi = wf_type(psi_old.orbital_indices(), cmos); + } else { + // Use trial wavefunction provided from initial guess + psi = psi0; + } - density_op_type rho_hat_new(psi_new.orbitals(), psi_new.occupations()); - chemist::braket::BraKet P_mn_new(aos, rho_hat_new, aos); - const auto& P_new = density_mod.run_as(P_mn_new); + // Step 2: Construct electronic density + density_op_type rho_hat(psi.orbitals(), psi.occupations()); + chemist::braket::BraKet P_mn(aos, rho_hat, aos); + const auto& P = density_mod.run_as(P_mn); + density_t rho(P, psi.orbitals()); - density_t rho_new(P_new, psi_new.orbitals()); + // Step 3: Construct new Fock matrix + const auto& f_hat = fock_mod.run_as(H, rho); + chemist::braket::BraKet f_mn(aos, f_hat, aos); + auto F = F_mod.run_as(f_mn); // Step 4: New electronic energy // Step 4a: New Fock operator to new electronic Hamiltonian - // TODO: Should just be H_core + F; - const auto& F_new = Fock_mod.run_as(H, rho_new); + // TODO: Should just be H_core + F_hat; + const auto& F_hat = Fock_mod.run_as(H, rho); 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()); + for(std::size_t i = 0; i < F_hat.size(); ++i) + H_new.emplace_back(F_hat.coefficient(i), + F_hat.get_operator(i).clone()); // Step 4b: New electronic hamiltonian to new electronic energy - chemist::braket::BraKet H_00(psi_new, H_new, psi_new); - auto e_new = egy_mod.run_as>(H_00); + chemist::braket::BraKet H_00(psi, H_new, psi); + auto e = egy_mod.run_as>(H_00); auto e_msg = "SCF iteration = " + std::to_string(iter) + ":"; - e_msg += " Electronic Energy = " + e_new.to_string(); + e_msg += " Electronic Energy = " + e.to_string(); logger.log(e_msg); - bool converged = false; // Step 5: Converged? + bool converged = false; if(iter > 0) { // Change in the energy - simde::type::tensor de; - de("") = e_new("") - e_old(""); + tensor_t de; + de("") = e("") - e_old(""); // Change in the density - simde::type::tensor dp; + tensor_t dp; const auto& P_old = rho_old.value(); - dp("m,n") = rho_new.value()("m,n") - P_old("m,n"); + dp("m,n") = rho.value()("m,n") - P_old("m,n"); auto dp_norm = tensorwrapper::operations::infinity_norm(dp); // Orbital gradient: FPS-SPF - // TODO: module satisfying BraKet(aos, Commutator(F,P), aos) - const auto& f_new = fock_mod.run_as(H, rho_new); - chemist::braket::BraKet F_mn(aos, f_new, aos); - const auto& F_matrix = F_mod.run_as(F_mn); - - auto grad = commutator(F_matrix, P_new, S); - - simde::type::tensor grad_norm; + auto grad = commutator(F, P, S); + tensor_t grad_norm; grad_norm("") = grad("m,n") * grad("n,m"); - Kernel k(get_runtime()); + // Log convergence metrics + logger.log(" dE = " + de.to_string()); + logger.log(" dP = " + dp_norm.to_string()); + logger.log(" dG = " + grad_norm.to_string()); - using tensorwrapper::utilities::floating_point_dispatch; + // Check for convergence + Kernel k(get_runtime()); 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); + if(e_conv && g_conv && dp_conv) converged = true; - logger.log(" dE = " + de.to_string()); - logger.log(" dP = " + dp_norm.to_string()); - logger.log(" dG = " + grad_norm.to_string()); + // If using DIIS and not converged, extrapolate new Fock matrix + if(diis_on && !converged) { F = diis.extrapolate(F, grad); } + } else if(diis_on) { + // For DIIS, still need to the orbital gradient + auto grad = commutator(F, P, S); - if(e_conv && g_conv && dp_conv) converged = true; + // If using DIIS, extrapolate new Fock matrix + F = diis.extrapolate(F, grad); } // Step 6: Not converged so reset - e_old = e_new; - psi_old = psi_new; - rho_old = rho_new; + e_old = e; + psi_old = psi; + rho_old = rho; + F_old = F; if(converged) break; ++iter; } if(iter == max_iter) throw std::runtime_error("SCF failed to converge"); - simde::type::tensor e_total; + tensor_t e_total; // e_nuclear is a double. This hack converts it to udouble (if needed) tensorwrapper::allocator::Eigen dalloc(get_runtime()); @@ -245,7 +277,7 @@ MODULE_RUN(SCFLoop) { tensorwrapper::allocator::Eigen ualloc(get_runtime()); if(ualloc.can_rebind(e_old.buffer())) { - simde::type::tensor temp(e_old); + tensor_t temp(e_old); auto val = dalloc.rebind(e_nuclear.buffer()).get_data(0); ualloc.rebind(temp.buffer()).set_data(0, val); e_nuclear = temp; diff --git a/tests/cxx/integration_tests/driver/scf_loop.cpp b/tests/cxx/integration_tests/driver/scf_loop.cpp index f786021..a38d451 100644 --- a/tests/cxx/integration_tests/driver/scf_loop.cpp +++ b/tests/cxx/integration_tests/driver/scf_loop.cpp @@ -41,25 +41,46 @@ TEMPLATE_LIST_TEST_CASE("SCFLoop", "", test_scf::float_types) { 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.template run_as>(H_00, psi0); - pcorr->set_elem({}, -1.1167592336); - tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); - REQUIRE(approximately_equal(corr, e, 1E-6)); + SECTION("Default") { + mod.change_input("DIIS", true); + + const auto& [e, psi] = mod.template run_as>(H_00, psi0); + pcorr->set_elem({}, -1.1167592336); + tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); + REQUIRE(approximately_equal(corr, e, 1E-6)); + } + + SECTION("With DIIS") { + const auto& [e, psi] = mod.template run_as>(H_00, psi0); + pcorr->set_elem({}, -1.1167592336); + tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); + REQUIRE(approximately_equal(corr, e, 1E-6)); + } } SECTION("He") { wf_type psi0(index_set{0}, test_scf::he_cmos()); auto H = test_scf::he_hamiltonian(); - chemist::braket::BraKet H_00(psi0, H, psi0); - const auto& [e, psi] = mod.template run_as>(H_00, psi0); - pcorr->set_elem({}, -2.807783957539); - tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); - REQUIRE(approximately_equal(corr, e, 1E-6)); + SECTION("Default") { + mod.change_input("DIIS", true); + const auto& [e, psi] = mod.template run_as>(H_00, psi0); + + pcorr->set_elem({}, -2.807783957539); + tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); + REQUIRE(approximately_equal(corr, e, 1E-6)); + } + + SECTION("With DIIS") { + const auto& [e, psi] = mod.template run_as>(H_00, psi0); + + pcorr->set_elem({}, -2.807783957539); + tensorwrapper::Tensor corr(shape_corr, std::move(pcorr)); + REQUIRE(approximately_equal(corr, e, 1E-6)); + } } } \ No newline at end of file diff --git a/tests/python/integration_tests/driver/test_scf_driver.py b/tests/python/integration_tests/driver/test_scf_driver.py index dbb1077..f7781e6 100644 --- a/tests/python/integration_tests/driver/test_scf_driver.py +++ b/tests/python/integration_tests/driver/test_scf_driver.py @@ -25,24 +25,21 @@ class TestSCFDriver(unittest.TestCase): def test_scf_driver(self): - # Property Types - mol_from_str = simde.MoleculeFromString() - mol_basis_set = simde.MolecularBasisSet() - ao_energy = simde.AOEnergy() - - # Inputs - water = self.mm.at("NWX Molecules").run_as(mol_from_str, "water") - aos = self.mm.at('STO-3G').run_as(mol_basis_set, water) - sys = chemist.ChemicalSystem(water) + self.mm.change_input("Loop", "DIIS", False) + egy = self.mm.run_as(self.ao_energy, "SCF Driver", self.aos, self.sys) + self.assertAlmostEqual(np.array(egy), -74.94208027122616, places=6) - # Run module - egy = self.mm.run_as(ao_energy, "SCF Driver", aos, sys) + def test_scf_driver_diis(self): + egy = self.mm.run_as(self.ao_energy, "SCF Driver", self.aos, self.sys) self.assertAlmostEqual(np.array(egy), -74.94208027122616, places=6) def setUp(self): + # Setup Module Manager self.mm = pp.ModuleManager(pz.runtime.RuntimeView()) nux.load_modules(self.mm) nwx.load_modules(self.mm) + + # Set Submods self.mm.change_submod("SCF Driver", "Hamiltonian", "Born-Oppenheimer Approximation") self.mm.change_submod("SCF integral driver", "Fundamental matrices", @@ -50,3 +47,13 @@ def setUp(self): self.mm.change_submod("Diagonalization Fock update", "Overlap matrix builder", "Overlap") self.mm.change_submod("Loop", "Overlap matrix builder", "Overlap") + + # Property Types + self.mol = simde.MoleculeFromString() + self.basis_set = simde.MolecularBasisSet() + self.ao_energy = simde.AOEnergy() + + # Inputs + self.water = self.mm.at("NWX Molecules").run_as(self.mol, "water") + self.aos = self.mm.at('STO-3G').run_as(self.basis_set, self.water) + self.sys = chemist.ChemicalSystem(self.water) diff --git a/uncertainty/eigen_generalized.py b/uncertainty/eigen_generalized.py index fc5cb39..2f02b49 100644 --- a/uncertainty/eigen_generalized.py +++ b/uncertainty/eigen_generalized.py @@ -1,10 +1,26 @@ +# 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. + import numpy as np from pluginplay import ModuleManager from nwchemex import load_modules from simde import GeneralizedEigenSolve from tensorwrapper import Tensor + class GeneralizedEigenSolverTester: + def __init__(self): self.mm = ModuleManager() load_modules(self.mm) @@ -15,25 +31,30 @@ def solve_gen_eigenproblem(self, A, B, verify=True): A = np.asarray(A, dtype=np.float64) B = np.asarray(B, dtype=np.float64) A_tensor, B_tensor = Tensor(A), Tensor(B) - + # Solve - λ, V = map(np.array, self.solver.run_as(GeneralizedEigenSolve(), A_tensor, B_tensor)) - + λ, V = map( + np.array, + self.solver.run_as(GeneralizedEigenSolve(), A_tensor, B_tensor)) + if verify: max_residual = self._verify_results(A, B, λ, V) print(f"Maximum residual norm: {max_residual:.2e}") - + return λ, V - + def _verify_results(self, A, B, λ, V, rtol=1e-8): """Internal verification method""" max_residual = 0 for i, ev in enumerate(λ): - res = A @ V[:,i] - ev * (B @ V[:,i]) + res = A @ V[:, i] - ev * (B @ V[:, i]) max_residual = max(max_residual, np.linalg.norm(res)) return max_residual - def analyze_uncertainty(self, num_trials=100, matrix_size=5, noise_level=1e-10): + def analyze_uncertainty(self, + num_trials=100, + matrix_size=5, + noise_level=1e-10): """ Analyze numerical uncertainty in generalized eigensolve by testing pertubed matricies @@ -47,9 +68,9 @@ def analyze_uncertainty(self, num_trials=100, matrix_size=5, noise_level=1e-10): eigenvector_errors = [] for _ in range(num_trials): - #Generate random symmetric matrices + #Generate random symmetric matrices A = np.random.rand(matrix_size, matrix_size) - A = (A + A.T) / 2 + A = (A + A.T) / 2 B = np.eye(matrix_size) A_perturbed = A + noise_level * np.random.randn(*A.shape) @@ -59,12 +80,15 @@ def analyze_uncertainty(self, num_trials=100, matrix_size=5, noise_level=1e-10): λ_p, V_p = self.solve_gen_eigenproblem(A_perturbed, B, False) eigenvalue_errors.append(np.abs(λ - λ_p)) - eigenvector_errors.append([np.linalg.norm(V[:,i] - V_p[:,i]) - for i in range(matrix_size)]) + eigenvector_errors.append([ + np.linalg.norm(V[:, i] - V_p[:, i]) for i in range(matrix_size) + ]) #results print("\n=== Uncertainty Analysis ===") - print(f"Tested {num_trials} random {matrix_size}*{matrix_size} matricies") + print( + f"Tested {num_trials} random {matrix_size}*{matrix_size} matricies" + ) print(f"Pertubation level: {noise_level:.1e}") print("\n Eigenvalue Error Statistics:") @@ -72,8 +96,12 @@ def analyze_uncertainty(self, num_trials=100, matrix_size=5, noise_level=1e-10): print(f"Max absolute error: {np.max(eigenvalue_errors):.2e}") print("\n Eigenvector Error Statistics:") - print(f"Mean angular differences: {np.mean(eigenvector_errors):.2e} radians") - print(f"Max angular difference: {np.max(eigenvector_errors):.2e} radians") + print( + f"Mean angular differences: {np.mean(eigenvector_errors):.2e} radians" + ) + print( + f"Max angular difference: {np.max(eigenvector_errors):.2e} radians" + ) def run_standard_example(self): """Run a fixed example""" @@ -82,18 +110,17 @@ def run_standard_example(self): λ, V = self.solve_gen_eigenproblem(A, B) print("Eigenvalues:\n", λ) + def main(): tester = GeneralizedEigenSolverTester() - tester.analyze_uncertainty( - num_trials=100, - matrix_size=5, - noise_level=1e-10 - ) + tester.analyze_uncertainty(num_trials=100, + matrix_size=5, + noise_level=1e-10) print("\n=== Standard Example ===") tester.run_standard_example() + if __name__ == "__main__": main() -