From 974186a5222b57325b3dfe7c52e0eb3bfa49c94b Mon Sep 17 00:00:00 2001 From: Jan Unsleber Date: Mon, 1 Jun 2026 18:29:03 +0200 Subject: [PATCH 01/13] Add nuclear derivative calculators --- cpp/include/qdk/chemistry.hpp | 3 + .../algorithms/nuclear_derivative.hpp | 192 ++++++++++ .../qdk/chemistry/data/nuclear_gradients.hpp | 142 ++++++++ .../qdk/chemistry/data/nuclear_hessian.hpp | 137 ++++++++ .../qdk/chemistry/algorithms/CMakeLists.txt | 1 + .../algorithms/algorithm_defaults.cpp | 2 + .../chemistry/algorithms/microsoft/scf.cpp | 41 ++- .../chemistry/algorithms/microsoft/scf.hpp | 26 ++ .../algorithms/nuclear_derivative.cpp | 328 ++++++++++++++++++ cpp/src/qdk/chemistry/data/CMakeLists.txt | 2 + .../qdk/chemistry/data/nuclear_gradients.cpp | 222 ++++++++++++ .../qdk/chemistry/data/nuclear_hessian.cpp | 212 +++++++++++ cpp/tests/test_nuclear_derivative.cpp | 127 +++++++ python/CMakeLists.txt | 3 + .../algorithms/nuclear_derivative.cpp | 118 +++++++ .../src/pybind11/data/nuclear_gradients.cpp | 125 +++++++ python/src/pybind11/data/nuclear_hessian.cpp | 120 +++++++ python/src/pybind11/module.cpp | 6 + .../src/qdk_chemistry/algorithms/__init__.py | 8 + .../algorithms/nuclear_derivative.py | 18 + .../src/qdk_chemistry/algorithms/registry.py | 2 + python/src/qdk_chemistry/data/__init__.py | 6 + python/tests/test_nuclear_derivative.py | 63 ++++ 23 files changed, 1899 insertions(+), 5 deletions(-) create mode 100644 cpp/include/qdk/chemistry/algorithms/nuclear_derivative.hpp create mode 100644 cpp/include/qdk/chemistry/data/nuclear_gradients.hpp create mode 100644 cpp/include/qdk/chemistry/data/nuclear_hessian.hpp create mode 100644 cpp/src/qdk/chemistry/algorithms/nuclear_derivative.cpp create mode 100644 cpp/src/qdk/chemistry/data/nuclear_gradients.cpp create mode 100644 cpp/src/qdk/chemistry/data/nuclear_hessian.cpp create mode 100644 cpp/tests/test_nuclear_derivative.cpp create mode 100644 python/src/pybind11/algorithms/nuclear_derivative.cpp create mode 100644 python/src/pybind11/data/nuclear_gradients.cpp create mode 100644 python/src/pybind11/data/nuclear_hessian.cpp create mode 100644 python/src/qdk_chemistry/algorithms/nuclear_derivative.py create mode 100644 python/tests/test_nuclear_derivative.py diff --git a/cpp/include/qdk/chemistry.hpp b/cpp/include/qdk/chemistry.hpp index c253148c7..9392a6a98 100644 --- a/cpp/include/qdk/chemistry.hpp +++ b/cpp/include/qdk/chemistry.hpp @@ -8,11 +8,14 @@ #include #include #include +#include #include #include #include #include #include +#include +#include #include #include #include diff --git a/cpp/include/qdk/chemistry/algorithms/nuclear_derivative.hpp b/cpp/include/qdk/chemistry/algorithms/nuclear_derivative.hpp new file mode 100644 index 000000000..f13bf8b1e --- /dev/null +++ b/cpp/include/qdk/chemistry/algorithms/nuclear_derivative.hpp @@ -0,0 +1,192 @@ +// Copyright (c) Microsoft Corporation. All rights reserved. +// Licensed under the MIT License. See LICENSE.txt in the project root for +// license information. + +#pragma once +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace qdk::chemistry::algorithms { + +/** + * @brief Initial basis, orbitals, or wavefunction seed for a derivative run. + */ +using NuclearDerivativeSeedType = + std::variant, + std::shared_ptr, + std::shared_ptr, std::string>; + +/** + * @brief Energy, gradients, optional Hessian, and optional wavefunction. + */ +using NuclearDerivativeResult = + std::tuple, + std::optional>, + std::optional>>; + +/** + * @class NuclearDerivativeSettings + * @brief Settings for nuclear derivative calculations. + */ +class NuclearDerivativeSettings : public data::Settings { + public: + /** + * @brief Construct settings with finite-difference defaults. + */ + NuclearDerivativeSettings() { + set_default("energy_calculator", data::AlgorithmRef("scf_solver", "qdk")); + set_default("orbital_solver", data::AlgorithmRef("scf_solver", "qdk")); + set_default("hamiltonian_constructor", + data::AlgorithmRef("hamiltonian_constructor", "qdk")); + set_default("compute_hessian", false); + set_default("finite_difference_step", 1.0e-3, + "Nuclear displacement step in Bohr", + data::BoundConstraint{1.0e-8, 1.0}); + set_default("symmetrize_hessian", true); + set_default( + "n_active_alpha_electrons", static_cast(0), + "Active alpha electrons for multi-configuration energy paths", + data::BoundConstraint{0, std::numeric_limits::max()}); + set_default( + "n_active_beta_electrons", static_cast(0), + "Active beta electrons for multi-configuration energy paths", + data::BoundConstraint{0, std::numeric_limits::max()}); + } +}; + +/** + * @class NuclearDerivativeCalculator + * @brief Base class for nuclear derivative algorithms. + * + * Implementations compute a total energy and nuclear gradients for a molecular + * structure. Hessians are returned when requested by settings, and a + * wavefunction is returned when the selected energy path naturally produces + * one. + */ +class NuclearDerivativeCalculator + : public Algorithm, int, int, + NuclearDerivativeSeedType> { + public: + /** + * @brief Construct a calculator with nuclear derivative settings. + */ + NuclearDerivativeCalculator() { + _settings = std::make_unique(); + } + virtual ~NuclearDerivativeCalculator() = default; + + using Algorithm::run; + + virtual std::string name() const = 0; + + /** + * @brief Return the factory type name for nuclear derivative calculators. + */ + std::string type_name() const final { + return "nuclear_derivative_calculator"; + } + + protected: + /** + * @brief Implementation hook for derived nuclear derivative calculators. + */ + virtual NuclearDerivativeResult _run_impl( + std::shared_ptr structure, int charge, + int spin_multiplicity, NuclearDerivativeSeedType seed) const = 0; +}; + +/** + * @brief Factory for nuclear derivative calculator implementations. + */ +struct NuclearDerivativeCalculatorFactory + : public AlgorithmFactory { + /** + * @brief Return the algorithm type name managed by this factory. + */ + static std::string algorithm_type_name() { + return "nuclear_derivative_calculator"; + } + + /** + * @brief Register built-in nuclear derivative calculator implementations. + */ + static void register_default_instances(); + + /** + * @brief Return the default nuclear derivative implementation name. + */ + static std::string default_algorithm_name() { return "finite_difference"; } +}; + +/** + * @class FiniteDifferenceNuclearDerivativeCalculator + * @brief Numeric nuclear derivative calculator using central finite + * differences. + */ +class FiniteDifferenceNuclearDerivativeCalculator + : public NuclearDerivativeCalculator { + public: + /** + * @brief Return the implementation name. + */ + std::string name() const final { return "finite_difference"; } + + /** + * @brief Return accepted factory aliases for this implementation. + */ + std::vector aliases() const final { + return {"finite_difference", "numeric"}; + } + + protected: + /** + * @brief Compute finite-difference nuclear derivatives. + */ + NuclearDerivativeResult _run_impl( + std::shared_ptr structure, int charge, + int spin_multiplicity, + NuclearDerivativeSeedType seed_or_basis) const override; +}; + +/** + * @class QdkNuclearDerivativeCalculator + * @brief Internal QDK derivative calculator using analytic SCF gradients. + */ +class QdkNuclearDerivativeCalculator : public NuclearDerivativeCalculator { + public: + /** + * @brief Return the implementation name. + */ + std::string name() const final { return "qdk"; } + + /** + * @brief Return accepted factory aliases for this implementation. + */ + std::vector aliases() const final { + return {"qdk", "analytical_gradient"}; + } + + protected: + /** + * @brief Compute analytic nuclear gradients with the internal SCF backend. + */ + NuclearDerivativeResult _run_impl( + std::shared_ptr structure, int charge, + int spin_multiplicity, + NuclearDerivativeSeedType seed_or_basis) const override; +}; + +} // namespace qdk::chemistry::algorithms \ No newline at end of file diff --git a/cpp/include/qdk/chemistry/data/nuclear_gradients.hpp b/cpp/include/qdk/chemistry/data/nuclear_gradients.hpp new file mode 100644 index 000000000..aa635dae9 --- /dev/null +++ b/cpp/include/qdk/chemistry/data/nuclear_gradients.hpp @@ -0,0 +1,142 @@ +// Copyright (c) Microsoft Corporation. All rights reserved. +// Licensed under the MIT License. See LICENSE.txt in the project root for +// license information. + +#pragma once +#include +#include +#include +#include +#include +#include +#include + +namespace qdk::chemistry::data { + +/** + * @class NuclearGradients + * @brief Nuclear energy gradients for a specific molecular structure. + * + * Stores the Cartesian derivative of the total energy with respect to nuclear + * coordinates. Values are stored atom-major as x, y, z components for each + * atom and are associated with the structure for which they were computed. + */ +class NuclearGradients : public DataClass, + public std::enable_shared_from_this { + public: + /** + * @brief Construct nuclear gradients for a structure. + * + * @param structure Molecular structure used to compute the gradients. + * @param values Atom-major gradient vector with length 3 * number of atoms. + * @throws std::invalid_argument If the structure is null or the vector size + * does not match the structure. + */ + NuclearGradients(std::shared_ptr structure, + const Eigen::VectorXd& values); + + /** + * @brief Get the molecular structure associated with these gradients. + */ + const std::shared_ptr get_structure() const; + + /** + * @brief Get the atom-major gradient vector in Hartree/Bohr. + */ + const Eigen::VectorXd& get_values() const { return values_; } + + /** + * @brief Return gradients as a num_atoms by 3 matrix. + */ + Eigen::MatrixXd as_matrix() const; + + /** + * @brief Get the number of atoms in the associated structure. + */ + size_t get_num_atoms() const; + + /** + * @brief Return the serialized data type name. + */ + std::string get_data_type_name() const override { + return DATACLASS_TO_SNAKE_CASE(NuclearGradients); + } + + /** + * @brief Return a short summary of the gradient data. + */ + std::string get_summary() const override; + + /** + * @brief Save gradients to a JSON or HDF5 file. + */ + void to_file(const std::string& filename, + const std::string& type) const override; + + /** + * @brief Convert gradients to JSON. + */ + nlohmann::json to_json() const override; + + /** + * @brief Save gradients to a JSON file. + */ + void to_json_file(const std::string& filename) const override; + + /** + * @brief Save gradients to an HDF5 group. + */ + void to_hdf5(H5::Group& group) const override; + + /** + * @brief Save gradients to an HDF5 file. + */ + void to_hdf5_file(const std::string& filename) const override; + + /** + * @brief Load gradients from a JSON or HDF5 file. + */ + static std::shared_ptr from_file( + const std::string& filename, const std::string& type); + + /** + * @brief Load gradients from a JSON file. + */ + static std::shared_ptr from_json_file( + const std::string& filename); + + /** + * @brief Load gradients from JSON. + */ + static std::shared_ptr from_json(const nlohmann::json& j); + + /** + * @brief Load gradients from an HDF5 file. + */ + static std::shared_ptr from_hdf5_file( + const std::string& filename); + + /** + * @brief Load gradients from an HDF5 group. + */ + static std::shared_ptr from_hdf5(H5::Group& group); + + private: + static constexpr const char* SERIALIZATION_VERSION = "0.1.0"; + + bool _is_valid() const; + void _to_json_file(const std::string& filename) const; + void _to_hdf5_file(const std::string& filename) const; + static std::shared_ptr _from_json_file( + const std::string& filename); + static std::shared_ptr _from_hdf5_file( + const std::string& filename); + + std::shared_ptr structure_; + Eigen::VectorXd values_; +}; + +static_assert(DataClassCompliant, + "NuclearGradients must implement the DataClass interface"); + +} // namespace qdk::chemistry::data \ No newline at end of file diff --git a/cpp/include/qdk/chemistry/data/nuclear_hessian.hpp b/cpp/include/qdk/chemistry/data/nuclear_hessian.hpp new file mode 100644 index 000000000..7939eea09 --- /dev/null +++ b/cpp/include/qdk/chemistry/data/nuclear_hessian.hpp @@ -0,0 +1,137 @@ +// Copyright (c) Microsoft Corporation. All rights reserved. +// Licensed under the MIT License. See LICENSE.txt in the project root for +// license information. + +#pragma once +#include +#include +#include +#include +#include +#include +#include + +namespace qdk::chemistry::data { + +/** + * @class NuclearHessian + * @brief Nuclear second derivatives for a specific molecular structure. + * + * Stores the Cartesian Hessian of the total energy with respect to nuclear + * coordinates. The matrix is ordered atom-major by x, y, z components and is + * associated with the structure for which it was computed. + */ +class NuclearHessian : public DataClass, + public std::enable_shared_from_this { + public: + /** + * @brief Construct a nuclear Hessian for a structure. + * + * @param structure Molecular structure used to compute the Hessian. + * @param matrix Square 3N by 3N Hessian matrix in Hartree/Bohr^2. + * @throws std::invalid_argument If the structure is null or the matrix shape + * does not match the structure. + */ + NuclearHessian(std::shared_ptr structure, + const Eigen::MatrixXd& matrix); + + /** + * @brief Get the molecular structure associated with this Hessian. + */ + const std::shared_ptr get_structure() const; + + /** + * @brief Get the Hessian matrix in Hartree/Bohr^2. + */ + const Eigen::MatrixXd& get_matrix() const { return matrix_; } + + /** + * @brief Get the number of atoms in the associated structure. + */ + size_t get_num_atoms() const; + + /** + * @brief Return the serialized data type name. + */ + std::string get_data_type_name() const override { + return DATACLASS_TO_SNAKE_CASE(NuclearHessian); + } + + /** + * @brief Return a short summary of the Hessian data. + */ + std::string get_summary() const override; + + /** + * @brief Save the Hessian to a JSON or HDF5 file. + */ + void to_file(const std::string& filename, + const std::string& type) const override; + + /** + * @brief Convert the Hessian to JSON. + */ + nlohmann::json to_json() const override; + + /** + * @brief Save the Hessian to a JSON file. + */ + void to_json_file(const std::string& filename) const override; + + /** + * @brief Save the Hessian to an HDF5 group. + */ + void to_hdf5(H5::Group& group) const override; + + /** + * @brief Save the Hessian to an HDF5 file. + */ + void to_hdf5_file(const std::string& filename) const override; + + /** + * @brief Load the Hessian from a JSON or HDF5 file. + */ + static std::shared_ptr from_file(const std::string& filename, + const std::string& type); + + /** + * @brief Load the Hessian from a JSON file. + */ + static std::shared_ptr from_json_file( + const std::string& filename); + + /** + * @brief Load the Hessian from JSON. + */ + static std::shared_ptr from_json(const nlohmann::json& j); + + /** + * @brief Load the Hessian from an HDF5 file. + */ + static std::shared_ptr from_hdf5_file( + const std::string& filename); + + /** + * @brief Load the Hessian from an HDF5 group. + */ + static std::shared_ptr from_hdf5(H5::Group& group); + + private: + static constexpr const char* SERIALIZATION_VERSION = "0.1.0"; + + bool _is_valid() const; + void _to_json_file(const std::string& filename) const; + void _to_hdf5_file(const std::string& filename) const; + static std::shared_ptr _from_json_file( + const std::string& filename); + static std::shared_ptr _from_hdf5_file( + const std::string& filename); + + std::shared_ptr structure_; + Eigen::MatrixXd matrix_; +}; + +static_assert(DataClassCompliant, + "NuclearHessian must implement the DataClass interface"); + +} // namespace qdk::chemistry::data \ No newline at end of file diff --git a/cpp/src/qdk/chemistry/algorithms/CMakeLists.txt b/cpp/src/qdk/chemistry/algorithms/CMakeLists.txt index f28aa2d61..c06f2abae 100644 --- a/cpp/src/qdk/chemistry/algorithms/CMakeLists.txt +++ b/cpp/src/qdk/chemistry/algorithms/CMakeLists.txt @@ -4,6 +4,7 @@ target_sources(chemistry PRIVATE hamiltonian.cpp localization.cpp mc.cpp + nuclear_derivative.cpp pmc.cpp dynamical_correlation_calculator.cpp scf.cpp diff --git a/cpp/src/qdk/chemistry/algorithms/algorithm_defaults.cpp b/cpp/src/qdk/chemistry/algorithms/algorithm_defaults.cpp index b77316bbf..cbc009ac6 100644 --- a/cpp/src/qdk/chemistry/algorithms/algorithm_defaults.cpp +++ b/cpp/src/qdk/chemistry/algorithms/algorithm_defaults.cpp @@ -9,6 +9,7 @@ #include #include #include +#include #include #include #include @@ -41,6 +42,7 @@ std::shared_ptr resolve_algorithm_defaults( REGISTER_FACTORY_SETTINGS_INIT(ProjectedMultiConfigurationCalculatorFactory) REGISTER_FACTORY_SETTINGS_INIT(DynamicalCorrelationCalculatorFactory) REGISTER_FACTORY_SETTINGS_INIT(MultiConfigurationScfFactory) + REGISTER_FACTORY_SETTINGS_INIT(NuclearDerivativeCalculatorFactory) REGISTER_FACTORY_SETTINGS_INIT(LocalizerFactory) REGISTER_FACTORY_SETTINGS_INIT(StabilityCheckerFactory) diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf.cpp index 885a5d66d..e6b7c9bfc 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf.cpp @@ -39,9 +39,9 @@ std::pair calculate_electron_counts(int nuclear_charge, int charge, return {n_alpha, n_beta}; } -std::pair> ScfSolver::_run_impl( +ScfCalculationResult ScfSolver::_run_with_options( std::shared_ptr structure, int charge, int multiplicity, - BasisOrGuessType basis_or_guess) const { + BasisOrGuessType basis_or_guess, bool require_gradient) const { QDK_LOG_TRACE_ENTERING(); // Initialize the backend if not already done utils::microsoft::initialize_backend(); @@ -182,7 +182,7 @@ std::pair> ScfSolver::_run_impl( // Create SCFConfig auto ms_scf_config = std::make_unique(); ms_scf_config->mpi = qcs::mpi_default_input(); - ms_scf_config->require_gradient = false; + ms_scf_config->require_gradient = require_gradient; ms_scf_config->require_polarizability = false; ms_scf_config->exc.xc_name = method; std::transform(ms_scf_config->exc.xc_name.begin(), @@ -502,7 +502,38 @@ std::pair> ScfSolver::_run_impl( // Return total energy double total_energy = context.result.scf_total_energy; - return std::make_pair(total_energy, std::make_shared( - std::move(wavefunction))); + std::optional gradient; + if (require_gradient) { + const auto& raw_gradient = context.result.scf_total_gradient; + const auto expected_size = 3 * structure->get_num_atoms(); + if (raw_gradient.size() != expected_size) { + throw std::runtime_error( + "Internal SCF did not return the requested analytic nuclear " + "gradient"); + } + gradient = Eigen::Map( + raw_gradient.data(), static_cast(raw_gradient.size())); + } + + auto wavefunction_ptr = + std::make_shared(std::move(wavefunction)); + + return {total_energy, wavefunction_ptr, gradient}; +} + +ScfCalculationResult ScfSolver::run_with_analytic_gradient( + std::shared_ptr structure, int charge, int multiplicity, + BasisOrGuessType basis_or_guess) const { + this->lock_settings(); + return _run_with_options(structure, charge, multiplicity, basis_or_guess, + true); +} + +std::pair> ScfSolver::_run_impl( + std::shared_ptr structure, int charge, int multiplicity, + BasisOrGuessType basis_or_guess) const { + auto result = + _run_with_options(structure, charge, multiplicity, basis_or_guess, false); + return {result.energy, result.wavefunction}; } } // namespace qdk::chemistry::algorithms::microsoft diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf.hpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf.hpp index 0a06b79ca..e70788714 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf.hpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf.hpp @@ -4,10 +4,20 @@ #pragma once +#include #include namespace qdk::chemistry::algorithms::microsoft { +/** + * @brief Internal SCF result including optional analytic nuclear gradients. + */ +struct ScfCalculationResult { + double energy; + std::shared_ptr wavefunction; + std::optional nuclear_gradient; +}; + /** * @class ScfSettings * @brief Settings container for the internal SCF solver implementation @@ -122,6 +132,16 @@ class ScfSolver : public qdk::chemistry::algorithms::ScfSolver { virtual std::string name() const final { return "qdk"; } + /** + * @brief Run the internal SCF solver and return analytic nuclear gradients. + * + * Settings are locked in the same way as the base run() API. The returned + * gradient vector is atom-major with length 3 * number of atoms. + */ + ScfCalculationResult run_with_analytic_gradient( + std::shared_ptr structure, int charge, + int spin_multiplicity, BasisOrGuessType basis_or_guess) const; + protected: /** * @brief Perform an SCF calculation on the given molecular structure @@ -154,6 +174,12 @@ class ScfSolver : public qdk::chemistry::algorithms::ScfSolver { std::pair> _run_impl( std::shared_ptr structure, int charge, int spin_multiplicity, BasisOrGuessType basis_or_guess) const override; + + private: + ScfCalculationResult _run_with_options( + std::shared_ptr structure, int charge, + int spin_multiplicity, BasisOrGuessType basis_or_guess, + bool require_gradient) const; }; } // namespace qdk::chemistry::algorithms::microsoft diff --git a/cpp/src/qdk/chemistry/algorithms/nuclear_derivative.cpp b/cpp/src/qdk/chemistry/algorithms/nuclear_derivative.cpp new file mode 100644 index 000000000..9947a7ca8 --- /dev/null +++ b/cpp/src/qdk/chemistry/algorithms/nuclear_derivative.cpp @@ -0,0 +1,328 @@ +// Copyright (c) Microsoft Corporation. All rights reserved. +// Licensed under the MIT License. See LICENSE.txt in the project root for +// license information. + +#include +#include +#include +#include +#include +#include +#include + +#include "microsoft/scf.hpp" + +namespace qdk::chemistry::algorithms { + +namespace { + +class ScopedLogLevel { + public: + explicit ScopedLogLevel(utils::LogLevel minimum_level) + : previous_level_(utils::Logger::get_global_level()) { + if (static_cast(previous_level_) < static_cast(minimum_level)) { + utils::Logger::set_global_level(minimum_level); + changed_ = true; + } + } + + ~ScopedLogLevel() { + if (changed_) { + utils::Logger::set_global_level(previous_level_); + } + } + + ScopedLogLevel(const ScopedLogLevel&) = delete; + ScopedLogLevel& operator=(const ScopedLogLevel&) = delete; + + private: + utils::LogLevel previous_level_; + bool changed_ = false; +}; + +struct EnergyEvaluation { + double energy = 0.0; + std::optional> wavefunction; +}; + +std::shared_ptr copy_structure( + const std::shared_ptr& structure) { + if (!structure) { + throw std::invalid_argument("Structure must not be null"); + } + return std::make_shared(*structure); +} + +std::shared_ptr displace_structure( + const std::shared_ptr& structure, Eigen::Index coordinate, + double displacement) { + Eigen::MatrixXd coordinates = structure->get_coordinates(); + const Eigen::Index atom = coordinate / 3; + const Eigen::Index component = coordinate % 3; + coordinates(atom, component) += displacement; + return std::make_shared( + coordinates, structure->get_elements(), structure->get_masses(), + structure->get_nuclear_charges()); +} + +BasisOrGuessType seed_to_scf_input(const NuclearDerivativeSeedType& seed, + bool allow_orbital_guess) { + if (std::holds_alternative(seed)) { + return std::get(seed); + } + if (std::holds_alternative>(seed)) { + auto basis = std::get>(seed); + if (!basis) { + throw std::invalid_argument("Basis set seed must not be null"); + } + if (allow_orbital_guess) { + return basis; + } + if (basis->get_name() == data::BasisSet::custom_name) { + throw std::invalid_argument( + "Custom basis set cannot be reused for displaced numeric " + "derivatives"); + } + return basis->get_name(); + } + std::shared_ptr orbitals; + if (std::holds_alternative>(seed)) { + orbitals = std::get>(seed); + } else if (std::holds_alternative>( + seed)) { + auto wavefunction = std::get>(seed); + if (!wavefunction) { + throw std::invalid_argument("Wavefunction seed must not be null"); + } + orbitals = wavefunction->get_orbitals(); + } + + if (!orbitals || !orbitals->get_basis_set()) { + throw std::invalid_argument( + "Orbital or wavefunction seed must include orbitals with a basis set"); + } + if (allow_orbital_guess) { + return orbitals; + } + + auto basis = orbitals->get_basis_set(); + if (basis->get_name() == data::BasisSet::custom_name) { + throw std::invalid_argument( + "Custom orbital basis cannot be reused for displaced numeric " + "derivatives"); + } + return basis->get_name(); +} + +template +typename Factory::return_type create_from_ref(const data::AlgorithmRef& ref) { + auto instance = Factory::create(ref.get_algorithm_name()); + if (ref.get_settings()) { + instance->settings().update(*ref.get_settings()); + } + return instance; +} + +unsigned int active_electrons(const data::Settings& settings, + const std::string& key) { + auto value = settings.get(key); + if (value <= 0) { + throw std::invalid_argument( + key + " must be set to a positive value for this energy calculator"); + } + return static_cast(value); +} + +EnergyEvaluation evaluate_energy(const data::Settings& settings, + std::shared_ptr structure, + int charge, int spin_multiplicity, + const NuclearDerivativeSeedType& seed, + bool allow_wavefunction_seed) { + const auto ref = settings.get("energy_calculator"); + const auto& algorithm_type = ref.get_algorithm_type(); + + if (algorithm_type == ScfSolverFactory::algorithm_type_name()) { + auto solver = create_from_ref(ref); + auto [energy, wavefunction] = + solver->run(structure, charge, spin_multiplicity, + seed_to_scf_input(seed, allow_wavefunction_seed)); + return {energy, wavefunction}; + } + + if (algorithm_type == MultiConfigurationScfFactory::algorithm_type_name()) { + auto orbital_solver = create_from_ref( + settings.get("orbital_solver")); + auto [_, reference_wavefunction] = + orbital_solver->run(structure, charge, spin_multiplicity, + seed_to_scf_input(seed, allow_wavefunction_seed)); + auto calculator = create_from_ref(ref); + auto [energy, wavefunction] = + calculator->run(reference_wavefunction->get_orbitals(), + active_electrons(settings, "n_active_alpha_electrons"), + active_electrons(settings, "n_active_beta_electrons")); + return {energy, wavefunction}; + } + + if (algorithm_type == + MultiConfigurationCalculatorFactory::algorithm_type_name()) { + auto orbital_solver = create_from_ref( + settings.get("orbital_solver")); + auto [_, reference_wavefunction] = + orbital_solver->run(structure, charge, spin_multiplicity, + seed_to_scf_input(seed, allow_wavefunction_seed)); + auto hamiltonian_constructor = + create_from_ref( + settings.get("hamiltonian_constructor")); + auto hamiltonian = + hamiltonian_constructor->run(reference_wavefunction->get_orbitals()); + auto calculator = create_from_ref(ref); + auto [energy, wavefunction] = calculator->run( + hamiltonian, active_electrons(settings, "n_active_alpha_electrons"), + active_electrons(settings, "n_active_beta_electrons")); + return {energy, wavefunction}; + } + + throw std::invalid_argument( + "Unsupported energy_calculator algorithm type for numeric nuclear " + "derivatives: " + + algorithm_type); +} + +} // namespace + +NuclearDerivativeResult FiniteDifferenceNuclearDerivativeCalculator::_run_impl( + std::shared_ptr structure, int charge, + int spin_multiplicity, NuclearDerivativeSeedType seed) const { + const ScopedLogLevel scoped_log_level(utils::LogLevel::error); + const double step = _settings->get("finite_difference_step"); + const bool compute_hessian = _settings->get("compute_hessian"); + const auto dimension = + static_cast(3 * structure->get_num_atoms()); + + auto central = evaluate_energy(*_settings, structure, charge, + spin_multiplicity, seed, true); + Eigen::VectorXd gradients = Eigen::VectorXd::Zero(dimension); + std::vector plus_energies(static_cast(dimension)); + std::vector minus_energies(static_cast(dimension)); + + for (Eigen::Index coordinate = 0; coordinate < dimension; ++coordinate) { + auto plus_structure = displace_structure(structure, coordinate, step); + auto minus_structure = displace_structure(structure, coordinate, -step); + plus_energies[coordinate] = + evaluate_energy(*_settings, plus_structure, charge, spin_multiplicity, + seed, false) + .energy; + minus_energies[coordinate] = + evaluate_energy(*_settings, minus_structure, charge, spin_multiplicity, + seed, false) + .energy; + gradients(coordinate) = + (plus_energies[coordinate] - minus_energies[coordinate]) / (2.0 * step); + } + + std::optional> hessian; + if (compute_hessian) { + Eigen::MatrixXd hessian_matrix = + Eigen::MatrixXd::Zero(dimension, dimension); + for (Eigen::Index coordinate = 0; coordinate < dimension; ++coordinate) { + hessian_matrix(coordinate, coordinate) = + (plus_energies[coordinate] - 2.0 * central.energy + + minus_energies[coordinate]) / + (step * step); + } + + for (Eigen::Index i = 0; i < dimension; ++i) { + for (Eigen::Index j = i + 1; j < dimension; ++j) { + auto pp = + displace_structure(displace_structure(structure, i, step), j, step); + auto pm = displace_structure(displace_structure(structure, i, step), j, + -step); + auto mp = displace_structure(displace_structure(structure, i, -step), j, + step); + auto mm = displace_structure(displace_structure(structure, i, -step), j, + -step); + const double value = (evaluate_energy(*_settings, pp, charge, + spin_multiplicity, seed, false) + .energy - + evaluate_energy(*_settings, pm, charge, + spin_multiplicity, seed, false) + .energy - + evaluate_energy(*_settings, mp, charge, + spin_multiplicity, seed, false) + .energy + + evaluate_energy(*_settings, mm, charge, + spin_multiplicity, seed, false) + .energy) / + (4.0 * step * step); + hessian_matrix(i, j) = value; + hessian_matrix(j, i) = value; + } + } + if (_settings->get("symmetrize_hessian")) { + hessian_matrix = 0.5 * (hessian_matrix + hessian_matrix.transpose()); + } + hessian = std::make_shared(copy_structure(structure), + hessian_matrix); + } + + return {central.energy, + std::make_shared(copy_structure(structure), + gradients), + hessian, central.wavefunction}; +} + +std::unique_ptr +make_finite_difference_nuclear_derivative_calculator() { + return std::make_unique(); +} + +std::unique_ptr +make_qdk_nuclear_derivative_calculator() { + return std::make_unique(); +} + +NuclearDerivativeResult QdkNuclearDerivativeCalculator::_run_impl( + std::shared_ptr structure, int charge, + int spin_multiplicity, NuclearDerivativeSeedType seed) const { + const ScopedLogLevel scoped_log_level(utils::LogLevel::error); + if (_settings->get("compute_hessian")) { + throw std::invalid_argument( + "The QDK analytic nuclear derivative calculator does not currently " + "provide Hessians. Use the finite_difference implementation for " + "numeric Hessians."); + } + + const auto ref = _settings->get("energy_calculator"); + if (ref.get_algorithm_type() != ScfSolverFactory::algorithm_type_name() || + ref.get_algorithm_name() != "qdk") { + throw std::invalid_argument( + "The QDK analytic nuclear derivative calculator requires " + "energy_calculator to reference scf_solver/qdk."); + } + + microsoft::ScfSolver solver; + if (ref.get_settings()) { + solver.settings().update(*ref.get_settings()); + } + + auto scf_result = solver.run_with_analytic_gradient( + structure, charge, spin_multiplicity, seed_to_scf_input(seed, true)); + if (!scf_result.nuclear_gradient.has_value()) { + throw std::runtime_error( + "Internal SCF did not return an analytic nuclear gradient"); + } + + return {scf_result.energy, + std::make_shared( + copy_structure(structure), *scf_result.nuclear_gradient), + std::nullopt, scf_result.wavefunction}; +} + +void NuclearDerivativeCalculatorFactory::register_default_instances() { + NuclearDerivativeCalculatorFactory::register_instance( + &make_finite_difference_nuclear_derivative_calculator); + NuclearDerivativeCalculatorFactory::register_instance( + &make_qdk_nuclear_derivative_calculator); +} + +} // namespace qdk::chemistry::algorithms \ No newline at end of file diff --git a/cpp/src/qdk/chemistry/data/CMakeLists.txt b/cpp/src/qdk/chemistry/data/CMakeLists.txt index 5d04df8d4..13a705136 100644 --- a/cpp/src/qdk/chemistry/data/CMakeLists.txt +++ b/cpp/src/qdk/chemistry/data/CMakeLists.txt @@ -12,6 +12,8 @@ target_sources(chemistry PRIVATE hamiltonian_containers/cholesky.cpp hamiltonian_containers/sparse.cpp orbitals.cpp + nuclear_gradients.cpp + nuclear_hessian.cpp pauli_operator.cpp settings.cpp stability_result.cpp diff --git a/cpp/src/qdk/chemistry/data/nuclear_gradients.cpp b/cpp/src/qdk/chemistry/data/nuclear_gradients.cpp new file mode 100644 index 000000000..33c7cbc5a --- /dev/null +++ b/cpp/src/qdk/chemistry/data/nuclear_gradients.cpp @@ -0,0 +1,222 @@ +// Copyright (c) Microsoft Corporation. All rights reserved. +// Licensed under the MIT License. See LICENSE.txt in the project root for +// license information. + +#include + +#include +#include +#include +#include +#include + +#include "filename_utils.hpp" +#include "hdf5_error_handling.hpp" +#include "hdf5_serialization.hpp" +#include "json_serialization.hpp" + +namespace qdk::chemistry::data { + +NuclearGradients::NuclearGradients(std::shared_ptr structure, + const Eigen::VectorXd& values) + : structure_(std::move(structure)), values_(values) { + if (!_is_valid()) { + throw std::invalid_argument( + "NuclearGradients require a non-null structure and 3*N values"); + } +} + +const std::shared_ptr NuclearGradients::get_structure() const { + if (!structure_) { + throw std::runtime_error("No structure is associated with these gradients"); + } + return structure_; +} + +Eigen::MatrixXd NuclearGradients::as_matrix() const { + Eigen::MatrixXd matrix(get_num_atoms(), 3); + for (Eigen::Index atom = 0; atom < static_cast(get_num_atoms()); + ++atom) { + matrix(atom, 0) = values_(3 * atom); + matrix(atom, 1) = values_(3 * atom + 1); + matrix(atom, 2) = values_(3 * atom + 2); + } + return matrix; +} + +size_t NuclearGradients::get_num_atoms() const { + return get_structure()->get_num_atoms(); +} + +std::string NuclearGradients::get_summary() const { + std::ostringstream oss; + oss << "NuclearGradients(" << get_num_atoms() + << " atoms, norm=" << values_.norm() << ")"; + return oss.str(); +} + +void NuclearGradients::to_file(const std::string& filename, + const std::string& type) const { + if (type == "json") { + to_json_file(filename); + } else if (type == "hdf5") { + to_hdf5_file(filename); + } else { + throw std::invalid_argument("Unsupported file type: " + type + + ". Supported types are: json, hdf5"); + } +} + +nlohmann::json NuclearGradients::to_json() const { + nlohmann::json j; + j["serialization_version"] = SERIALIZATION_VERSION; + j["type"] = "NuclearGradients"; + j["structure"] = get_structure()->to_json(); + j["values"] = vector_to_json(values_); + j["units"] = "hartree/bohr"; + return j; +} + +void NuclearGradients::to_json_file(const std::string& filename) const { + auto validated_filename = DataTypeFilename::validate_write_suffix( + filename, DATACLASS_TO_SNAKE_CASE(NuclearGradients)); + _to_json_file(validated_filename); +} + +void NuclearGradients::to_hdf5(H5::Group& group) const { + auto scalar_space = H5::DataSpace(H5S_SCALAR); + auto str_type = H5::StrType(H5::PredType::C_S1, H5T_VARIABLE); + + auto version_attr = + group.createAttribute("serialization_version", str_type, scalar_space); + std::string version_str(SERIALIZATION_VERSION); + version_attr.write(str_type, version_str); + + auto type_attr = group.createAttribute("type", str_type, scalar_space); + std::string type_str = "NuclearGradients"; + type_attr.write(str_type, type_str); + + auto units_attr = group.createAttribute("units", str_type, scalar_space); + std::string units_str = "hartree/bohr"; + units_attr.write(str_type, units_str); + + auto structure_group = group.createGroup("structure"); + get_structure()->to_hdf5(structure_group); + save_vector_to_group(group, "values", values_); +} + +void NuclearGradients::to_hdf5_file(const std::string& filename) const { + auto validated_filename = DataTypeFilename::validate_write_suffix( + filename, DATACLASS_TO_SNAKE_CASE(NuclearGradients)); + _to_hdf5_file(validated_filename); +} + +std::shared_ptr NuclearGradients::from_file( + const std::string& filename, const std::string& type) { + if (type == "json") { + return from_json_file(filename); + } else if (type == "hdf5") { + return from_hdf5_file(filename); + } + throw std::invalid_argument("Unsupported file type: " + type + + ". Supported types are: json, hdf5"); +} + +std::shared_ptr NuclearGradients::from_json_file( + const std::string& filename) { + auto validated_filename = + DataTypeFilename::validate_read_suffix(filename, "nuclear_gradients"); + return _from_json_file(validated_filename); +} + +std::shared_ptr NuclearGradients::from_json( + const nlohmann::json& j) { + if (j.contains("serialization_version")) { + validate_serialization_version( + SERIALIZATION_VERSION, j["serialization_version"].get()); + } + if (j.contains("type") && + j["type"].get() != "NuclearGradients") { + throw std::runtime_error("Invalid type in JSON data"); + } + if (!j.contains("structure") || !j.contains("values")) { + throw std::runtime_error( + "NuclearGradients JSON requires structure and values fields"); + } + + return std::make_shared( + Structure::from_json(j["structure"]), json_to_vector(j["values"])); +} + +std::shared_ptr NuclearGradients::from_hdf5_file( + const std::string& filename) { + auto validated_filename = + DataTypeFilename::validate_read_suffix(filename, "nuclear_gradients"); + return _from_hdf5_file(validated_filename); +} + +std::shared_ptr NuclearGradients::from_hdf5( + H5::Group& group) { + if (group.attrExists("serialization_version")) { + auto attr = group.openAttribute("serialization_version"); + auto str_type = H5::StrType(H5::PredType::C_S1, H5T_VARIABLE); + std::string version; + attr.read(str_type, version); + validate_serialization_version(SERIALIZATION_VERSION, version); + } + if (group.attrExists("type")) { + auto attr = group.openAttribute("type"); + auto str_type = H5::StrType(H5::PredType::C_S1, H5T_VARIABLE); + std::string type; + attr.read(str_type, type); + if (type != "NuclearGradients") { + throw std::runtime_error("Invalid type in HDF5 data"); + } + } + auto structure_group = group.openGroup("structure"); + return std::make_shared( + Structure::from_hdf5(structure_group), + load_vector_from_group(group, "values")); +} + +bool NuclearGradients::_is_valid() const { + return structure_ != nullptr && + values_.size() == + static_cast(3 * structure_->get_num_atoms()); +} + +void NuclearGradients::_to_json_file(const std::string& filename) const { + std::ofstream file(filename); + if (!file.is_open()) { + throw std::runtime_error("Failed to open file for writing: " + filename); + } + file << to_json().dump(2); +} + +void NuclearGradients::_to_hdf5_file(const std::string& filename) const { + H5::H5File file(filename, H5F_ACC_TRUNC); + to_hdf5(file); +} + +std::shared_ptr NuclearGradients::_from_json_file( + const std::string& filename) { + std::ifstream file(filename); + if (!file.is_open()) { + throw std::runtime_error("Unable to open NuclearGradients JSON file: " + + filename); + } + nlohmann::json j; + file >> j; + return from_json(j); +} + +std::shared_ptr NuclearGradients::_from_hdf5_file( + const std::string& filename) { + if (hdf5_errors_should_be_suppressed()) { + H5::Exception::dontPrint(); + } + H5::H5File file(filename, H5F_ACC_RDONLY); + return from_hdf5(file); +} + +} // namespace qdk::chemistry::data \ No newline at end of file diff --git a/cpp/src/qdk/chemistry/data/nuclear_hessian.cpp b/cpp/src/qdk/chemistry/data/nuclear_hessian.cpp new file mode 100644 index 000000000..b8ed4d606 --- /dev/null +++ b/cpp/src/qdk/chemistry/data/nuclear_hessian.cpp @@ -0,0 +1,212 @@ +// Copyright (c) Microsoft Corporation. All rights reserved. +// Licensed under the MIT License. See LICENSE.txt in the project root for +// license information. + +#include + +#include +#include +#include +#include +#include + +#include "filename_utils.hpp" +#include "hdf5_error_handling.hpp" +#include "hdf5_serialization.hpp" +#include "json_serialization.hpp" + +namespace qdk::chemistry::data { + +NuclearHessian::NuclearHessian(std::shared_ptr structure, + const Eigen::MatrixXd& matrix) + : structure_(std::move(structure)), matrix_(matrix) { + if (!_is_valid()) { + throw std::invalid_argument( + "NuclearHessian requires a non-null structure and a 3*N by 3*N matrix"); + } +} + +const std::shared_ptr NuclearHessian::get_structure() const { + if (!structure_) { + throw std::runtime_error("No structure is associated with this Hessian"); + } + return structure_; +} + +size_t NuclearHessian::get_num_atoms() const { + return get_structure()->get_num_atoms(); +} + +std::string NuclearHessian::get_summary() const { + std::ostringstream oss; + oss << "NuclearHessian(" << get_num_atoms() << " atoms, " << matrix_.rows() + << "x" << matrix_.cols() << ")"; + return oss.str(); +} + +void NuclearHessian::to_file(const std::string& filename, + const std::string& type) const { + if (type == "json") { + to_json_file(filename); + } else if (type == "hdf5") { + to_hdf5_file(filename); + } else { + throw std::invalid_argument("Unsupported file type: " + type + + ". Supported types are: json, hdf5"); + } +} + +nlohmann::json NuclearHessian::to_json() const { + nlohmann::json j; + j["serialization_version"] = SERIALIZATION_VERSION; + j["type"] = "NuclearHessian"; + j["structure"] = get_structure()->to_json(); + j["matrix"] = matrix_to_json(matrix_); + j["units"] = "hartree/bohr^2"; + return j; +} + +void NuclearHessian::to_json_file(const std::string& filename) const { + auto validated_filename = DataTypeFilename::validate_write_suffix( + filename, DATACLASS_TO_SNAKE_CASE(NuclearHessian)); + _to_json_file(validated_filename); +} + +void NuclearHessian::to_hdf5(H5::Group& group) const { + auto scalar_space = H5::DataSpace(H5S_SCALAR); + auto str_type = H5::StrType(H5::PredType::C_S1, H5T_VARIABLE); + + auto version_attr = + group.createAttribute("serialization_version", str_type, scalar_space); + std::string version_str(SERIALIZATION_VERSION); + version_attr.write(str_type, version_str); + + auto type_attr = group.createAttribute("type", str_type, scalar_space); + std::string type_str = "NuclearHessian"; + type_attr.write(str_type, type_str); + + auto units_attr = group.createAttribute("units", str_type, scalar_space); + std::string units_str = "hartree/bohr^2"; + units_attr.write(str_type, units_str); + + auto structure_group = group.createGroup("structure"); + get_structure()->to_hdf5(structure_group); + save_matrix_to_group(group, "matrix", matrix_); +} + +void NuclearHessian::to_hdf5_file(const std::string& filename) const { + auto validated_filename = DataTypeFilename::validate_write_suffix( + filename, DATACLASS_TO_SNAKE_CASE(NuclearHessian)); + _to_hdf5_file(validated_filename); +} + +std::shared_ptr NuclearHessian::from_file( + const std::string& filename, const std::string& type) { + if (type == "json") { + return from_json_file(filename); + } else if (type == "hdf5") { + return from_hdf5_file(filename); + } + throw std::invalid_argument("Unsupported file type: " + type + + ". Supported types are: json, hdf5"); +} + +std::shared_ptr NuclearHessian::from_json_file( + const std::string& filename) { + auto validated_filename = + DataTypeFilename::validate_read_suffix(filename, "nuclear_hessian"); + return _from_json_file(validated_filename); +} + +std::shared_ptr NuclearHessian::from_json( + const nlohmann::json& j) { + if (j.contains("serialization_version")) { + validate_serialization_version( + SERIALIZATION_VERSION, j["serialization_version"].get()); + } + if (j.contains("type") && j["type"].get() != "NuclearHessian") { + throw std::runtime_error("Invalid type in JSON data"); + } + if (!j.contains("structure") || !j.contains("matrix")) { + throw std::runtime_error( + "NuclearHessian JSON requires structure and matrix fields"); + } + + return std::make_shared(Structure::from_json(j["structure"]), + json_to_matrix(j["matrix"])); +} + +std::shared_ptr NuclearHessian::from_hdf5_file( + const std::string& filename) { + auto validated_filename = + DataTypeFilename::validate_read_suffix(filename, "nuclear_hessian"); + return _from_hdf5_file(validated_filename); +} + +std::shared_ptr NuclearHessian::from_hdf5(H5::Group& group) { + if (group.attrExists("serialization_version")) { + auto attr = group.openAttribute("serialization_version"); + auto str_type = H5::StrType(H5::PredType::C_S1, H5T_VARIABLE); + std::string version; + attr.read(str_type, version); + validate_serialization_version(SERIALIZATION_VERSION, version); + } + if (group.attrExists("type")) { + auto attr = group.openAttribute("type"); + auto str_type = H5::StrType(H5::PredType::C_S1, H5T_VARIABLE); + std::string type; + attr.read(str_type, type); + if (type != "NuclearHessian") { + throw std::runtime_error("Invalid type in HDF5 data"); + } + } + auto structure_group = group.openGroup("structure"); + return std::make_shared( + Structure::from_hdf5(structure_group), + load_matrix_from_group(group, "matrix")); +} + +bool NuclearHessian::_is_valid() const { + if (!structure_ || matrix_.rows() != matrix_.cols()) { + return false; + } + const auto expected = + static_cast(3 * structure_->get_num_atoms()); + return matrix_.rows() == expected && matrix_.cols() == expected; +} + +void NuclearHessian::_to_json_file(const std::string& filename) const { + std::ofstream file(filename); + if (!file.is_open()) { + throw std::runtime_error("Failed to open file for writing: " + filename); + } + file << to_json().dump(2); +} + +void NuclearHessian::_to_hdf5_file(const std::string& filename) const { + H5::H5File file(filename, H5F_ACC_TRUNC); + to_hdf5(file); +} + +std::shared_ptr NuclearHessian::_from_json_file( + const std::string& filename) { + std::ifstream file(filename); + if (!file.is_open()) { + throw std::runtime_error("Unable to open NuclearHessian JSON file: " + + filename); + } + nlohmann::json j; + file >> j; + return from_json(j); +} + +std::shared_ptr NuclearHessian::_from_hdf5_file( + const std::string& filename) { + if (hdf5_errors_should_be_suppressed()) { + H5::Exception::dontPrint(); + } + H5::H5File file(filename, H5F_ACC_RDONLY); + return from_hdf5(file); +} + +} // namespace qdk::chemistry::data \ No newline at end of file diff --git a/cpp/tests/test_nuclear_derivative.cpp b/cpp/tests/test_nuclear_derivative.cpp new file mode 100644 index 000000000..f2f05cff0 --- /dev/null +++ b/cpp/tests/test_nuclear_derivative.cpp @@ -0,0 +1,127 @@ +// Copyright (c) Microsoft Corporation. All rights reserved. +// Licensed under the MIT License. See LICENSE.txt in the project root for +// license information. + +#include + +#include +#include +#include +#include +#include +#include + +using namespace qdk::chemistry::algorithms; +using namespace qdk::chemistry::data; + +namespace { + +std::shared_ptr create_h2_structure() { + std::vector coordinates = {{0.0, 0.0, 0.0}, {0.0, 0.0, 1.4}}; + std::vector elements = {Element::H, Element::H}; + return std::make_shared(coordinates, elements); +} + +void expect_same_structure(const std::shared_ptr& left, + const std::shared_ptr& right) { + ASSERT_NE(left, nullptr); + ASSERT_NE(right, nullptr); + ASSERT_EQ(left->get_num_atoms(), right->get_num_atoms()); + EXPECT_TRUE( + left->get_coordinates().isApprox(right->get_coordinates(), 1.0e-12)); + EXPECT_EQ(left->get_elements(), right->get_elements()); +} + +} // namespace + +TEST(NuclearDerivativeDataTest, GradientsReferenceStructureAndRoundTripJson) { + auto structure = create_h2_structure(); + Eigen::VectorXd values(6); + values << 0.0, 0.1, 0.2, 0.3, 0.4, 0.5; + + NuclearGradients gradients(structure, values); + + EXPECT_EQ(gradients.get_num_atoms(), 2); + EXPECT_TRUE(gradients.get_values().isApprox(values)); + EXPECT_EQ(gradients.as_matrix().rows(), 2); + EXPECT_EQ(gradients.as_matrix().cols(), 3); + expect_same_structure(structure, gradients.get_structure()); + + auto loaded = NuclearGradients::from_json(gradients.to_json()); + EXPECT_TRUE(loaded->get_values().isApprox(values)); + expect_same_structure(structure, loaded->get_structure()); +} + +TEST(NuclearDerivativeDataTest, HessianReferencesStructureAndRoundTripsJson) { + auto structure = create_h2_structure(); + Eigen::MatrixXd matrix = Eigen::MatrixXd::Identity(6, 6); + + NuclearHessian hessian(structure, matrix); + + EXPECT_EQ(hessian.get_num_atoms(), 2); + EXPECT_TRUE(hessian.get_matrix().isApprox(matrix)); + expect_same_structure(structure, hessian.get_structure()); + + auto loaded = NuclearHessian::from_json(hessian.to_json()); + EXPECT_TRUE(loaded->get_matrix().isApprox(matrix)); + expect_same_structure(structure, loaded->get_structure()); +} + +TEST(NuclearDerivativeCalculatorTest, FiniteDifferenceRunsRealScfForH2) { + auto calculator = NuclearDerivativeCalculatorFactory::create(); + calculator->settings().set("compute_hessian", true); + calculator->settings().set("finite_difference_step", 1.0e-3); + + auto [energy, gradients, hessian, wavefunction] = + calculator->run(create_h2_structure(), 0, 1, std::string("sto-3g")); + + ASSERT_TRUE(std::isfinite(energy)); + EXPECT_LT(energy, 0.0); + ASSERT_NE(gradients, nullptr); + EXPECT_EQ(gradients->get_num_atoms(), 2); + EXPECT_EQ(gradients->get_values().size(), 6); + EXPECT_TRUE(gradients->get_values().allFinite()); + + auto gradient_matrix = gradients->as_matrix(); + EXPECT_NEAR(gradient_matrix.col(0).sum(), 0.0, 1.0e-5); + EXPECT_NEAR(gradient_matrix.col(1).sum(), 0.0, 1.0e-5); + EXPECT_NEAR(gradient_matrix.col(2).sum(), 0.0, 1.0e-5); + + ASSERT_TRUE(hessian.has_value()); + ASSERT_NE(*hessian, nullptr); + EXPECT_EQ((*hessian)->get_num_atoms(), 2); + EXPECT_EQ((*hessian)->get_matrix().rows(), 6); + EXPECT_EQ((*hessian)->get_matrix().cols(), 6); + EXPECT_TRUE((*hessian)->get_matrix().allFinite()); + EXPECT_TRUE((*hessian)->get_matrix().isApprox( + (*hessian)->get_matrix().transpose(), 1.0e-8)); + + ASSERT_TRUE(wavefunction.has_value()); + EXPECT_NE(*wavefunction, nullptr); +} + +TEST(NuclearDerivativeCalculatorTest, QdkAnalyticGradientRunsRealDftForH2) { + auto calculator = NuclearDerivativeCalculatorFactory::create("qdk"); + AlgorithmRef scf_ref("scf_solver", "qdk"); + scf_ref.get_settings()->set("method", std::string("pbe")); + calculator->settings().set("energy_calculator", scf_ref); + + auto [energy, gradients, hessian, wavefunction] = + calculator->run(create_h2_structure(), 0, 1, std::string("sto-3g")); + + ASSERT_TRUE(std::isfinite(energy)); + EXPECT_LT(energy, 0.0); + ASSERT_NE(gradients, nullptr); + EXPECT_EQ(gradients->get_num_atoms(), 2); + EXPECT_EQ(gradients->get_values().size(), 6); + EXPECT_TRUE(gradients->get_values().allFinite()); + + auto gradient_matrix = gradients->as_matrix(); + EXPECT_NEAR(gradient_matrix.col(0).sum(), 0.0, 1.0e-8); + EXPECT_NEAR(gradient_matrix.col(1).sum(), 0.0, 1.0e-8); + EXPECT_NEAR(gradient_matrix.col(2).sum(), 0.0, 1.0e-8); + + EXPECT_FALSE(hessian.has_value()); + ASSERT_TRUE(wavefunction.has_value()); + EXPECT_NE(*wavefunction, nullptr); +} diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index c0a3d5ddd..e000fed5e 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -128,6 +128,8 @@ pybind11_add_module(_core src/pybind11/data/configuration_set.cpp src/pybind11/data/ansatz.cpp src/pybind11/data/stability_result.cpp + src/pybind11/data/nuclear_gradients.cpp + src/pybind11/data/nuclear_hessian.cpp src/pybind11/data/pauli_operator.cpp src/pybind11/data/serialization.cpp src/pybind11/data/lattice_graph.cpp @@ -138,6 +140,7 @@ pybind11_add_module(_core src/pybind11/algorithms/mcscf.cpp src/pybind11/algorithms/hamiltonian.cpp src/pybind11/algorithms/scf.cpp + src/pybind11/algorithms/nuclear_derivative.cpp src/pybind11/algorithms/active_space.cpp src/pybind11/algorithms/dynamical_correlation_calculator.cpp src/pybind11/algorithms/davidson_solver.cpp diff --git a/python/src/pybind11/algorithms/nuclear_derivative.cpp b/python/src/pybind11/algorithms/nuclear_derivative.cpp new file mode 100644 index 000000000..49ac3e7dd --- /dev/null +++ b/python/src/pybind11/algorithms/nuclear_derivative.cpp @@ -0,0 +1,118 @@ +// Copyright (c) Microsoft Corporation. All rights reserved. +// Licensed under the MIT License. See LICENSE.txt in the project root for +// license information. + +#include +#include + +#include + +#include "factory_bindings.hpp" + +namespace py = pybind11; +using namespace qdk::chemistry::algorithms; +using namespace qdk::chemistry::data; + +class NuclearDerivativeCalculatorBase + : public NuclearDerivativeCalculator, + public pybind11::trampoline_self_life_support { + public: + std::string name() const override { + PYBIND11_OVERRIDE_PURE(std::string, NuclearDerivativeCalculator, name); + } + + std::vector aliases() const override { + PYBIND11_OVERRIDE(std::vector, NuclearDerivativeCalculator, + aliases); + } + + void replace_settings( + std::unique_ptr new_settings) { + this->_settings = std::move(new_settings); + } + + protected: + NuclearDerivativeResult _run_impl( + std::shared_ptr structure, int charge, int spin_multiplicity, + NuclearDerivativeSeedType seed) const override { + PYBIND11_OVERRIDE_PURE(NuclearDerivativeResult, NuclearDerivativeCalculator, + _run_impl, structure, charge, spin_multiplicity, + seed); + } +}; + +void bind_nuclear_derivative(py::module& m) { + py::class_ + calculator(m, "NuclearDerivativeCalculator", + R"( + Base class for nuclear derivative algorithms. + + Calculators return the total energy, nuclear gradients, an optional Hessian, + and an optional wavefunction for a molecular structure. + )"); + + calculator.def(py::init<>(), R"(Create a nuclear derivative calculator.)"); + calculator.def( + "run", + [](const NuclearDerivativeCalculator& self, + std::shared_ptr structure, int charge, + int spin_multiplicity, NuclearDerivativeSeedType seed) { + return self.run(structure, charge, spin_multiplicity, seed); + }, + py::arg("structure"), py::arg("charge"), py::arg("spin_multiplicity"), + py::arg("seed_or_basis"), + R"( +Compute nuclear derivatives for a molecular structure. + +Args: + structure: Molecular structure to evaluate. + charge: Total molecular charge. + spin_multiplicity: Spin multiplicity of the molecular system. + seed_or_basis: Basis name, basis set, orbitals, or wavefunction seed. + +Returns: + tuple: ``(energy, gradients, hessian, wavefunction)``. +)"); + calculator.def("settings", &NuclearDerivativeCalculator::settings, + py::return_value_policy::reference_internal, + R"(Return the calculator settings.)"); + calculator.def_property( + "_settings", + [](NuclearDerivativeCalculatorBase& self) -> Settings& { + return self.settings(); + }, + [](NuclearDerivativeCalculatorBase& self, + std::unique_ptr new_settings) { + self.replace_settings(std::move(new_settings)); + }, + py::return_value_policy::reference_internal, + R"(Internal settings replacement hook for Python subclasses.)"); + calculator.def("name", &NuclearDerivativeCalculator::name, + R"(Return the implementation name.)"); + calculator.def("type_name", &NuclearDerivativeCalculator::type_name, + R"(Return the algorithm type name.)"); + calculator.def("__repr__", [](const NuclearDerivativeCalculator& self) { + return ""; + }); + + qdk::chemistry::python::bind_algorithm_factory< + NuclearDerivativeCalculatorFactory, NuclearDerivativeCalculator, + NuclearDerivativeCalculatorBase>(m, "NuclearDerivativeCalculatorFactory"); + qdk::chemistry::python::bind_create_nested(calculator); + + py::class_( + m, "FiniteDifferenceNuclearDerivativeCalculator", + R"(Numeric nuclear derivative calculator using central finite differences.)") + .def(py::init<>(), + R"(Create a finite-difference derivative calculator.)"); + + py::class_( + m, "QdkNuclearDerivativeCalculator", + R"(QDK nuclear derivative calculator using analytic internal SCF gradients.)") + .def(py::init<>(), + R"(Create a QDK analytic nuclear derivative calculator.)"); +} \ No newline at end of file diff --git a/python/src/pybind11/data/nuclear_gradients.cpp b/python/src/pybind11/data/nuclear_gradients.cpp new file mode 100644 index 000000000..665dedf61 --- /dev/null +++ b/python/src/pybind11/data/nuclear_gradients.cpp @@ -0,0 +1,125 @@ +// Copyright (c) Microsoft Corporation. All rights reserved. +// Licensed under the MIT License. See LICENSE.txt in the project root for +// license information. + +#include +#include +#include + +#include +#include +#include + +#include "path_utils.hpp" + +namespace py = pybind11; +using namespace qdk::chemistry::data; + +namespace { + +void nuclear_gradients_to_file(NuclearGradients& self, + const py::object& filename, + const std::string& format_type) { + self.to_file(qdk::chemistry::python::utils::to_string_path(filename), + format_type); +} + +void nuclear_gradients_to_json_file(NuclearGradients& self, + const py::object& filename) { + self.to_json_file(qdk::chemistry::python::utils::to_string_path(filename)); +} + +void nuclear_gradients_to_hdf5_file(NuclearGradients& self, + const py::object& filename) { + self.to_hdf5_file(qdk::chemistry::python::utils::to_string_path(filename)); +} + +std::shared_ptr nuclear_gradients_from_file( + const py::object& filename, const std::string& format_type) { + return NuclearGradients::from_file( + qdk::chemistry::python::utils::to_string_path(filename), format_type); +} + +std::shared_ptr nuclear_gradients_from_json_file( + const py::object& filename) { + return NuclearGradients::from_json_file( + qdk::chemistry::python::utils::to_string_path(filename)); +} + +std::shared_ptr nuclear_gradients_from_hdf5_file( + const py::object& filename) { + return NuclearGradients::from_hdf5_file( + qdk::chemistry::python::utils::to_string_path(filename)); +} + +} // namespace + +void bind_nuclear_gradients(py::module& m) { + py::class_(m, + "NuclearGradients", + R"( +Nuclear energy gradients for a molecular structure. + +Values are stored as an atom-major vector with x, y, z components for each +atom. The associated structure records the geometry for which the gradients +were computed. +)") + .def(py::init, const Eigen::VectorXd&>(), + R"( +Create nuclear gradients for a structure. + +Args: + structure: Molecular structure used to compute the gradients. + values: Atom-major gradient vector with length ``3 * num_atoms``. +)", + py::arg("structure"), py::arg("values")) + .def_property_readonly( + "structure", &NuclearGradients::get_structure, + R"(Molecular structure associated with the gradients.)") + .def_property_readonly( + "values", + [](const NuclearGradients& gradients) -> const Eigen::VectorXd& { + return gradients.get_values(); + }, + py::return_value_policy::reference_internal, + R"(Atom-major gradient vector in Hartree/Bohr.)") + .def("get_structure", &NuclearGradients::get_structure, + R"(Return the molecular structure associated with the gradients.)") + .def("get_values", &NuclearGradients::get_values, + py::return_value_policy::reference_internal, + R"(Return the atom-major gradient vector in Hartree/Bohr.)") + .def("as_matrix", &NuclearGradients::as_matrix, + R"(Return gradients as an ``(num_atoms, 3)`` matrix.)") + .def("get_num_atoms", &NuclearGradients::get_num_atoms, + R"(Return the number of atoms in the associated structure.)") + .def("get_data_type_name", &NuclearGradients::get_data_type_name, + R"(Return the serialized data type name.)") + .def("get_summary", &NuclearGradients::get_summary, + R"(Return a short human-readable summary.)") + .def("to_file", nuclear_gradients_to_file, py::arg("filename"), + py::arg("format_type"), R"(Save gradients to a JSON or HDF5 file.)") + .def( + "to_json", + [](const NuclearGradients& self) { return self.to_json().dump(2); }, + R"(Serialize gradients to a JSON string.)") + .def("to_json_file", nuclear_gradients_to_json_file, py::arg("filename"), + R"(Save gradients to a JSON file.)") + .def("to_hdf5_file", nuclear_gradients_to_hdf5_file, py::arg("filename"), + R"(Save gradients to an HDF5 file.)") + .def_static("from_file", nuclear_gradients_from_file, py::arg("filename"), + py::arg("format_type"), + R"(Load gradients from a JSON or HDF5 file.)") + .def_static( + "from_json", + [](const std::string& json_str) { + return NuclearGradients::from_json(nlohmann::json::parse(json_str)); + }, + py::arg("json_str"), R"(Load gradients from a JSON string.)") + .def_static("from_json_file", nuclear_gradients_from_json_file, + py::arg("filename"), R"(Load gradients from a JSON file.)") + .def_static("from_hdf5_file", nuclear_gradients_from_hdf5_file, + py::arg("filename"), R"(Load gradients from an HDF5 file.)") + .def("__repr__", [](const NuclearGradients& gradients) { + return gradients.get_summary(); + }); +} \ No newline at end of file diff --git a/python/src/pybind11/data/nuclear_hessian.cpp b/python/src/pybind11/data/nuclear_hessian.cpp new file mode 100644 index 000000000..d4cd2dbea --- /dev/null +++ b/python/src/pybind11/data/nuclear_hessian.cpp @@ -0,0 +1,120 @@ +// Copyright (c) Microsoft Corporation. All rights reserved. +// Licensed under the MIT License. See LICENSE.txt in the project root for +// license information. + +#include +#include +#include + +#include +#include +#include + +#include "path_utils.hpp" + +namespace py = pybind11; +using namespace qdk::chemistry::data; + +namespace { + +void nuclear_hessian_to_file(NuclearHessian& self, const py::object& filename, + const std::string& format_type) { + self.to_file(qdk::chemistry::python::utils::to_string_path(filename), + format_type); +} + +void nuclear_hessian_to_json_file(NuclearHessian& self, + const py::object& filename) { + self.to_json_file(qdk::chemistry::python::utils::to_string_path(filename)); +} + +void nuclear_hessian_to_hdf5_file(NuclearHessian& self, + const py::object& filename) { + self.to_hdf5_file(qdk::chemistry::python::utils::to_string_path(filename)); +} + +std::shared_ptr nuclear_hessian_from_file( + const py::object& filename, const std::string& format_type) { + return NuclearHessian::from_file( + qdk::chemistry::python::utils::to_string_path(filename), format_type); +} + +std::shared_ptr nuclear_hessian_from_json_file( + const py::object& filename) { + return NuclearHessian::from_json_file( + qdk::chemistry::python::utils::to_string_path(filename)); +} + +std::shared_ptr nuclear_hessian_from_hdf5_file( + const py::object& filename) { + return NuclearHessian::from_hdf5_file( + qdk::chemistry::python::utils::to_string_path(filename)); +} + +} // namespace + +void bind_nuclear_hessian(py::module& m) { + py::class_(m, "NuclearHessian", + R"( +Nuclear second derivatives for a molecular structure. + +The matrix is ordered atom-major by x, y, z components and records the +geometry for which the Hessian was computed. +)") + .def(py::init, const Eigen::MatrixXd&>(), + R"( +Create a nuclear Hessian for a structure. + +Args: + structure: Molecular structure used to compute the Hessian. + matrix: Square ``3N x 3N`` Hessian matrix in Hartree/Bohr^2. +)", + py::arg("structure"), py::arg("matrix")) + .def_property_readonly( + "structure", &NuclearHessian::get_structure, + R"(Molecular structure associated with the Hessian.)") + .def_property_readonly( + "matrix", + [](const NuclearHessian& hessian) -> const Eigen::MatrixXd& { + return hessian.get_matrix(); + }, + py::return_value_policy::reference_internal, + R"(Hessian matrix in Hartree/Bohr^2.)") + .def("get_structure", &NuclearHessian::get_structure, + R"(Return the molecular structure associated with the Hessian.)") + .def("get_matrix", &NuclearHessian::get_matrix, + py::return_value_policy::reference_internal, + R"(Return the Hessian matrix in Hartree/Bohr^2.)") + .def("get_num_atoms", &NuclearHessian::get_num_atoms, + R"(Return the number of atoms in the associated structure.)") + .def("get_data_type_name", &NuclearHessian::get_data_type_name, + R"(Return the serialized data type name.)") + .def("get_summary", &NuclearHessian::get_summary, + R"(Return a short human-readable summary.)") + .def("to_file", nuclear_hessian_to_file, py::arg("filename"), + py::arg("format_type"), + R"(Save the Hessian to a JSON or HDF5 file.)") + .def( + "to_json", + [](const NuclearHessian& self) { return self.to_json().dump(2); }, + R"(Serialize the Hessian to a JSON string.)") + .def("to_json_file", nuclear_hessian_to_json_file, py::arg("filename"), + R"(Save the Hessian to a JSON file.)") + .def("to_hdf5_file", nuclear_hessian_to_hdf5_file, py::arg("filename"), + R"(Save the Hessian to an HDF5 file.)") + .def_static("from_file", nuclear_hessian_from_file, py::arg("filename"), + py::arg("format_type"), + R"(Load a Hessian from a JSON or HDF5 file.)") + .def_static( + "from_json", + [](const std::string& json_str) { + return NuclearHessian::from_json(nlohmann::json::parse(json_str)); + }, + py::arg("json_str"), R"(Load a Hessian from a JSON string.)") + .def_static("from_json_file", nuclear_hessian_from_json_file, + py::arg("filename"), R"(Load a Hessian from a JSON file.)") + .def_static("from_hdf5_file", nuclear_hessian_from_hdf5_file, + py::arg("filename"), R"(Load a Hessian from an HDF5 file.)") + .def("__repr__", + [](const NuclearHessian& hessian) { return hessian.get_summary(); }); +} \ No newline at end of file diff --git a/python/src/pybind11/module.cpp b/python/src/pybind11/module.cpp index 23ba9ddfe..6c8ad1490 100644 --- a/python/src/pybind11/module.cpp +++ b/python/src/pybind11/module.cpp @@ -18,6 +18,8 @@ void bind_configuration_set(py::module& m); void bind_localizer(py::module& m); void bind_stability(py::module& m); void bind_stability_result(py::module& m); +void bind_nuclear_gradients(py::module& m); +void bind_nuclear_hessian(py::module& m); void bind_settings(py::module& m); void bind_structure(py::module& m); void bind_basis_set(py::module& m); @@ -26,6 +28,7 @@ void bind_mc(py::module& m); void bind_mcscf(py::module& m); void bind_hamiltonian_constructor(py::module& m); void bind_scf(py::module& m); +void bind_nuclear_derivative(py::module& m); void bind_active_space(py::module& m); void bind_constants(py::module& m); void bind_pmc(py::module& m); @@ -68,6 +71,8 @@ PYBIND11_MODULE(_core, m) { bind_wavefunction(data); bind_ansatz(data); bind_stability_result(data); + bind_nuclear_gradients(data); + bind_nuclear_hessian(data); bind_pauli_operator(data); bind_serialization(data); @@ -76,6 +81,7 @@ PYBIND11_MODULE(_core, m) { bind_mcscf(algorithms); bind_hamiltonian_constructor(algorithms); bind_scf(algorithms); + bind_nuclear_derivative(algorithms); bind_active_space(algorithms); bind_dynamical_correlation_calculator(algorithms); bind_pmc(algorithms); diff --git a/python/src/qdk_chemistry/algorithms/__init__.py b/python/src/qdk_chemistry/algorithms/__init__.py index c57d8b10d..cfb0c8d10 100644 --- a/python/src/qdk_chemistry/algorithms/__init__.py +++ b/python/src/qdk_chemistry/algorithms/__init__.py @@ -37,6 +37,11 @@ QdkMacisCas, ) from qdk_chemistry.algorithms.multi_configuration_scf import MultiConfigurationScf +from qdk_chemistry.algorithms.nuclear_derivative import ( + FiniteDifferenceNuclearDerivativeCalculator, + NuclearDerivativeCalculator, + QdkNuclearDerivativeCalculator, +) from qdk_chemistry.algorithms.orbital_localizer import ( OrbitalLocalizer, QdkMP2NaturalOrbitalLocalizer, @@ -63,10 +68,12 @@ "ControlledCircuitMapper", "DynamicalCorrelationCalculator", "EnergyEstimator", + "FiniteDifferenceNuclearDerivativeCalculator", "HamiltonianConstructor", "HamiltonianUnitaryBuilder", "MultiConfigurationCalculator", "MultiConfigurationScf", + "NuclearDerivativeCalculator", "OrbitalLocalizer", "PhaseEstimation", "ProjectedMultiConfigurationCalculator", @@ -80,6 +87,7 @@ "QdkMacisPmc", "QdkOccupationActiveSpaceSelector", "QdkPipekMezeyLocalizer", + "QdkNuclearDerivativeCalculator", "QdkQubitMapper", "QdkScfSolver", "QdkStabilityChecker", diff --git a/python/src/qdk_chemistry/algorithms/nuclear_derivative.py b/python/src/qdk_chemistry/algorithms/nuclear_derivative.py new file mode 100644 index 000000000..1758d15ea --- /dev/null +++ b/python/src/qdk_chemistry/algorithms/nuclear_derivative.py @@ -0,0 +1,18 @@ +"""Public entry point for nuclear derivative calculators.""" + +# -------------------------------------------------------------------------------------------- +# Copyright (c) Microsoft Corporation. All rights reserved. +# Licensed under the MIT License. See LICENSE.txt in the project root for license information. +# -------------------------------------------------------------------------------------------- + +from qdk_chemistry._core._algorithms import ( # noqa: F401 + FiniteDifferenceNuclearDerivativeCalculator, + NuclearDerivativeCalculator, + QdkNuclearDerivativeCalculator, +) + +__all__ = [ + "FiniteDifferenceNuclearDerivativeCalculator", + "NuclearDerivativeCalculator", + "QdkNuclearDerivativeCalculator", +] \ No newline at end of file diff --git a/python/src/qdk_chemistry/algorithms/registry.py b/python/src/qdk_chemistry/algorithms/registry.py index 49d55fc2b..b3fa9211e 100644 --- a/python/src/qdk_chemistry/algorithms/registry.py +++ b/python/src/qdk_chemistry/algorithms/registry.py @@ -477,6 +477,7 @@ def _register_cpp_factories(): LocalizerFactory, MultiConfigurationCalculatorFactory, MultiConfigurationScfFactory, + NuclearDerivativeCalculatorFactory, ProjectedMultiConfigurationCalculatorFactory, ScfSolverFactory, StabilityCheckerFactory, @@ -487,6 +488,7 @@ def _register_cpp_factories(): register_factory(LocalizerFactory) register_factory(MultiConfigurationCalculatorFactory) register_factory(MultiConfigurationScfFactory) + register_factory(NuclearDerivativeCalculatorFactory) register_factory(ProjectedMultiConfigurationCalculatorFactory) register_factory(DynamicalCorrelationCalculatorFactory) register_factory(ScfSolverFactory) diff --git a/python/src/qdk_chemistry/data/__init__.py b/python/src/qdk_chemistry/data/__init__.py index d7d17dbdc..861314471 100644 --- a/python/src/qdk_chemistry/data/__init__.py +++ b/python/src/qdk_chemistry/data/__init__.py @@ -29,6 +29,8 @@ - :class:`SparseHamiltonianContainer`: Container for lattice model Hamiltonians with sparse internal storage. - :class:`ModelOrbitals`: Simple orbital representation for model systems without full basis set information. - :class:`MP2Container`: Container for MP2 wavefunction with Hamiltonian reference and optional amplitudes. +- :class:`NuclearGradients`: Nuclear gradient values associated with a molecular structure. +- :class:`NuclearHessian`: Nuclear Hessian matrix associated with a molecular structure. - :class:`Orbitals`: Molecular orbital information and properties. - :class:`OrbitalType`: Enumeration of orbital angular momentum types (s, p, d, f, etc.). - :class:`PauliOperator`: Pauli operator (I, X, Y, Z) for quantum operator expressions with arithmetic support. @@ -85,6 +87,8 @@ LatticeGraph, ModelOrbitals, MP2Container, + NuclearGradients, + NuclearHessian, Orbitals, OrbitalType, PauliOperator, @@ -163,6 +167,8 @@ "MP2Container", "MeasurementData", "ModelOrbitals", + "NuclearGradients", + "NuclearHessian", "OrbitalType", "Orbitals", "PauliOperator", diff --git a/python/tests/test_nuclear_derivative.py b/python/tests/test_nuclear_derivative.py new file mode 100644 index 000000000..104f51ecf --- /dev/null +++ b/python/tests/test_nuclear_derivative.py @@ -0,0 +1,63 @@ +"""Tests for nuclear derivative data and algorithm bindings.""" + +# -------------------------------------------------------------------------------------------- +# Copyright (c) Microsoft Corporation. All rights reserved. +# Licensed under the MIT License. See LICENSE.txt in the project root for license information. +# -------------------------------------------------------------------------------------------- + +import numpy as np + +from qdk_chemistry import algorithms, data + + +def _h2_structure(): + return data.Structure([[0.0, 0.0, 0.0], [0.0, 0.0, 1.4]], [1, 1]) + + +def test_nuclear_gradients_roundtrip_json(): + structure = _h2_structure() + values = np.arange(6, dtype=float) + + gradients = data.NuclearGradients(structure, values) + + assert gradients.get_num_atoms() == 2 + np.testing.assert_allclose(gradients.get_values(), values) + assert gradients.as_matrix().shape == (2, 3) + + loaded = data.NuclearGradients.from_json(gradients.to_json()) + np.testing.assert_allclose(loaded.get_values(), values) + np.testing.assert_allclose( + loaded.get_structure().get_coordinates(), + structure.get_coordinates(), + ) + + +def test_nuclear_hessian_roundtrip_json(): + structure = _h2_structure() + matrix = np.eye(6) + + hessian = data.NuclearHessian(structure, matrix) + + assert hessian.get_num_atoms() == 2 + np.testing.assert_allclose(hessian.get_matrix(), matrix) + + loaded = data.NuclearHessian.from_json(hessian.to_json()) + np.testing.assert_allclose(loaded.get_matrix(), matrix) + np.testing.assert_allclose( + loaded.get_structure().get_coordinates(), + structure.get_coordinates(), + ) + + +def test_nuclear_derivative_factory_registered(): + calculator = algorithms.create("nuclear_derivative_calculator") + + assert isinstance(calculator, algorithms.NuclearDerivativeCalculator) + assert calculator.name() == "finite_difference" + + +def test_qdk_nuclear_derivative_factory_registered(): + calculator = algorithms.create("nuclear_derivative_calculator", "qdk") + + assert isinstance(calculator, algorithms.QdkNuclearDerivativeCalculator) + assert calculator.name() == "qdk" From 3f7fbd4683f7c68b725d1d5d6d1e5c1f7111f769 Mon Sep 17 00:00:00 2001 From: Jan Unsleber Date: Tue, 2 Jun 2026 12:14:53 +0200 Subject: [PATCH 02/13] add qdk path and rs-hybrid support --- .../algorithms/nuclear_derivative.hpp | 2 +- .../qdk/chemistry/data/nuclear_gradients.hpp | 2 +- .../qdk/chemistry/data/nuclear_hessian.hpp | 2 +- .../src/eri/LIBINT2_DIRECT/libint2_direct.cpp | 231 +++++++++++++++++- .../algorithms/nuclear_derivative.cpp | 2 +- .../qdk/chemistry/data/nuclear_gradients.cpp | 2 +- .../qdk/chemistry/data/nuclear_hessian.cpp | 2 +- cpp/tests/test_nuclear_derivative.cpp | 108 ++++++-- .../algorithms/nuclear_derivative.cpp | 2 +- .../src/pybind11/data/nuclear_gradients.cpp | 2 +- python/src/pybind11/data/nuclear_hessian.cpp | 2 +- .../src/qdk_chemistry/algorithms/__init__.py | 2 +- .../algorithms/nuclear_derivative.py | 4 +- 13 files changed, 320 insertions(+), 43 deletions(-) diff --git a/cpp/include/qdk/chemistry/algorithms/nuclear_derivative.hpp b/cpp/include/qdk/chemistry/algorithms/nuclear_derivative.hpp index f13bf8b1e..d9dd8cdc7 100644 --- a/cpp/include/qdk/chemistry/algorithms/nuclear_derivative.hpp +++ b/cpp/include/qdk/chemistry/algorithms/nuclear_derivative.hpp @@ -189,4 +189,4 @@ class QdkNuclearDerivativeCalculator : public NuclearDerivativeCalculator { NuclearDerivativeSeedType seed_or_basis) const override; }; -} // namespace qdk::chemistry::algorithms \ No newline at end of file +} // namespace qdk::chemistry::algorithms diff --git a/cpp/include/qdk/chemistry/data/nuclear_gradients.hpp b/cpp/include/qdk/chemistry/data/nuclear_gradients.hpp index aa635dae9..d90a4ab35 100644 --- a/cpp/include/qdk/chemistry/data/nuclear_gradients.hpp +++ b/cpp/include/qdk/chemistry/data/nuclear_gradients.hpp @@ -139,4 +139,4 @@ class NuclearGradients : public DataClass, static_assert(DataClassCompliant, "NuclearGradients must implement the DataClass interface"); -} // namespace qdk::chemistry::data \ No newline at end of file +} // namespace qdk::chemistry::data diff --git a/cpp/include/qdk/chemistry/data/nuclear_hessian.hpp b/cpp/include/qdk/chemistry/data/nuclear_hessian.hpp index 7939eea09..a7118191b 100644 --- a/cpp/include/qdk/chemistry/data/nuclear_hessian.hpp +++ b/cpp/include/qdk/chemistry/data/nuclear_hessian.hpp @@ -134,4 +134,4 @@ class NuclearHessian : public DataClass, static_assert(DataClassCompliant, "NuclearHessian must implement the DataClass interface"); -} // namespace qdk::chemistry::data \ No newline at end of file +} // namespace qdk::chemistry::data diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/eri/LIBINT2_DIRECT/libint2_direct.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/eri/LIBINT2_DIRECT/libint2_direct.cpp index 10c2c71f9..79ec3ef2f 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/eri/LIBINT2_DIRECT/libint2_direct.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/eri/LIBINT2_DIRECT/libint2_direct.cpp @@ -170,9 +170,11 @@ class ERI { bool use_thread_local_buffers_; ///< Use thread-local buffers (true) or ///< atomic ops (false) double eri_threshold_; ///< Integral screening threshold + size_t n_atoms_; ///< Number of atoms in the molecule ::libint2::BasisSet obs_; ///< Libint2 orbital basis set representation std::vector shell2bf_; ///< Mapping from shell index to first atomic orbital + std::vector sh2atom_; ///< Mapping from shell index to atom index shellpair_list_t splist_; ///< Pre-computed shell pair list for screening shellpair_data_t spdata_; ///< Shell pair geometric data and bounds RowMajorMatrix K_schwarz_; ///< Schwarz screening matrix for integral bounds @@ -203,10 +205,15 @@ class ERI { : spin_density_factor_(spin_density_factor), use_thread_local_buffers_(!use_atomics), eri_threshold_(eri_threshold), + n_atoms_(basis_set.mol->n_atoms), obs_(libint2_util::convert_to_libint_basisset(basis_set)) { QDK_LOG_TRACE_ENTERING(); shell2bf_ = obs_.shell2bf(); + sh2atom_.reserve(basis_set.shells.size()); + for (const auto& shell : basis_set.shells) { + sh2atom_.push_back(static_cast(shell.atom_index)); + } // Compute Shell Pairs std::tie(splist_, spdata_) = compute_shellpairs(obs_, shell_pair_threshold); @@ -256,8 +263,8 @@ class ERI { const size_t mat_size = num_density_matrices * num_atomic_orbitals * num_atomic_orbitals; const bool is_rsx = std::abs(omega) > 1e-12; - - if (is_rsx) throw std::runtime_error("RSX + LIBINT2_DIRECT NYI"); + const double exchange_scale = alpha + beta; + const bool need_erf_exchange = K && is_rsx && std::abs(beta) > 1e-12; // Compute shell block norm of P const auto P_shnrm = compute_shellblock_norm(obs_, P, num_atomic_orbitals); @@ -291,6 +298,17 @@ class ERI { engines_coulomb[0].set_precision(engine_precision); for (int i = 1; i < nthreads; ++i) engines_coulomb[i] = engines_coulomb[0]; + std::vector<::libint2::Engine> engines_erf; + if (need_erf_exchange) { + engines_erf.resize(nthreads); + engines_erf[0] = ::libint2::Engine(::libint2::Operator::erf_coulomb, + obs_.max_nprim(), obs_.max_l(), 0); + engines_erf[0].set_params(omega); + engines_erf[0].set(::libint2::ScreeningMethod::Original); + engines_erf[0].set_precision(engine_precision); + for (int i = 1; i < nthreads; ++i) engines_erf[i] = engines_erf[0]; + } + if (J) std::memset(J, 0, mat_size * sizeof(double)); if (K) std::memset(K, 0, mat_size * sizeof(double)); @@ -324,6 +342,7 @@ class ERI { #endif auto& engine = engines_coulomb[thread_id]; const auto& buf = engine.results(); + auto* engine_erf = need_erf_exchange ? &engines_erf[thread_id] : nullptr; // Get pointers to thread-local buffers double* J_thread = nullptr; @@ -395,11 +414,19 @@ class ERI { obs_[s1], obs_[s2], obs_[s3], obs_[s4], sp12_data, sp34_data); // Coarse integral screening - const auto buf_1234 = buf[0]; - if (buf_1234 == nullptr) continue; + const auto* buf_1234 = buf[0]; + const double* buf_erf_1234 = nullptr; + if (need_erf_exchange) { + engine_erf->compute2<::libint2::Operator::erf_coulomb, + ::libint2::BraKet::xx_xx, 0>( + obs_[s1], obs_[s2], obs_[s3], obs_[s4], sp12_data, + sp34_data); + buf_erf_1234 = engine_erf->results()[0]; + } + if (buf_1234 == nullptr && buf_erf_1234 == nullptr) continue; // Contract shell quartet (J) - if (J) + if (J && buf_1234) for (size_t idm = 0; idm < num_density_matrices; ++idm) { auto* J_cur = use_thread_local_buffers_ @@ -420,7 +447,13 @@ class ERI { for (size_t l = 0; l < n4; ++l, ++ijkl) { const size_t bf4 = bf4_st + l; - const auto value = buf_1234[ijkl] * s1234_deg; + const auto coulomb_value = + buf_1234 ? buf_1234[ijkl] : 0.0; + const auto erf_value = + buf_erf_1234 ? buf_erf_1234[ijkl] : 0.0; + const auto value = (exchange_scale * coulomb_value - + beta * erf_value) * + s1234_deg; // J contractions J_ij += @@ -566,7 +599,8 @@ class ERI { j * num_atomic_orbitals + i] = J_ij; } - // Symmetrize K + scale by alpha/beta + // Symmetrize K. Exchange scaling is applied during contraction so that + // range-separated exchange can combine Coulomb and erf operators. if (K) for (size_t idm = 0; idm < num_density_matrices; ++idm) for (size_t i = 0; i < num_atomic_orbitals; ++i) @@ -575,7 +609,7 @@ class ERI { i * num_atomic_orbitals + j]; auto K_ji = K[idm * num_atomic_orbitals * num_atomic_orbitals + j * num_atomic_orbitals + i]; - K_ij = (alpha + beta) * 0.5 * (K_ij + K_ji); + K_ij = 0.5 * (K_ij + K_ji); K[idm * num_atomic_orbitals * num_atomic_orbitals + i * num_atomic_orbitals + j] = K_ij; K[idm * num_atomic_orbitals * num_atomic_orbitals + @@ -599,15 +633,190 @@ class ERI { * @param beta Scaling factor for DFT exchange * @param omega Range-separation parameter * - * @throws std::runtime_error Always - energy gradients not yet implemented - * * @note These are energy derivatives, not matrix element derivatives */ void get_gradients(const double* P, double* dJ, double* dK, double alpha, double beta, double omega) { QDK_LOG_TRACE_ENTERING(); - throw std::runtime_error("LIBINT2_DIRECT + Gradients Not Yet Implemented"); + AutoTimer t("ERI::get_gradients"); + const size_t num_atomic_orbitals = obs_.nbf(); + const size_t num_atomic_orbitals2 = + num_atomic_orbitals * num_atomic_orbitals; + const size_t nsh = obs_.size(); + const bool is_rsx = std::abs(omega) > 1e-12; + const double exchange_scale = alpha + beta; + const bool need_coulomb_exchange = + dK != nullptr && std::abs(exchange_scale) > 1e-12; + const bool need_erf_exchange = + dK != nullptr && is_rsx && std::abs(beta) > 1e-12; + const bool need_coulomb = dJ != nullptr || need_coulomb_exchange; + const bool need_J = dJ != nullptr; + const bool need_K = need_coulomb_exchange || need_erf_exchange; + + std::vector dJ_buffer; + std::vector dK_buffer; + auto* dJ_out = dJ; + auto* dK_out = dK; + if (!dJ_out) { + dJ_buffer.resize(3 * n_atoms_, 0.0); + dJ_out = dJ_buffer.data(); + } + if (!dK_out) { + dK_buffer.resize(3 * n_atoms_, 0.0); + dK_out = dK_buffer.data(); + } + + std::fill_n(dJ_out, 3 * n_atoms_, 0.0); + std::fill_n(dK_out, 3 * n_atoms_, 0.0); + if (!need_J && !need_K) return; + + std::vector P_total; + const double* P_J = P; + if (spin_density_factor_ > 1) { + P_total.resize(num_atomic_orbitals2, 0.0); + for (size_t idx = 0; idx < num_atomic_orbitals2; ++idx) { + for (size_t idm = 0; idm < spin_density_factor_; ++idm) { + P_total[idx] += P[idm * num_atomic_orbitals2 + idx]; + } + } + P_J = P_total.data(); + } + +#ifdef _OPENMP + const int nthreads = omp_get_max_threads(); +#else + const int nthreads = 1; +#endif + std::vector<::libint2::Engine> engines( + nthreads, ::libint2::Engine(::libint2::Operator::coulomb, + obs_.max_nprim(), obs_.max_l(), 1)); + for (auto& engine : engines) { + engine.set(::libint2::BraKet::xx_xx); + engine.set(::libint2::ScreeningMethod::Original); + engine.set_precision(max_engine_precision); + } + + std::vector<::libint2::Engine> engines_erf; + if (need_erf_exchange) { + engines_erf.resize(nthreads); + engines_erf[0] = ::libint2::Engine(::libint2::Operator::erf_coulomb, + obs_.max_nprim(), obs_.max_l(), 1); + engines_erf[0].set_params(omega); + engines_erf[0].set(::libint2::BraKet::xx_xx); + engines_erf[0].set(::libint2::ScreeningMethod::Original); + engines_erf[0].set_precision(max_engine_precision); + for (int i = 1; i < nthreads; ++i) engines_erf[i] = engines_erf[0]; + } + +#ifdef _OPENMP +#pragma omp parallel reduction(+ : dJ_out[ : 3 * n_atoms_], \ + dK_out[ : 3 * n_atoms_]) +#endif + { +#ifdef _OPENMP + const int thread_id = omp_get_thread_num(); +#else + const int thread_id = 0; +#endif + auto& engine = engines[thread_id]; + const auto& buf = engine.results(); + auto* engine_erf = need_erf_exchange ? &engines_erf[thread_id] : nullptr; + + for (size_t s1 = 0, s1234 = 0; s1 < nsh; ++s1) { + const auto bf1_st = shell2bf_[s1]; + const auto n1 = obs_[s1].size(); + for (size_t s2 = 0; s2 < nsh; ++s2) { + const auto bf2_st = shell2bf_[s2]; + const auto n2 = obs_[s2].size(); + for (size_t s3 = 0; s3 < nsh; ++s3) { + const auto bf3_st = shell2bf_[s3]; + const auto n3 = obs_[s3].size(); + for (size_t s4 = 0; s4 < nsh; ++s4, ++s1234) { + if (s1234 % nthreads != static_cast(thread_id)) continue; + + if (K_schwarz_(s1, s2) * K_schwarz_(s3, s4) < eri_threshold_) + continue; + + const auto bf4_st = shell2bf_[s4]; + const auto n4 = obs_[s4].size(); + + if (need_coulomb) { + engine.compute2<::libint2::Operator::coulomb, + ::libint2::BraKet::xx_xx, 1>( + obs_[s1], obs_[s2], obs_[s3], obs_[s4]); + } + if (need_erf_exchange) { + engine_erf->compute2<::libint2::Operator::erf_coulomb, + ::libint2::BraKet::xx_xx, 1>( + obs_[s1], obs_[s2], obs_[s3], obs_[s4]); + } + + const bool has_coulomb = need_coulomb && buf[0] != nullptr; + const bool has_erf = + need_erf_exchange && engine_erf->results()[0] != nullptr; + if (!has_coulomb && !has_erf) continue; + + const auto& buf_erf = + need_erf_exchange ? engine_erf->results() : engine.results(); + + for (int d = 0; d < 12; ++d) { + const int center = d / 3; + const int xyz = d % 3; + const int atom = center == 0 ? sh2atom_[s1] + : center == 1 ? sh2atom_[s2] + : center == 2 ? sh2atom_[s3] + : sh2atom_[s4]; + const size_t coord = static_cast(atom) + xyz * n_atoms_; + const auto* shset = has_coulomb ? buf[d] : nullptr; + const auto* shset_erf = has_erf ? buf_erf[d] : nullptr; + if (shset == nullptr && shset_erf == nullptr) continue; + + double dJ_coord = 0.0; + double dK_coord = 0.0; + for (size_t i = 0, ijkl = 0; i < n1; ++i) { + const size_t bf1 = bf1_st + i; + for (size_t j = 0; j < n2; ++j) { + const size_t bf2 = bf2_st + j; + for (size_t k = 0; k < n3; ++k) { + const size_t bf3 = bf3_st + k; + for (size_t l = 0; l < n4; ++l, ++ijkl) { + const size_t bf4 = bf4_st + l; + const double coulomb_value = shset ? shset[ijkl] : 0.0; + const double erf_value = + shset_erf ? shset_erf[ijkl] : 0.0; + if (need_J && shset) { + dJ_coord += P_J[bf1 * num_atomic_orbitals + bf2] * + P_J[bf3 * num_atomic_orbitals + bf4] * + coulomb_value; + } + if (need_K) { + for (size_t idm = 0; idm < spin_density_factor_; + ++idm) { + const auto* P_cur = P + idm * num_atomic_orbitals2; + const double scaled_value = + exchange_scale * coulomb_value - + beta * erf_value; + dK_coord += P_cur[bf1 * num_atomic_orbitals + bf3] * + P_cur[bf2 * num_atomic_orbitals + bf4] * + scaled_value; + } + } + } + } + } + } + + if (need_J) dJ_out[coord] -= 0.25 * dJ_coord; + if (need_K) { + dK_out[coord] += 0.5 * dK_coord; + } + } + } + } + } + } + } } /** diff --git a/cpp/src/qdk/chemistry/algorithms/nuclear_derivative.cpp b/cpp/src/qdk/chemistry/algorithms/nuclear_derivative.cpp index 9947a7ca8..458271dcd 100644 --- a/cpp/src/qdk/chemistry/algorithms/nuclear_derivative.cpp +++ b/cpp/src/qdk/chemistry/algorithms/nuclear_derivative.cpp @@ -325,4 +325,4 @@ void NuclearDerivativeCalculatorFactory::register_default_instances() { &make_qdk_nuclear_derivative_calculator); } -} // namespace qdk::chemistry::algorithms \ No newline at end of file +} // namespace qdk::chemistry::algorithms diff --git a/cpp/src/qdk/chemistry/data/nuclear_gradients.cpp b/cpp/src/qdk/chemistry/data/nuclear_gradients.cpp index 33c7cbc5a..8d9aeeac9 100644 --- a/cpp/src/qdk/chemistry/data/nuclear_gradients.cpp +++ b/cpp/src/qdk/chemistry/data/nuclear_gradients.cpp @@ -219,4 +219,4 @@ std::shared_ptr NuclearGradients::_from_hdf5_file( return from_hdf5(file); } -} // namespace qdk::chemistry::data \ No newline at end of file +} // namespace qdk::chemistry::data diff --git a/cpp/src/qdk/chemistry/data/nuclear_hessian.cpp b/cpp/src/qdk/chemistry/data/nuclear_hessian.cpp index b8ed4d606..d7d7c0902 100644 --- a/cpp/src/qdk/chemistry/data/nuclear_hessian.cpp +++ b/cpp/src/qdk/chemistry/data/nuclear_hessian.cpp @@ -209,4 +209,4 @@ std::shared_ptr NuclearHessian::_from_hdf5_file( return from_hdf5(file); } -} // namespace qdk::chemistry::data \ No newline at end of file +} // namespace qdk::chemistry::data diff --git a/cpp/tests/test_nuclear_derivative.cpp b/cpp/tests/test_nuclear_derivative.cpp index f2f05cff0..499db20d4 100644 --- a/cpp/tests/test_nuclear_derivative.cpp +++ b/cpp/tests/test_nuclear_derivative.cpp @@ -32,6 +32,70 @@ void expect_same_structure(const std::shared_ptr& left, EXPECT_EQ(left->get_elements(), right->get_elements()); } +void expect_qdk_analytic_gradient_runs_for_functional( + const std::string& functional) { + auto calculator = NuclearDerivativeCalculatorFactory::create("qdk"); + AlgorithmRef scf_ref("scf_solver", "qdk"); + scf_ref.get_settings()->set("method", functional); + calculator->settings().set("energy_calculator", scf_ref); + + auto [energy, gradients, hessian, wavefunction] = + calculator->run(create_h2_structure(), 0, 1, std::string("sto-3g")); + + ASSERT_TRUE(std::isfinite(energy)); + EXPECT_LT(energy, 0.0); + ASSERT_NE(gradients, nullptr); + EXPECT_EQ(gradients->get_num_atoms(), 2); + EXPECT_EQ(gradients->get_values().size(), 6); + EXPECT_TRUE(gradients->get_values().allFinite()); + + auto gradient_matrix = gradients->as_matrix(); + EXPECT_NEAR(gradient_matrix.col(0).sum(), 0.0, 1.0e-8); + EXPECT_NEAR(gradient_matrix.col(1).sum(), 0.0, 1.0e-8); + EXPECT_NEAR(gradient_matrix.col(2).sum(), 0.0, 1.0e-8); + + EXPECT_FALSE(hessian.has_value()); + ASSERT_TRUE(wavefunction.has_value()); + EXPECT_NE(*wavefunction, nullptr); +} + +NuclearDerivativeResult run_qdk_derivative_for_functional( + const std::string& calculator_name, const std::string& functional) { + auto calculator = NuclearDerivativeCalculatorFactory::create(calculator_name); + AlgorithmRef scf_ref("scf_solver", "qdk"); + scf_ref.get_settings()->set("method", functional); + calculator->settings().set("energy_calculator", scf_ref); + calculator->settings().set("finite_difference_step", 1.0e-2); + + return calculator->run(create_h2_structure(), 0, 1, std::string("sto-3g")); +} + +void expect_qdk_analytic_gradient_matches_numeric( + const std::string& functional) { + auto [numeric_energy, numeric_gradients, numeric_hessian, + numeric_wavefunction] = + run_qdk_derivative_for_functional("finite_difference", functional); + auto [analytic_energy, analytic_gradients, analytic_hessian, + analytic_wavefunction] = + run_qdk_derivative_for_functional("qdk", functional); + + EXPECT_NEAR(analytic_energy, numeric_energy, 1.0e-8); + ASSERT_NE(numeric_gradients, nullptr); + ASSERT_NE(analytic_gradients, nullptr); + EXPECT_FALSE(numeric_hessian.has_value()); + EXPECT_FALSE(analytic_hessian.has_value()); + ASSERT_TRUE(numeric_wavefunction.has_value()); + ASSERT_TRUE(analytic_wavefunction.has_value()); + + const auto& numeric_values = numeric_gradients->get_values(); + const auto& analytic_values = analytic_gradients->get_values(); + ASSERT_EQ(analytic_values.size(), numeric_values.size()); + for (Eigen::Index i = 0; i < analytic_values.size(); ++i) { + EXPECT_NEAR(analytic_values(i), numeric_values(i), 1.0e-3) + << "gradient component " << i; + } +} + } // namespace TEST(NuclearDerivativeDataTest, GradientsReferenceStructureAndRoundTripJson) { @@ -100,28 +164,32 @@ TEST(NuclearDerivativeCalculatorTest, FiniteDifferenceRunsRealScfForH2) { EXPECT_NE(*wavefunction, nullptr); } -TEST(NuclearDerivativeCalculatorTest, QdkAnalyticGradientRunsRealDftForH2) { - auto calculator = NuclearDerivativeCalculatorFactory::create("qdk"); - AlgorithmRef scf_ref("scf_solver", "qdk"); - scf_ref.get_settings()->set("method", std::string("pbe")); - calculator->settings().set("energy_calculator", scf_ref); +TEST(NuclearDerivativeCalculatorTest, + QdkAnalyticGradientRunsRealNonHybridDftForH2) { + expect_qdk_analytic_gradient_runs_for_functional("pbe"); +} - auto [energy, gradients, hessian, wavefunction] = - calculator->run(create_h2_structure(), 0, 1, std::string("sto-3g")); +TEST(NuclearDerivativeCalculatorTest, + QdkAnalyticGradientRunsRealHybridDftForH2) { + expect_qdk_analytic_gradient_runs_for_functional("b3lyp"); +} - ASSERT_TRUE(std::isfinite(energy)); - EXPECT_LT(energy, 0.0); - ASSERT_NE(gradients, nullptr); - EXPECT_EQ(gradients->get_num_atoms(), 2); - EXPECT_EQ(gradients->get_values().size(), 6); - EXPECT_TRUE(gradients->get_values().allFinite()); +TEST(NuclearDerivativeCalculatorTest, + QdkAnalyticGradientRunsRealRangeSeparatedHybridDftForH2) { + expect_qdk_analytic_gradient_runs_for_functional("wB97x"); +} - auto gradient_matrix = gradients->as_matrix(); - EXPECT_NEAR(gradient_matrix.col(0).sum(), 0.0, 1.0e-8); - EXPECT_NEAR(gradient_matrix.col(1).sum(), 0.0, 1.0e-8); - EXPECT_NEAR(gradient_matrix.col(2).sum(), 0.0, 1.0e-8); +TEST(NuclearDerivativeCalculatorTest, + QdkAnalyticGradientMatchesNumericHybridGgaDftForH2) { + expect_qdk_analytic_gradient_matches_numeric("b3lyp"); +} - EXPECT_FALSE(hessian.has_value()); - ASSERT_TRUE(wavefunction.has_value()); - EXPECT_NE(*wavefunction, nullptr); +TEST(NuclearDerivativeCalculatorTest, + QdkAnalyticGradientMatchesNumericGgaDftForH2) { + expect_qdk_analytic_gradient_matches_numeric("pbe"); +} + +TEST(NuclearDerivativeCalculatorTest, + QdkAnalyticGradientMatchesNumericRangeSeparatedHybridDftForH2) { + expect_qdk_analytic_gradient_matches_numeric("wB97x"); } diff --git a/python/src/pybind11/algorithms/nuclear_derivative.cpp b/python/src/pybind11/algorithms/nuclear_derivative.cpp index 49ac3e7dd..d972a883a 100644 --- a/python/src/pybind11/algorithms/nuclear_derivative.cpp +++ b/python/src/pybind11/algorithms/nuclear_derivative.cpp @@ -115,4 +115,4 @@ Compute nuclear derivatives for a molecular structure. R"(QDK nuclear derivative calculator using analytic internal SCF gradients.)") .def(py::init<>(), R"(Create a QDK analytic nuclear derivative calculator.)"); -} \ No newline at end of file +} diff --git a/python/src/pybind11/data/nuclear_gradients.cpp b/python/src/pybind11/data/nuclear_gradients.cpp index 665dedf61..b99c8e5e9 100644 --- a/python/src/pybind11/data/nuclear_gradients.cpp +++ b/python/src/pybind11/data/nuclear_gradients.cpp @@ -122,4 +122,4 @@ Create nuclear gradients for a structure. .def("__repr__", [](const NuclearGradients& gradients) { return gradients.get_summary(); }); -} \ No newline at end of file +} diff --git a/python/src/pybind11/data/nuclear_hessian.cpp b/python/src/pybind11/data/nuclear_hessian.cpp index d4cd2dbea..9dba729e9 100644 --- a/python/src/pybind11/data/nuclear_hessian.cpp +++ b/python/src/pybind11/data/nuclear_hessian.cpp @@ -117,4 +117,4 @@ Create a nuclear Hessian for a structure. py::arg("filename"), R"(Load a Hessian from an HDF5 file.)") .def("__repr__", [](const NuclearHessian& hessian) { return hessian.get_summary(); }); -} \ No newline at end of file +} diff --git a/python/src/qdk_chemistry/algorithms/__init__.py b/python/src/qdk_chemistry/algorithms/__init__.py index cfb0c8d10..6c26ea585 100644 --- a/python/src/qdk_chemistry/algorithms/__init__.py +++ b/python/src/qdk_chemistry/algorithms/__init__.py @@ -85,9 +85,9 @@ "QdkMacisAsci", "QdkMacisCas", "QdkMacisPmc", + "QdkNuclearDerivativeCalculator", "QdkOccupationActiveSpaceSelector", "QdkPipekMezeyLocalizer", - "QdkNuclearDerivativeCalculator", "QdkQubitMapper", "QdkScfSolver", "QdkStabilityChecker", diff --git a/python/src/qdk_chemistry/algorithms/nuclear_derivative.py b/python/src/qdk_chemistry/algorithms/nuclear_derivative.py index 1758d15ea..1dfbf3d34 100644 --- a/python/src/qdk_chemistry/algorithms/nuclear_derivative.py +++ b/python/src/qdk_chemistry/algorithms/nuclear_derivative.py @@ -5,7 +5,7 @@ # Licensed under the MIT License. See LICENSE.txt in the project root for license information. # -------------------------------------------------------------------------------------------- -from qdk_chemistry._core._algorithms import ( # noqa: F401 +from qdk_chemistry._core._algorithms import ( FiniteDifferenceNuclearDerivativeCalculator, NuclearDerivativeCalculator, QdkNuclearDerivativeCalculator, @@ -15,4 +15,4 @@ "FiniteDifferenceNuclearDerivativeCalculator", "NuclearDerivativeCalculator", "QdkNuclearDerivativeCalculator", -] \ No newline at end of file +] From e7d2162f79123286b6a7e4daf88171d50e8eb5e6 Mon Sep 17 00:00:00 2001 From: Jan Unsleber Date: Wed, 3 Jun 2026 11:17:42 +0200 Subject: [PATCH 03/13] fixes --- .../src/eri/LIBINT2_DIRECT/libint2_direct.cpp | 20 +++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/eri/LIBINT2_DIRECT/libint2_direct.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/eri/LIBINT2_DIRECT/libint2_direct.cpp index 79ec3ef2f..db04e5768 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/eri/LIBINT2_DIRECT/libint2_direct.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/eri/LIBINT2_DIRECT/libint2_direct.cpp @@ -447,13 +447,7 @@ class ERI { for (size_t l = 0; l < n4; ++l, ++ijkl) { const size_t bf4 = bf4_st + l; - const auto coulomb_value = - buf_1234 ? buf_1234[ijkl] : 0.0; - const auto erf_value = - buf_erf_1234 ? buf_erf_1234[ijkl] : 0.0; - const auto value = (exchange_scale * coulomb_value - - beta * erf_value) * - s1234_deg; + const auto value = buf_1234[ijkl] * s1234_deg; // J contractions J_ij += @@ -510,7 +504,13 @@ class ERI { for (size_t l = 0; l < n4; ++l, ++ijkl) { const size_t bf4 = bf4_st + l; - const auto value = buf_1234[ijkl] * s1234_deg; + const auto coulomb_value = + buf_1234 ? buf_1234[ijkl] : 0.0; + const auto erf_value = + buf_erf_1234 ? buf_erf_1234[ijkl] : 0.0; + const auto value = (exchange_scale * coulomb_value - + beta * erf_value) * + s1234_deg; // K contractions K_ik += 0.25 * @@ -807,9 +807,9 @@ class ERI { } } - if (need_J) dJ_out[coord] -= 0.25 * dJ_coord; + if (need_J) dJ_out[coord] += 0.5 * dJ_coord; if (need_K) { - dK_out[coord] += 0.5 * dK_coord; + dK_out[coord] -= 0.25 * dK_coord; } } } From 0c0a8e93fcbc597c7899d99953c3c8ed81d0eb0b Mon Sep 17 00:00:00 2001 From: Jan Unsleber Date: Wed, 3 Jun 2026 16:14:28 +0200 Subject: [PATCH 04/13] fix mcscf/casci num-grads and add settings descriptions --- .../algorithms/nuclear_derivative.hpp | 59 ++++- cpp/include/qdk/chemistry/data/settings.hpp | 18 +- .../algorithms/nuclear_derivative.cpp | 215 ++++++++++++++++-- cpp/src/qdk/chemistry/data/settings.cpp | 4 +- cpp/tests/test_nuclear_derivative.cpp | 111 +++++++++ 5 files changed, 384 insertions(+), 23 deletions(-) diff --git a/cpp/include/qdk/chemistry/algorithms/nuclear_derivative.hpp b/cpp/include/qdk/chemistry/algorithms/nuclear_derivative.hpp index d9dd8cdc7..e8aa2646f 100644 --- a/cpp/include/qdk/chemistry/algorithms/nuclear_derivative.hpp +++ b/cpp/include/qdk/chemistry/algorithms/nuclear_derivative.hpp @@ -45,22 +45,65 @@ class NuclearDerivativeSettings : public data::Settings { * @brief Construct settings with finite-difference defaults. */ NuclearDerivativeSettings() { - set_default("energy_calculator", data::AlgorithmRef("scf_solver", "qdk")); - set_default("orbital_solver", data::AlgorithmRef("scf_solver", "qdk")); + set_default( + "energy_calculator", data::AlgorithmRef("scf_solver", "qdk"), + "Algorithm used for each energy evaluation. Use an scf_solver for " + "direct SCF finite differences, a multi_configuration_scf solver for " + "MCSCF energies, or a multi_configuration_calculator such as CASCI or " + "ASCI for Hamiltonian-based multi-reference energies."); + allow_algorithm_ref_type_change("energy_calculator"); + set_default( + "orbital_solver", data::AlgorithmRef("scf_solver", "qdk"), + "SCF solver used to generate reference orbitals for multi-reference " + "energy paths. This setting is ignored for direct SCF energy paths and " + "is skipped when the derivative input seed already provides usable " + "orbitals for the current geometry."); set_default("hamiltonian_constructor", - data::AlgorithmRef("hamiltonian_constructor", "qdk")); - set_default("compute_hessian", false); + data::AlgorithmRef("hamiltonian_constructor", "qdk"), + "Hamiltonian constructor used when energy_calculator is a " + "multi_configuration_calculator. It builds the active-space " + "Hamiltonian from the reference orbitals before the MR energy " + "calculation."); + set_default("compute_hessian", false, + "Whether to compute a nuclear Hessian in addition to energy " + "and gradients. Currently supported by the finite_difference " + "implementation only."); set_default("finite_difference_step", 1.0e-3, - "Nuclear displacement step in Bohr", + "Central finite-difference nuclear displacement step in Bohr. " + "Used by the finite_difference implementation for numeric " + "gradients and Hessians.", data::BoundConstraint{1.0e-8, 1.0}); - set_default("symmetrize_hessian", true); + set_default("symmetrize_hessian", true, + "Whether to replace the finite-difference Hessian by the " + "average of itself and its transpose before returning it."); + set_default( + "reuse_seed_active_space", true, + "For multi-reference energy paths, reuse active and inactive orbital " + "indices from an Orbitals or Wavefunction seed when fresh reference " + "orbitals are generated for another geometry. Orbital coefficients are " + "not reused for displaced finite-difference geometries."); + set_default( + "localize_reference_orbitals", false, + "Whether to localize reference orbitals before multi-reference energy " + "evaluations. Localization is applied only to MR energy paths and uses " + "the current active-space orbital indices as the localization subset."); + set_default( + "orbital_localizer", + data::AlgorithmRef("orbital_localizer", "qdk_pipek_mezey"), + "Orbital localizer used when localize_reference_orbitals is true. The " + "localizer runs after reference orbitals are obtained and active-space " + "metadata from the seed has been reapplied."); set_default( "n_active_alpha_electrons", static_cast(0), - "Active alpha electrons for multi-configuration energy paths", + "Number of alpha electrons in the active space for " + "multi_configuration_scf and multi_configuration_calculator energy " + "paths. Must be positive when those energy paths are selected.", data::BoundConstraint{0, std::numeric_limits::max()}); set_default( "n_active_beta_electrons", static_cast(0), - "Active beta electrons for multi-configuration energy paths", + "Number of beta electrons in the active space for " + "multi_configuration_scf and multi_configuration_calculator energy " + "paths. Must be positive when those energy paths are selected.", data::BoundConstraint{0, std::numeric_limits::max()}); } }; diff --git a/cpp/include/qdk/chemistry/data/settings.hpp b/cpp/include/qdk/chemistry/data/settings.hpp index fbe60adc8..90e583b9c 100644 --- a/cpp/include/qdk/chemistry/data/settings.hpp +++ b/cpp/include/qdk/chemistry/data/settings.hpp @@ -17,6 +17,7 @@ #include #include #include +#include #include #include #include @@ -421,7 +422,8 @@ class Settings : public DataClass, * @throws SettingTypeMismatch if @p value type does not match the existing * type for @p key * @throws std::invalid_argument if an AlgorithmRef's algorithm_type is - * changed, or if a constrained value is out of range / not in the + * changed for a setting that is not marked as a cross-type dispatch + * setting, or if a constrained value is out of range / not in the * allowed set */ void set(const std::string& key, const SettingValue& value); @@ -923,6 +925,19 @@ class Settings : public DataClass, bool show_undocumented = false) const; protected: + /** + * @brief Allow an AlgorithmRef setting to change algorithm_type. + * + * Most nested algorithm settings are intentionally fixed to one factory type. + * Some settings are dispatch points across factory types and can opt in to + * accepting AlgorithmRef values with different algorithm_type fields. + * + * @param key Existing AlgorithmRef setting key to mark as cross-type. + */ + void allow_algorithm_ref_type_change(const std::string& key) { + algorithm_ref_type_change_allowed_.insert(key); + } + /** * @brief Set a default value for a setting (only if not already set) - * variant version @@ -1159,6 +1174,7 @@ class Settings : public DataClass, std::map descriptions_; std::map limits_; std::map documented_; + std::set algorithm_ref_type_change_allowed_; /// Flag to indicate if settings are locked mutable bool _locked = false; diff --git a/cpp/src/qdk/chemistry/algorithms/nuclear_derivative.cpp b/cpp/src/qdk/chemistry/algorithms/nuclear_derivative.cpp index 458271dcd..ecbfcc020 100644 --- a/cpp/src/qdk/chemistry/algorithms/nuclear_derivative.cpp +++ b/cpp/src/qdk/chemistry/algorithms/nuclear_derivative.cpp @@ -2,11 +2,14 @@ // Licensed under the MIT License. See LICENSE.txt in the project root for // license information. +#include #include +#include #include #include #include #include +#include #include #include @@ -45,6 +48,11 @@ struct EnergyEvaluation { std::optional> wavefunction; }; +struct ReferenceOrbitals { + std::shared_ptr orbitals; + std::shared_ptr wavefunction; +}; + std::shared_ptr copy_structure( const std::shared_ptr& structure) { if (!structure) { @@ -114,6 +122,192 @@ BasisOrGuessType seed_to_scf_input(const NuclearDerivativeSeedType& seed, return basis->get_name(); } +template +typename Factory::return_type create_from_ref(const data::AlgorithmRef& ref); + +unsigned int active_electrons(const data::Settings& settings, + const std::string& key); + +std::shared_ptr seed_to_orbitals( + const NuclearDerivativeSeedType& seed) { + std::shared_ptr orbitals; + if (std::holds_alternative>(seed)) { + orbitals = std::get>(seed); + if (!orbitals) { + throw std::invalid_argument("Orbital seed must not be null"); + } + } else if (std::holds_alternative>( + seed)) { + auto wavefunction = std::get>(seed); + if (!wavefunction) { + throw std::invalid_argument("Wavefunction seed must not be null"); + } + orbitals = wavefunction->get_orbitals(); + } + + if (!orbitals) { + return nullptr; + } + if (!orbitals->get_basis_set()) { + throw std::invalid_argument( + "Orbital or wavefunction seed must include orbitals with a basis set"); + } + return orbitals; +} + +std::shared_ptr seed_to_wavefunction( + const NuclearDerivativeSeedType& seed) { + if (!std::holds_alternative>(seed)) { + return nullptr; + } + auto wavefunction = std::get>(seed); + if (!wavefunction) { + throw std::invalid_argument("Wavefunction seed must not be null"); + } + return wavefunction; +} + +std::shared_ptr copy_active_space_metadata( + const std::shared_ptr& orbitals, + const std::shared_ptr& metadata_source) { + if (!orbitals || !metadata_source) { + return orbitals; + } + + const auto& [active_a, active_b] = + metadata_source->get_active_space_indices(); + const auto& [inactive_a, inactive_b] = + metadata_source->get_inactive_space_indices(); + std::optional ao_overlap; + if (orbitals->has_overlap_matrix()) { + ao_overlap = orbitals->get_overlap_matrix(); + } + std::shared_ptr basis_set = orbitals->get_basis_set(); + + if (orbitals->is_restricted()) { + std::optional energies; + if (orbitals->has_energies()) { + energies = orbitals->get_energies().first; + } + return std::make_shared( + orbitals->get_coefficients().first, energies, ao_overlap, basis_set, + std::make_tuple( + std::vector(active_a.begin(), active_a.end()), + std::vector(inactive_a.begin(), inactive_a.end()))); + } + + std::optional energies_a; + std::optional energies_b; + if (orbitals->has_energies()) { + auto [source_energies_a, source_energies_b] = orbitals->get_energies(); + energies_a = source_energies_a; + energies_b = source_energies_b; + } + return std::make_shared( + orbitals->get_coefficients().first, orbitals->get_coefficients().second, + energies_a, energies_b, ao_overlap, basis_set, + std::make_tuple( + std::vector(active_a.begin(), active_a.end()), + std::vector(active_b.begin(), active_b.end()), + std::vector(inactive_a.begin(), inactive_a.end()), + std::vector(inactive_b.begin(), inactive_b.end()))); +} + +std::string active_determinant_string(size_t n_active_orbitals, + unsigned int n_alpha, + unsigned int n_beta) { + if (n_alpha > n_active_orbitals || n_beta > n_active_orbitals) { + throw std::invalid_argument( + "Active electron count exceeds the number of active orbitals"); + } + + std::string determinant(n_active_orbitals, '0'); + size_t alpha_remaining = n_alpha; + size_t beta_remaining = n_beta; + for (auto& occupation : determinant) { + if (alpha_remaining > 0 && beta_remaining > 0) { + occupation = '2'; + --alpha_remaining; + --beta_remaining; + } else if (alpha_remaining > 0) { + occupation = 'u'; + --alpha_remaining; + } else if (beta_remaining > 0) { + occupation = 'd'; + --beta_remaining; + } + } + return determinant; +} + +std::shared_ptr wavefunction_from_orbitals( + std::shared_ptr orbitals, unsigned int n_active_alpha, + unsigned int n_active_beta) { + const auto& [active_a, active_b] = orbitals->get_active_space_indices(); + if (active_a.size() != active_b.size()) { + throw std::invalid_argument( + "Reference orbital localization requires matching alpha and beta " + "active spaces"); + } + auto determinant = data::Configuration(active_determinant_string( + active_a.size(), n_active_alpha, n_active_beta)); + auto container = + std::make_unique(determinant, orbitals); + return std::make_shared(std::move(container)); +} + +ReferenceOrbitals localize_reference_orbitals(const data::Settings& settings, + ReferenceOrbitals reference) { + if (!settings.get("localize_reference_orbitals")) { + return reference; + } + + if (!reference.wavefunction) { + reference.wavefunction = wavefunction_from_orbitals( + reference.orbitals, + active_electrons(settings, "n_active_alpha_electrons"), + active_electrons(settings, "n_active_beta_electrons")); + } + + auto localizer = create_from_ref( + settings.get("orbital_localizer")); + auto [loc_indices_a, loc_indices_b] = + reference.orbitals->get_active_space_indices(); + reference.wavefunction = + localizer->run(reference.wavefunction, loc_indices_a, loc_indices_b); + reference.orbitals = reference.wavefunction->get_orbitals(); + return reference; +} + +ReferenceOrbitals reference_orbitals_for_mr_energy( + const data::Settings& settings, std::shared_ptr structure, + int charge, int spin_multiplicity, const NuclearDerivativeSeedType& seed, + bool allow_orbital_seed) { + auto seed_orbitals = seed_to_orbitals(seed); + if (allow_orbital_seed && seed_orbitals) { + return localize_reference_orbitals( + settings, {seed_orbitals, seed_to_wavefunction(seed)}); + } + + auto orbital_solver = create_from_ref( + settings.get("orbital_solver")); + auto [_, reference_wavefunction] = + orbital_solver->run(structure, charge, spin_multiplicity, + seed_to_scf_input(seed, allow_orbital_seed)); + auto reference_orbitals = reference_wavefunction->get_orbitals(); + if (settings.get("reuse_seed_active_space")) { + reference_orbitals = + copy_active_space_metadata(reference_orbitals, seed_orbitals); + if (reference_orbitals != reference_wavefunction->get_orbitals()) { + reference_wavefunction = + detail::new_wavefunction(reference_wavefunction, reference_orbitals); + } + } + + return localize_reference_orbitals( + settings, {reference_orbitals, reference_wavefunction}); +} + template typename Factory::return_type create_from_ref(const data::AlgorithmRef& ref) { auto instance = Factory::create(ref.get_algorithm_name()); @@ -150,14 +344,12 @@ EnergyEvaluation evaluate_energy(const data::Settings& settings, } if (algorithm_type == MultiConfigurationScfFactory::algorithm_type_name()) { - auto orbital_solver = create_from_ref( - settings.get("orbital_solver")); - auto [_, reference_wavefunction] = - orbital_solver->run(structure, charge, spin_multiplicity, - seed_to_scf_input(seed, allow_wavefunction_seed)); + auto reference = reference_orbitals_for_mr_energy( + settings, structure, charge, spin_multiplicity, seed, + allow_wavefunction_seed); auto calculator = create_from_ref(ref); auto [energy, wavefunction] = - calculator->run(reference_wavefunction->get_orbitals(), + calculator->run(reference.orbitals, active_electrons(settings, "n_active_alpha_electrons"), active_electrons(settings, "n_active_beta_electrons")); return {energy, wavefunction}; @@ -165,16 +357,13 @@ EnergyEvaluation evaluate_energy(const data::Settings& settings, if (algorithm_type == MultiConfigurationCalculatorFactory::algorithm_type_name()) { - auto orbital_solver = create_from_ref( - settings.get("orbital_solver")); - auto [_, reference_wavefunction] = - orbital_solver->run(structure, charge, spin_multiplicity, - seed_to_scf_input(seed, allow_wavefunction_seed)); + auto reference = reference_orbitals_for_mr_energy( + settings, structure, charge, spin_multiplicity, seed, + allow_wavefunction_seed); auto hamiltonian_constructor = create_from_ref( settings.get("hamiltonian_constructor")); - auto hamiltonian = - hamiltonian_constructor->run(reference_wavefunction->get_orbitals()); + auto hamiltonian = hamiltonian_constructor->run(reference.orbitals); auto calculator = create_from_ref(ref); auto [energy, wavefunction] = calculator->run( hamiltonian, active_electrons(settings, "n_active_alpha_electrons"), diff --git a/cpp/src/qdk/chemistry/data/settings.cpp b/cpp/src/qdk/chemistry/data/settings.cpp index 2ab4f191b..25ee47459 100644 --- a/cpp/src/qdk/chemistry/data/settings.cpp +++ b/cpp/src/qdk/chemistry/data/settings.cpp @@ -94,7 +94,9 @@ void Settings::set(const std::string& key, const SettingValue& value) { if (std::holds_alternative(value)) { const auto& existing = std::get(settings_[key]); const auto& incoming = std::get(value); - if (incoming.get_algorithm_type() != existing.get_algorithm_type()) { + if (incoming.get_algorithm_type() != existing.get_algorithm_type() && + algorithm_ref_type_change_allowed_.find(key) == + algorithm_ref_type_change_allowed_.end()) { throw std::invalid_argument( "Algorithm type for setting '" + key + "' is fixed to '" + existing.get_algorithm_type() + "' and cannot be changed to '" + diff --git a/cpp/tests/test_nuclear_derivative.cpp b/cpp/tests/test_nuclear_derivative.cpp index 499db20d4..fd70ffcf3 100644 --- a/cpp/tests/test_nuclear_derivative.cpp +++ b/cpp/tests/test_nuclear_derivative.cpp @@ -6,16 +6,68 @@ #include #include +#include #include +#include #include #include +#include #include +#include "ut_common.hpp" + using namespace qdk::chemistry::algorithms; using namespace qdk::chemistry::data; namespace { +std::vector> recorded_active_spaces; +std::vector> recorded_inactive_spaces; + +std::string determinant_string(size_t n_active_orbitals, unsigned int n_alpha, + unsigned int n_beta) { + std::string determinant(n_active_orbitals, '0'); + for (auto& occupation : determinant) { + if (n_alpha > 0 && n_beta > 0) { + occupation = '2'; + --n_alpha; + --n_beta; + } else if (n_alpha > 0) { + occupation = 'u'; + --n_alpha; + } else if (n_beta > 0) { + occupation = 'd'; + --n_beta; + } + } + return determinant; +} + +class RecordingMultiConfigurationCalculator + : public MultiConfigurationCalculator { + public: + std::string name() const override { + return "_test_nuclear_derivative_recording_mc"; + } + + protected: + std::pair> _run_impl( + std::shared_ptr hamiltonian, unsigned int n_alpha, + unsigned int n_beta) const override { + auto orbitals = hamiltonian->get_orbitals(); + recorded_active_spaces.push_back( + orbitals->get_active_space_indices().first); + recorded_inactive_spaces.push_back( + orbitals->get_inactive_space_indices().first); + + auto determinant = Configuration(determinant_string( + orbitals->get_active_space_indices().first.size(), n_alpha, n_beta)); + auto container = + std::make_unique(determinant, orbitals); + return {0.0, std::make_shared(std::move(container))}; + } +}; + std::shared_ptr create_h2_structure() { std::vector coordinates = {{0.0, 0.0, 0.0}, {0.0, 0.0, 1.4}}; std::vector elements = {Element::H, Element::H}; @@ -98,6 +150,19 @@ void expect_qdk_analytic_gradient_matches_numeric( } // namespace +TEST(NuclearDerivativeSettingsTest, SettingsHaveUserFacingDescriptions) { + NuclearDerivativeSettings settings; + for (const auto& key : + {"energy_calculator", "orbital_solver", "hamiltonian_constructor", + "compute_hessian", "finite_difference_step", "symmetrize_hessian", + "reuse_seed_active_space", "localize_reference_orbitals", + "orbital_localizer", "n_active_alpha_electrons", + "n_active_beta_electrons"}) { + EXPECT_TRUE(settings.has_description(key)) << key; + EXPECT_FALSE(settings.get_description(key).empty()) << key; + } +} + TEST(NuclearDerivativeDataTest, GradientsReferenceStructureAndRoundTripJson) { auto structure = create_h2_structure(); Eigen::VectorXd values(6); @@ -164,6 +229,52 @@ TEST(NuclearDerivativeCalculatorTest, FiniteDifferenceRunsRealScfForH2) { EXPECT_NE(*wavefunction, nullptr); } +TEST(NuclearDerivativeCalculatorTest, + MultiReferenceFiniteDifferenceReusesSeedActiveSpace) { + recorded_active_spaces.clear(); + recorded_inactive_spaces.clear(); + MultiConfigurationCalculatorFactory::unregister_instance( + "_test_nuclear_derivative_recording_mc"); + MultiConfigurationCalculatorFactory::register_instance( + []() -> MultiConfigurationCalculatorFactory::return_type { + return std::make_unique(); + }); + + auto structure = create_h2_structure(); + auto scf_solver = ScfSolverFactory::create(); + auto [_, wavefunction] = scf_solver->run(structure, 0, 1, "sto-3g"); + auto seed_orbitals = + testing::with_active_space(wavefunction->get_orbitals(), + std::vector{0}, std::vector{}); + + auto calculator = NuclearDerivativeCalculatorFactory::create(); + calculator->settings().set( + "energy_calculator", + AlgorithmRef("multi_configuration_calculator", + "_test_nuclear_derivative_recording_mc")); + calculator->settings().set("n_active_alpha_electrons", 1); + calculator->settings().set("n_active_beta_electrons", 1); + calculator->settings().set("finite_difference_step", 1.0e-2); + + auto [energy, gradients, hessian, result_wavefunction] = + calculator->run(structure, 0, 1, seed_orbitals); + + EXPECT_TRUE(std::isfinite(energy)); + ASSERT_NE(gradients, nullptr); + EXPECT_FALSE(hessian.has_value()); + ASSERT_TRUE(result_wavefunction.has_value()); + ASSERT_FALSE(recorded_active_spaces.empty()); + for (const auto& active_space : recorded_active_spaces) { + EXPECT_EQ(active_space, std::vector{0}); + } + for (const auto& inactive_space : recorded_inactive_spaces) { + EXPECT_TRUE(inactive_space.empty()); + } + + MultiConfigurationCalculatorFactory::unregister_instance( + "_test_nuclear_derivative_recording_mc"); +} + TEST(NuclearDerivativeCalculatorTest, QdkAnalyticGradientRunsRealNonHybridDftForH2) { expect_qdk_analytic_gradient_runs_for_functional("pbe"); From 23d5f3dc0ba21d3751c62c2b5eccfc2ceb467cff Mon Sep 17 00:00:00 2001 From: Jan Unsleber Date: Fri, 5 Jun 2026 21:36:05 +0200 Subject: [PATCH 05/13] comments --- .../algorithms/nuclear_derivative.hpp | 47 ++++++-- .../qdk/chemistry/data/nuclear_gradients.hpp | 10 +- .../qdk/chemistry/data/nuclear_hessian.hpp | 5 - .../chemistry/algorithms/microsoft/scf.cpp | 2 +- .../algorithms/nuclear_derivative.cpp | 18 +-- .../qdk/chemistry/data/nuclear_gradients.cpp | 17 ++- .../qdk/chemistry/data/nuclear_hessian.cpp | 8 +- cpp/tests/test_nuclear_derivative.cpp | 109 +++++++++++------- 8 files changed, 129 insertions(+), 87 deletions(-) diff --git a/cpp/include/qdk/chemistry/algorithms/nuclear_derivative.hpp b/cpp/include/qdk/chemistry/algorithms/nuclear_derivative.hpp index e8aa2646f..9238c6b16 100644 --- a/cpp/include/qdk/chemistry/algorithms/nuclear_derivative.hpp +++ b/cpp/include/qdk/chemistry/algorithms/nuclear_derivative.hpp @@ -20,7 +20,11 @@ namespace qdk::chemistry::algorithms { /** - * @brief Initial basis, orbitals, or wavefunction seed for a derivative run. + * @brief Initial basis-set name, basis, orbitals, or wavefunction seed for a + * derivative run. + * + * A string seed is interpreted as the basis-set name, such as "sto-3g", to use + * for SCF calculations. */ using NuclearDerivativeSeedType = std::variant, @@ -42,7 +46,7 @@ using NuclearDerivativeResult = class NuclearDerivativeSettings : public data::Settings { public: /** - * @brief Construct settings with finite-difference defaults. + * @brief Construct settings with shared nuclear derivative defaults. */ NuclearDerivativeSettings() { set_default( @@ -66,16 +70,7 @@ class NuclearDerivativeSettings : public data::Settings { "calculation."); set_default("compute_hessian", false, "Whether to compute a nuclear Hessian in addition to energy " - "and gradients. Currently supported by the finite_difference " - "implementation only."); - set_default("finite_difference_step", 1.0e-3, - "Central finite-difference nuclear displacement step in Bohr. " - "Used by the finite_difference implementation for numeric " - "gradients and Hessians.", - data::BoundConstraint{1.0e-8, 1.0}); - set_default("symmetrize_hessian", true, - "Whether to replace the finite-difference Hessian by the " - "average of itself and its transpose before returning it."); + "and gradients."); set_default( "reuse_seed_active_space", true, "For multi-reference energy paths, reuse active and inactive orbital " @@ -108,6 +103,27 @@ class NuclearDerivativeSettings : public data::Settings { } }; +/** + * @class FiniteDifferenceNuclearDerivativeSettings + * @brief Settings for finite-difference nuclear derivative calculations. + */ +class FiniteDifferenceNuclearDerivativeSettings + : public NuclearDerivativeSettings { + public: + /** + * @brief Construct settings with finite-difference defaults. + */ + FiniteDifferenceNuclearDerivativeSettings() : NuclearDerivativeSettings() { + set_default("finite_difference_step", 1.0e-3, + "Central finite-difference nuclear displacement step in Bohr. " + "Used for numeric gradients and Hessians.", + data::BoundConstraint{1.0e-8, 1.0}); + set_default("symmetrize_hessian", true, + "Whether to replace the finite-difference Hessian by the " + "average of itself and its transpose before returning it."); + } +}; + /** * @class NuclearDerivativeCalculator * @brief Base class for nuclear derivative algorithms. @@ -182,6 +198,13 @@ struct NuclearDerivativeCalculatorFactory class FiniteDifferenceNuclearDerivativeCalculator : public NuclearDerivativeCalculator { public: + /** + * @brief Construct a calculator with finite-difference settings. + */ + FiniteDifferenceNuclearDerivativeCalculator() { + _settings = std::make_unique(); + } + /** * @brief Return the implementation name. */ diff --git a/cpp/include/qdk/chemistry/data/nuclear_gradients.hpp b/cpp/include/qdk/chemistry/data/nuclear_gradients.hpp index d90a4ab35..a6b376eb8 100644 --- a/cpp/include/qdk/chemistry/data/nuclear_gradients.hpp +++ b/cpp/include/qdk/chemistry/data/nuclear_gradients.hpp @@ -28,12 +28,13 @@ class NuclearGradients : public DataClass, * @brief Construct nuclear gradients for a structure. * * @param structure Molecular structure used to compute the gradients. - * @param values Atom-major gradient vector with length 3 * number of atoms. + * @param gradient_values Atom-major gradient vector with length 3 * number of + * atoms. * @throws std::invalid_argument If the structure is null or the vector size * does not match the structure. */ NuclearGradients(std::shared_ptr structure, - const Eigen::VectorXd& values); + const Eigen::VectorXd& gradient_values); /** * @brief Get the molecular structure associated with these gradients. @@ -50,11 +51,6 @@ class NuclearGradients : public DataClass, */ Eigen::MatrixXd as_matrix() const; - /** - * @brief Get the number of atoms in the associated structure. - */ - size_t get_num_atoms() const; - /** * @brief Return the serialized data type name. */ diff --git a/cpp/include/qdk/chemistry/data/nuclear_hessian.hpp b/cpp/include/qdk/chemistry/data/nuclear_hessian.hpp index a7118191b..18bd3460f 100644 --- a/cpp/include/qdk/chemistry/data/nuclear_hessian.hpp +++ b/cpp/include/qdk/chemistry/data/nuclear_hessian.hpp @@ -45,11 +45,6 @@ class NuclearHessian : public DataClass, */ const Eigen::MatrixXd& get_matrix() const { return matrix_; } - /** - * @brief Get the number of atoms in the associated structure. - */ - size_t get_num_atoms() const; - /** * @brief Return the serialized data type name. */ diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf.cpp index e6b7c9bfc..79dcef526 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf.cpp @@ -503,7 +503,7 @@ ScfCalculationResult ScfSolver::_run_with_options( double total_energy = context.result.scf_total_energy; std::optional gradient; - if (require_gradient) { + if (require_gradient && ms_scf_config->mpi.world_rank == 0) { const auto& raw_gradient = context.result.scf_total_gradient; const auto expected_size = 3 * structure->get_num_atoms(); if (raw_gradient.size() != expected_size) { diff --git a/cpp/src/qdk/chemistry/algorithms/nuclear_derivative.cpp b/cpp/src/qdk/chemistry/algorithms/nuclear_derivative.cpp index ecbfcc020..dc2cea53b 100644 --- a/cpp/src/qdk/chemistry/algorithms/nuclear_derivative.cpp +++ b/cpp/src/qdk/chemistry/algorithms/nuclear_derivative.cpp @@ -383,6 +383,9 @@ NuclearDerivativeResult FiniteDifferenceNuclearDerivativeCalculator::_run_impl( std::shared_ptr structure, int charge, int spin_multiplicity, NuclearDerivativeSeedType seed) const { const ScopedLogLevel scoped_log_level(utils::LogLevel::error); + if (!structure) { + throw std::invalid_argument("Structure must not be null"); + } const double step = _settings->get("finite_difference_step"); const bool compute_hessian = _settings->get("compute_hessian"); const auto dimension = @@ -474,6 +477,9 @@ NuclearDerivativeResult QdkNuclearDerivativeCalculator::_run_impl( std::shared_ptr structure, int charge, int spin_multiplicity, NuclearDerivativeSeedType seed) const { const ScopedLogLevel scoped_log_level(utils::LogLevel::error); + if (!structure) { + throw std::invalid_argument("Structure must not be null"); + } if (_settings->get("compute_hessian")) { throw std::invalid_argument( "The QDK analytic nuclear derivative calculator does not currently " @@ -496,15 +502,13 @@ NuclearDerivativeResult QdkNuclearDerivativeCalculator::_run_impl( auto scf_result = solver.run_with_analytic_gradient( structure, charge, spin_multiplicity, seed_to_scf_input(seed, true)); - if (!scf_result.nuclear_gradient.has_value()) { - throw std::runtime_error( - "Internal SCF did not return an analytic nuclear gradient"); + std::shared_ptr gradients; + if (scf_result.nuclear_gradient.has_value()) { + gradients = std::make_shared( + copy_structure(structure), *scf_result.nuclear_gradient); } - return {scf_result.energy, - std::make_shared( - copy_structure(structure), *scf_result.nuclear_gradient), - std::nullopt, scf_result.wavefunction}; + return {scf_result.energy, gradients, std::nullopt, scf_result.wavefunction}; } void NuclearDerivativeCalculatorFactory::register_default_instances() { diff --git a/cpp/src/qdk/chemistry/data/nuclear_gradients.cpp b/cpp/src/qdk/chemistry/data/nuclear_gradients.cpp index 8d9aeeac9..6314fb649 100644 --- a/cpp/src/qdk/chemistry/data/nuclear_gradients.cpp +++ b/cpp/src/qdk/chemistry/data/nuclear_gradients.cpp @@ -18,11 +18,11 @@ namespace qdk::chemistry::data { NuclearGradients::NuclearGradients(std::shared_ptr structure, - const Eigen::VectorXd& values) - : structure_(std::move(structure)), values_(values) { + const Eigen::VectorXd& gradient_values) + : structure_(std::move(structure)), values_(gradient_values) { if (!_is_valid()) { throw std::invalid_argument( - "NuclearGradients require a non-null structure and 3*N values"); + "NuclearGradients requires a non-null structure and 3*N values"); } } @@ -34,8 +34,9 @@ const std::shared_ptr NuclearGradients::get_structure() const { } Eigen::MatrixXd NuclearGradients::as_matrix() const { - Eigen::MatrixXd matrix(get_num_atoms(), 3); - for (Eigen::Index atom = 0; atom < static_cast(get_num_atoms()); + const auto num_atoms = get_structure()->get_num_atoms(); + Eigen::MatrixXd matrix(num_atoms, 3); + for (Eigen::Index atom = 0; atom < static_cast(num_atoms); ++atom) { matrix(atom, 0) = values_(3 * atom); matrix(atom, 1) = values_(3 * atom + 1); @@ -44,13 +45,9 @@ Eigen::MatrixXd NuclearGradients::as_matrix() const { return matrix; } -size_t NuclearGradients::get_num_atoms() const { - return get_structure()->get_num_atoms(); -} - std::string NuclearGradients::get_summary() const { std::ostringstream oss; - oss << "NuclearGradients(" << get_num_atoms() + oss << "NuclearGradients(" << get_structure()->get_num_atoms() << " atoms, norm=" << values_.norm() << ")"; return oss.str(); } diff --git a/cpp/src/qdk/chemistry/data/nuclear_hessian.cpp b/cpp/src/qdk/chemistry/data/nuclear_hessian.cpp index d7d7c0902..998dfa134 100644 --- a/cpp/src/qdk/chemistry/data/nuclear_hessian.cpp +++ b/cpp/src/qdk/chemistry/data/nuclear_hessian.cpp @@ -33,14 +33,10 @@ const std::shared_ptr NuclearHessian::get_structure() const { return structure_; } -size_t NuclearHessian::get_num_atoms() const { - return get_structure()->get_num_atoms(); -} - std::string NuclearHessian::get_summary() const { std::ostringstream oss; - oss << "NuclearHessian(" << get_num_atoms() << " atoms, " << matrix_.rows() - << "x" << matrix_.cols() << ")"; + oss << "NuclearHessian(" << get_structure()->get_num_atoms() << " atoms, " + << matrix_.rows() << "x" << matrix_.cols() << ")"; return oss.str(); } diff --git a/cpp/tests/test_nuclear_derivative.cpp b/cpp/tests/test_nuclear_derivative.cpp index fd70ffcf3..36e27b092 100644 --- a/cpp/tests/test_nuclear_derivative.cpp +++ b/cpp/tests/test_nuclear_derivative.cpp @@ -12,6 +12,7 @@ #include #include #include +#include #include #include "ut_common.hpp" @@ -68,12 +69,6 @@ class RecordingMultiConfigurationCalculator } }; -std::shared_ptr create_h2_structure() { - std::vector coordinates = {{0.0, 0.0, 0.0}, {0.0, 0.0, 1.4}}; - std::vector elements = {Element::H, Element::H}; - return std::make_shared(coordinates, elements); -} - void expect_same_structure(const std::shared_ptr& left, const std::shared_ptr& right) { ASSERT_NE(left, nullptr); @@ -91,14 +86,14 @@ void expect_qdk_analytic_gradient_runs_for_functional( scf_ref.get_settings()->set("method", functional); calculator->settings().set("energy_calculator", scf_ref); - auto [energy, gradients, hessian, wavefunction] = - calculator->run(create_h2_structure(), 0, 1, std::string("sto-3g")); + auto [energy, gradients, hessian, wavefunction] = calculator->run( + testing::create_water_structure(), 0, 1, std::string("sto-3g")); ASSERT_TRUE(std::isfinite(energy)); EXPECT_LT(energy, 0.0); ASSERT_NE(gradients, nullptr); - EXPECT_EQ(gradients->get_num_atoms(), 2); - EXPECT_EQ(gradients->get_values().size(), 6); + EXPECT_EQ(gradients->get_structure()->get_num_atoms(), 3); + EXPECT_EQ(gradients->get_values().size(), 9); EXPECT_TRUE(gradients->get_values().allFinite()); auto gradient_matrix = gradients->as_matrix(); @@ -117,9 +112,12 @@ NuclearDerivativeResult run_qdk_derivative_for_functional( AlgorithmRef scf_ref("scf_solver", "qdk"); scf_ref.get_settings()->set("method", functional); calculator->settings().set("energy_calculator", scf_ref); - calculator->settings().set("finite_difference_step", 1.0e-2); + if (calculator->settings().has("finite_difference_step")) { + calculator->settings().set("finite_difference_step", 1.0e-2); + } - return calculator->run(create_h2_structure(), 0, 1, std::string("sto-3g")); + return calculator->run(testing::create_water_structure(), 0, 1, + std::string("sto-3g")); } void expect_qdk_analytic_gradient_matches_numeric( @@ -154,25 +152,45 @@ TEST(NuclearDerivativeSettingsTest, SettingsHaveUserFacingDescriptions) { NuclearDerivativeSettings settings; for (const auto& key : {"energy_calculator", "orbital_solver", "hamiltonian_constructor", - "compute_hessian", "finite_difference_step", "symmetrize_hessian", - "reuse_seed_active_space", "localize_reference_orbitals", - "orbital_localizer", "n_active_alpha_electrons", - "n_active_beta_electrons"}) { + "compute_hessian", "reuse_seed_active_space", + "localize_reference_orbitals", "orbital_localizer", + "n_active_alpha_electrons", "n_active_beta_electrons"}) { + EXPECT_TRUE(settings.has_description(key)) << key; + EXPECT_FALSE(settings.get_description(key).empty()) << key; + } + EXPECT_FALSE(settings.has("finite_difference_step")); + EXPECT_FALSE(settings.has("symmetrize_hessian")); +} + +TEST(NuclearDerivativeSettingsTest, + FiniteDifferenceSettingsHaveUserFacingDescriptions) { + FiniteDifferenceNuclearDerivativeSettings settings; + for (const auto& key : {"finite_difference_step", "symmetrize_hessian"}) { EXPECT_TRUE(settings.has_description(key)) << key; EXPECT_FALSE(settings.get_description(key).empty()) << key; } } +TEST(NuclearDerivativeSettingsTest, NumericStepBelongsToFiniteDifferenceOnly) { + auto finite_difference = NuclearDerivativeCalculatorFactory::create(); + auto analytic = NuclearDerivativeCalculatorFactory::create("qdk"); + + EXPECT_TRUE(finite_difference->settings().has("finite_difference_step")); + EXPECT_TRUE(finite_difference->settings().has("symmetrize_hessian")); + EXPECT_FALSE(analytic->settings().has("finite_difference_step")); + EXPECT_FALSE(analytic->settings().has("symmetrize_hessian")); +} + TEST(NuclearDerivativeDataTest, GradientsReferenceStructureAndRoundTripJson) { - auto structure = create_h2_structure(); - Eigen::VectorXd values(6); - values << 0.0, 0.1, 0.2, 0.3, 0.4, 0.5; + auto structure = testing::create_water_structure(); + Eigen::VectorXd values(9); + values << 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8; NuclearGradients gradients(structure, values); - EXPECT_EQ(gradients.get_num_atoms(), 2); + EXPECT_EQ(gradients.get_structure()->get_num_atoms(), 3); EXPECT_TRUE(gradients.get_values().isApprox(values)); - EXPECT_EQ(gradients.as_matrix().rows(), 2); + EXPECT_EQ(gradients.as_matrix().rows(), 3); EXPECT_EQ(gradients.as_matrix().cols(), 3); expect_same_structure(structure, gradients.get_structure()); @@ -182,12 +200,12 @@ TEST(NuclearDerivativeDataTest, GradientsReferenceStructureAndRoundTripJson) { } TEST(NuclearDerivativeDataTest, HessianReferencesStructureAndRoundTripsJson) { - auto structure = create_h2_structure(); - Eigen::MatrixXd matrix = Eigen::MatrixXd::Identity(6, 6); + auto structure = testing::create_water_structure(); + Eigen::MatrixXd matrix = Eigen::MatrixXd::Identity(9, 9); NuclearHessian hessian(structure, matrix); - EXPECT_EQ(hessian.get_num_atoms(), 2); + EXPECT_EQ(hessian.get_structure()->get_num_atoms(), 3); EXPECT_TRUE(hessian.get_matrix().isApprox(matrix)); expect_same_structure(structure, hessian.get_structure()); @@ -196,19 +214,19 @@ TEST(NuclearDerivativeDataTest, HessianReferencesStructureAndRoundTripsJson) { expect_same_structure(structure, loaded->get_structure()); } -TEST(NuclearDerivativeCalculatorTest, FiniteDifferenceRunsRealScfForH2) { +TEST(NuclearDerivativeCalculatorTest, FiniteDifferenceRunsRealScfForWater) { auto calculator = NuclearDerivativeCalculatorFactory::create(); calculator->settings().set("compute_hessian", true); calculator->settings().set("finite_difference_step", 1.0e-3); - auto [energy, gradients, hessian, wavefunction] = - calculator->run(create_h2_structure(), 0, 1, std::string("sto-3g")); + auto [energy, gradients, hessian, wavefunction] = calculator->run( + testing::create_water_structure(), 0, 1, std::string("sto-3g")); ASSERT_TRUE(std::isfinite(energy)); EXPECT_LT(energy, 0.0); ASSERT_NE(gradients, nullptr); - EXPECT_EQ(gradients->get_num_atoms(), 2); - EXPECT_EQ(gradients->get_values().size(), 6); + EXPECT_EQ(gradients->get_structure()->get_num_atoms(), 3); + EXPECT_EQ(gradients->get_values().size(), 9); EXPECT_TRUE(gradients->get_values().allFinite()); auto gradient_matrix = gradients->as_matrix(); @@ -218,9 +236,9 @@ TEST(NuclearDerivativeCalculatorTest, FiniteDifferenceRunsRealScfForH2) { ASSERT_TRUE(hessian.has_value()); ASSERT_NE(*hessian, nullptr); - EXPECT_EQ((*hessian)->get_num_atoms(), 2); - EXPECT_EQ((*hessian)->get_matrix().rows(), 6); - EXPECT_EQ((*hessian)->get_matrix().cols(), 6); + EXPECT_EQ((*hessian)->get_structure()->get_num_atoms(), 3); + EXPECT_EQ((*hessian)->get_matrix().rows(), 9); + EXPECT_EQ((*hessian)->get_matrix().cols(), 9); EXPECT_TRUE((*hessian)->get_matrix().allFinite()); EXPECT_TRUE((*hessian)->get_matrix().isApprox( (*hessian)->get_matrix().transpose(), 1.0e-8)); @@ -229,6 +247,19 @@ TEST(NuclearDerivativeCalculatorTest, FiniteDifferenceRunsRealScfForH2) { EXPECT_NE(*wavefunction, nullptr); } +TEST(NuclearDerivativeCalculatorTest, NullStructureThrowsInvalidArgument) { + for (const auto& calculator_name : {"finite_difference", "qdk"}) { + auto calculator = + NuclearDerivativeCalculatorFactory::create(calculator_name); + try { + calculator->run(nullptr, 0, 1, std::string("sto-3g")); + FAIL() << calculator_name << " accepted a null structure"; + } catch (const std::invalid_argument& ex) { + EXPECT_STREQ(ex.what(), "Structure must not be null"); + } + } +} + TEST(NuclearDerivativeCalculatorTest, MultiReferenceFiniteDifferenceReusesSeedActiveSpace) { recorded_active_spaces.clear(); @@ -240,7 +271,7 @@ TEST(NuclearDerivativeCalculatorTest, return std::make_unique(); }); - auto structure = create_h2_structure(); + auto structure = testing::create_water_structure(); auto scf_solver = ScfSolverFactory::create(); auto [_, wavefunction] = scf_solver->run(structure, 0, 1, "sto-3g"); auto seed_orbitals = @@ -276,31 +307,31 @@ TEST(NuclearDerivativeCalculatorTest, } TEST(NuclearDerivativeCalculatorTest, - QdkAnalyticGradientRunsRealNonHybridDftForH2) { + QdkAnalyticGradientRunsRealNonHybridDftForWater) { expect_qdk_analytic_gradient_runs_for_functional("pbe"); } TEST(NuclearDerivativeCalculatorTest, - QdkAnalyticGradientRunsRealHybridDftForH2) { + QdkAnalyticGradientRunsRealHybridDftForWater) { expect_qdk_analytic_gradient_runs_for_functional("b3lyp"); } TEST(NuclearDerivativeCalculatorTest, - QdkAnalyticGradientRunsRealRangeSeparatedHybridDftForH2) { + QdkAnalyticGradientRunsRealRangeSeparatedHybridDftForWater) { expect_qdk_analytic_gradient_runs_for_functional("wB97x"); } TEST(NuclearDerivativeCalculatorTest, - QdkAnalyticGradientMatchesNumericHybridGgaDftForH2) { + QdkAnalyticGradientMatchesNumericHybridGgaDftForWater) { expect_qdk_analytic_gradient_matches_numeric("b3lyp"); } TEST(NuclearDerivativeCalculatorTest, - QdkAnalyticGradientMatchesNumericGgaDftForH2) { + QdkAnalyticGradientMatchesNumericGgaDftForWater) { expect_qdk_analytic_gradient_matches_numeric("pbe"); } TEST(NuclearDerivativeCalculatorTest, - QdkAnalyticGradientMatchesNumericRangeSeparatedHybridDftForH2) { + QdkAnalyticGradientMatchesNumericRangeSeparatedHybridDftForWater) { expect_qdk_analytic_gradient_matches_numeric("wB97x"); } From 6cb999259dad572b56fbe6d09339b1d89be6c62e Mon Sep 17 00:00:00 2001 From: Jan Unsleber Date: Fri, 5 Jun 2026 21:44:47 +0200 Subject: [PATCH 06/13] resolve comments --- .../qdk/chemistry/data/nuclear_hessian.hpp | 4 +- .../algorithms/nuclear_derivative.cpp | 8 +- .../qdk/chemistry/data/nuclear_hessian.cpp | 4 +- cpp/tests/test_nuclear_derivative.cpp | 82 +++++++++++++++++++ 4 files changed, 90 insertions(+), 8 deletions(-) diff --git a/cpp/include/qdk/chemistry/data/nuclear_hessian.hpp b/cpp/include/qdk/chemistry/data/nuclear_hessian.hpp index 18bd3460f..963901d93 100644 --- a/cpp/include/qdk/chemistry/data/nuclear_hessian.hpp +++ b/cpp/include/qdk/chemistry/data/nuclear_hessian.hpp @@ -28,12 +28,12 @@ class NuclearHessian : public DataClass, * @brief Construct a nuclear Hessian for a structure. * * @param structure Molecular structure used to compute the Hessian. - * @param matrix Square 3N by 3N Hessian matrix in Hartree/Bohr^2. + * @param hessian_matrix Square 3N by 3N Hessian matrix in Hartree/Bohr^2. * @throws std::invalid_argument If the structure is null or the matrix shape * does not match the structure. */ NuclearHessian(std::shared_ptr structure, - const Eigen::MatrixXd& matrix); + const Eigen::MatrixXd& hessian_matrix); /** * @brief Get the molecular structure associated with this Hessian. diff --git a/cpp/src/qdk/chemistry/algorithms/nuclear_derivative.cpp b/cpp/src/qdk/chemistry/algorithms/nuclear_derivative.cpp index dc2cea53b..e59c7ea65 100644 --- a/cpp/src/qdk/chemistry/algorithms/nuclear_derivative.cpp +++ b/cpp/src/qdk/chemistry/algorithms/nuclear_derivative.cpp @@ -213,9 +213,9 @@ std::shared_ptr copy_active_space_metadata( std::vector(inactive_b.begin(), inactive_b.end()))); } -std::string active_determinant_string(size_t n_active_orbitals, - unsigned int n_alpha, - unsigned int n_beta) { +std::string hartree_fock_determinant_string(size_t n_active_orbitals, + unsigned int n_alpha, + unsigned int n_beta) { if (n_alpha > n_active_orbitals || n_beta > n_active_orbitals) { throw std::invalid_argument( "Active electron count exceeds the number of active orbitals"); @@ -249,7 +249,7 @@ std::shared_ptr wavefunction_from_orbitals( "Reference orbital localization requires matching alpha and beta " "active spaces"); } - auto determinant = data::Configuration(active_determinant_string( + auto determinant = data::Configuration(hartree_fock_determinant_string( active_a.size(), n_active_alpha, n_active_beta)); auto container = std::make_unique(determinant, orbitals); diff --git a/cpp/src/qdk/chemistry/data/nuclear_hessian.cpp b/cpp/src/qdk/chemistry/data/nuclear_hessian.cpp index 998dfa134..61227e742 100644 --- a/cpp/src/qdk/chemistry/data/nuclear_hessian.cpp +++ b/cpp/src/qdk/chemistry/data/nuclear_hessian.cpp @@ -18,8 +18,8 @@ namespace qdk::chemistry::data { NuclearHessian::NuclearHessian(std::shared_ptr structure, - const Eigen::MatrixXd& matrix) - : structure_(std::move(structure)), matrix_(matrix) { + const Eigen::MatrixXd& hessian_matrix) + : structure_(std::move(structure)), matrix_(hessian_matrix) { if (!_is_valid()) { throw std::invalid_argument( "NuclearHessian requires a non-null structure and a 3*N by 3*N matrix"); diff --git a/cpp/tests/test_nuclear_derivative.cpp b/cpp/tests/test_nuclear_derivative.cpp index 36e27b092..1bc4a05c9 100644 --- a/cpp/tests/test_nuclear_derivative.cpp +++ b/cpp/tests/test_nuclear_derivative.cpp @@ -6,6 +6,7 @@ #include #include +#include #include #include #include @@ -24,6 +25,8 @@ namespace { std::vector> recorded_active_spaces; std::vector> recorded_inactive_spaces; +std::vector> recorded_localized_alpha_spaces; +std::vector> recorded_localized_beta_spaces; std::string determinant_string(size_t n_active_orbitals, unsigned int n_alpha, unsigned int n_beta) { @@ -69,6 +72,23 @@ class RecordingMultiConfigurationCalculator } }; +class RecordingLocalizer : public Localizer { + public: + std::string name() const override { + return "_test_nuclear_derivative_recording_localizer"; + } + + protected: + std::shared_ptr _run_impl( + std::shared_ptr wavefunction, + const std::vector& loc_indices_a, + const std::vector& loc_indices_b) const override { + recorded_localized_alpha_spaces.push_back(loc_indices_a); + recorded_localized_beta_spaces.push_back(loc_indices_b); + return wavefunction; + } +}; + void expect_same_structure(const std::shared_ptr& left, const std::shared_ptr& right) { ASSERT_NE(left, nullptr); @@ -306,6 +326,68 @@ TEST(NuclearDerivativeCalculatorTest, "_test_nuclear_derivative_recording_mc"); } +TEST(NuclearDerivativeCalculatorTest, + MultiReferenceFiniteDifferenceLocalizesReferenceOrbitals) { + recorded_active_spaces.clear(); + recorded_inactive_spaces.clear(); + recorded_localized_alpha_spaces.clear(); + recorded_localized_beta_spaces.clear(); + MultiConfigurationCalculatorFactory::unregister_instance( + "_test_nuclear_derivative_recording_mc"); + MultiConfigurationCalculatorFactory::register_instance( + []() -> MultiConfigurationCalculatorFactory::return_type { + return std::make_unique(); + }); + LocalizerFactory::unregister_instance( + "_test_nuclear_derivative_recording_localizer"); + LocalizerFactory::register_instance([]() -> LocalizerFactory::return_type { + return std::make_unique(); + }); + + auto structure = testing::create_hydrogen_structure(); + auto scf_solver = ScfSolverFactory::create(); + auto [_, wavefunction] = scf_solver->run(structure, 0, 2, "sto-3g"); + auto seed_orbitals = + testing::with_active_space(wavefunction->get_orbitals(), + std::vector{0}, std::vector{}); + + auto calculator = NuclearDerivativeCalculatorFactory::create(); + calculator->settings().set( + "energy_calculator", + AlgorithmRef("multi_configuration_calculator", + "_test_nuclear_derivative_recording_mc")); + calculator->settings().set( + "orbital_localizer", + AlgorithmRef("orbital_localizer", + "_test_nuclear_derivative_recording_localizer")); + calculator->settings().set("localize_reference_orbitals", true); + calculator->settings().set("n_active_alpha_electrons", 1); + calculator->settings().set("n_active_beta_electrons", 1); + calculator->settings().set("finite_difference_step", 1.0e-2); + + auto [energy, gradients, hessian, result_wavefunction] = + calculator->run(structure, 0, 2, seed_orbitals); + + EXPECT_TRUE(std::isfinite(energy)); + ASSERT_NE(gradients, nullptr); + EXPECT_FALSE(hessian.has_value()); + ASSERT_TRUE(result_wavefunction.has_value()); + ASSERT_FALSE(recorded_localized_alpha_spaces.empty()); + ASSERT_EQ(recorded_localized_alpha_spaces.size(), + recorded_localized_beta_spaces.size()); + for (const auto& localized_space : recorded_localized_alpha_spaces) { + EXPECT_EQ(localized_space, std::vector{0}); + } + for (const auto& localized_space : recorded_localized_beta_spaces) { + EXPECT_EQ(localized_space, std::vector{0}); + } + + LocalizerFactory::unregister_instance( + "_test_nuclear_derivative_recording_localizer"); + MultiConfigurationCalculatorFactory::unregister_instance( + "_test_nuclear_derivative_recording_mc"); +} + TEST(NuclearDerivativeCalculatorTest, QdkAnalyticGradientRunsRealNonHybridDftForWater) { expect_qdk_analytic_gradient_runs_for_functional("pbe"); From 13dc463fbe36cbca7253dd714a5fa2d541901d82 Mon Sep 17 00:00:00 2001 From: Jan Unsleber Date: Sat, 6 Jun 2026 17:32:55 +0200 Subject: [PATCH 07/13] cleaning --- cpp/include/qdk/chemistry/utils/logger.hpp | 14 ++++ .../chemistry/algorithms/microsoft/scf.hpp | 10 ++- .../algorithms/nuclear_derivative.cpp | 14 ++-- cpp/src/qdk/chemistry/utils/logger.cpp | 27 +++++--- cpp/tests/test_logger.cpp | 14 ++++ cpp/tests/test_nuclear_derivative.cpp | 66 +++++-------------- cpp/tests/ut_common.hpp | 18 +++++ 7 files changed, 97 insertions(+), 66 deletions(-) diff --git a/cpp/include/qdk/chemistry/utils/logger.hpp b/cpp/include/qdk/chemistry/utils/logger.hpp index b32a16bac..153ed1677 100644 --- a/cpp/include/qdk/chemistry/utils/logger.hpp +++ b/cpp/include/qdk/chemistry/utils/logger.hpp @@ -108,6 +108,20 @@ class Logger { */ static void set_global_level(LogLevel level); + /** + * @brief Restore the global log level if it has not changed. + * + * Atomically compares the current global log level with expected_current. If + * they match, restores restored_level and returns true. Otherwise leaves the + * current level unchanged and returns false. + * + * @param expected_current The level that must still be current + * @param restored_level The level to restore when expected_current matches + * @return Whether restored_level was applied + */ + static bool restore_global_level_if_unchanged(LogLevel expected_current, + LogLevel restored_level); + /** * @brief Get the current global log level * diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf.hpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf.hpp index e70788714..416b159a9 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf.hpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf.hpp @@ -133,10 +133,14 @@ class ScfSolver : public qdk::chemistry::algorithms::ScfSolver { virtual std::string name() const final { return "qdk"; } /** - * @brief Run the internal SCF solver and return analytic nuclear gradients. + * @brief Run the internal SCF solver and return analytic nuclear gradients + * when available. * - * Settings are locked in the same way as the base run() API. The returned - * gradient vector is atom-major with length 3 * number of atoms. + * Settings are locked in the same way as the base run() API. The gradient is + * returned atom-major with length 3 * number of atoms. + * + * @note When MPI is enabled, the analytic gradient is only populated on rank + * 0. See ScfCalculationResult::nuclear_gradient. */ ScfCalculationResult run_with_analytic_gradient( std::shared_ptr structure, int charge, diff --git a/cpp/src/qdk/chemistry/algorithms/nuclear_derivative.cpp b/cpp/src/qdk/chemistry/algorithms/nuclear_derivative.cpp index e59c7ea65..2a38b567c 100644 --- a/cpp/src/qdk/chemistry/algorithms/nuclear_derivative.cpp +++ b/cpp/src/qdk/chemistry/algorithms/nuclear_derivative.cpp @@ -25,13 +25,15 @@ class ScopedLogLevel { : previous_level_(utils::Logger::get_global_level()) { if (static_cast(previous_level_) < static_cast(minimum_level)) { utils::Logger::set_global_level(minimum_level); + scoped_level_ = minimum_level; changed_ = true; } } ~ScopedLogLevel() { if (changed_) { - utils::Logger::set_global_level(previous_level_); + utils::Logger::restore_global_level_if_unchanged(scoped_level_, + previous_level_); } } @@ -40,6 +42,7 @@ class ScopedLogLevel { private: utils::LogLevel previous_level_; + utils::LogLevel scoped_level_ = utils::LogLevel::off; bool changed_ = false; }; @@ -502,11 +505,12 @@ NuclearDerivativeResult QdkNuclearDerivativeCalculator::_run_impl( auto scf_result = solver.run_with_analytic_gradient( structure, charge, spin_multiplicity, seed_to_scf_input(seed, true)); - std::shared_ptr gradients; - if (scf_result.nuclear_gradient.has_value()) { - gradients = std::make_shared( - copy_structure(structure), *scf_result.nuclear_gradient); + if (!scf_result.nuclear_gradient.has_value()) { + throw std::runtime_error( + "Internal SCF did not return the requested analytic nuclear gradient"); } + auto gradients = std::make_shared( + copy_structure(structure), *scf_result.nuclear_gradient); return {scf_result.energy, gradients, std::nullopt, scf_result.wavefunction}; } diff --git a/cpp/src/qdk/chemistry/utils/logger.cpp b/cpp/src/qdk/chemistry/utils/logger.cpp index ec87e7915..63ed87feb 100644 --- a/cpp/src/qdk/chemistry/utils/logger.cpp +++ b/cpp/src/qdk/chemistry/utils/logger.cpp @@ -195,6 +195,12 @@ static void apply_spdlog_global_level_and_flush_policy( spdlog::flush_on(level); } +static void set_global_level_locked(spdlog::level::level_enum level) { + g_global_level = level; + apply_logger_instance_level_and_flush_policy(g_logger, level); + apply_spdlog_global_level_and_flush_policy(level); +} + static void init_global_logger() { try { g_logger = spdlog::stdout_color_mt("qdk-chemistry"); @@ -251,17 +257,22 @@ std::string Logger::get_source_context(const std::source_location& location) { void Logger::set_global_level(LogLevel level) { auto spdlog_level = to_spdlog_level(level); - std::shared_ptr logger; + std::lock_guard lock(g_level_mutex); + set_global_level_locked(spdlog_level); +} - // Update both our tracked level and spdlog's global level - { - std::lock_guard lock(g_level_mutex); - g_global_level = spdlog_level; - logger = g_logger; +bool Logger::restore_global_level_if_unchanged(LogLevel expected_current, + LogLevel restored_level) { + auto expected_spdlog_level = to_spdlog_level(expected_current); + auto restored_spdlog_level = to_spdlog_level(restored_level); + + std::lock_guard lock(g_level_mutex); + if (g_global_level != expected_spdlog_level) { + return false; } - apply_logger_instance_level_and_flush_policy(logger, spdlog_level); - apply_spdlog_global_level_and_flush_policy(spdlog_level); + set_global_level_locked(restored_spdlog_level); + return true; } LogLevel Logger::get_global_level() { diff --git a/cpp/tests/test_logger.cpp b/cpp/tests/test_logger.cpp index a9e854447..28e545b99 100644 --- a/cpp/tests/test_logger.cpp +++ b/cpp/tests/test_logger.cpp @@ -88,6 +88,20 @@ TEST_F(LoggerTest, GlobalLevelControl) { EXPECT_NO_THROW({ QDK_LOGGER().critical("This should be suppressed"); }); } +TEST_F(LoggerTest, ConditionalGlobalLevelRestoreOnlyRestoresExpectedLevel) { + Logger::set_global_level(LogLevel::error); + + EXPECT_TRUE(Logger::restore_global_level_if_unchanged(LogLevel::error, + LogLevel::trace)); + EXPECT_EQ(Logger::get_global_level(), LogLevel::trace); + + Logger::set_global_level(LogLevel::critical); + + EXPECT_FALSE(Logger::restore_global_level_if_unchanged(LogLevel::error, + LogLevel::trace)); + EXPECT_EQ(Logger::get_global_level(), LogLevel::critical); +} + TEST_F(LoggerTest, FormattedLogging) { EXPECT_NO_THROW({ QDK_LOGGER().info("Testing formatted message: value = {}", 42); diff --git a/cpp/tests/test_nuclear_derivative.cpp b/cpp/tests/test_nuclear_derivative.cpp index 1bc4a05c9..082d1dadc 100644 --- a/cpp/tests/test_nuclear_derivative.cpp +++ b/cpp/tests/test_nuclear_derivative.cpp @@ -27,6 +27,7 @@ std::vector> recorded_active_spaces; std::vector> recorded_inactive_spaces; std::vector> recorded_localized_alpha_spaces; std::vector> recorded_localized_beta_spaces; +constexpr double kDftAnalyticNumericGradientTolerance = 5.0e-5; std::string determinant_string(size_t n_active_orbitals, unsigned int n_alpha, unsigned int n_beta) { @@ -99,33 +100,6 @@ void expect_same_structure(const std::shared_ptr& left, EXPECT_EQ(left->get_elements(), right->get_elements()); } -void expect_qdk_analytic_gradient_runs_for_functional( - const std::string& functional) { - auto calculator = NuclearDerivativeCalculatorFactory::create("qdk"); - AlgorithmRef scf_ref("scf_solver", "qdk"); - scf_ref.get_settings()->set("method", functional); - calculator->settings().set("energy_calculator", scf_ref); - - auto [energy, gradients, hessian, wavefunction] = calculator->run( - testing::create_water_structure(), 0, 1, std::string("sto-3g")); - - ASSERT_TRUE(std::isfinite(energy)); - EXPECT_LT(energy, 0.0); - ASSERT_NE(gradients, nullptr); - EXPECT_EQ(gradients->get_structure()->get_num_atoms(), 3); - EXPECT_EQ(gradients->get_values().size(), 9); - EXPECT_TRUE(gradients->get_values().allFinite()); - - auto gradient_matrix = gradients->as_matrix(); - EXPECT_NEAR(gradient_matrix.col(0).sum(), 0.0, 1.0e-8); - EXPECT_NEAR(gradient_matrix.col(1).sum(), 0.0, 1.0e-8); - EXPECT_NEAR(gradient_matrix.col(2).sum(), 0.0, 1.0e-8); - - EXPECT_FALSE(hessian.has_value()); - ASSERT_TRUE(wavefunction.has_value()); - EXPECT_NE(*wavefunction, nullptr); -} - NuclearDerivativeResult run_qdk_derivative_for_functional( const std::string& calculator_name, const std::string& functional) { auto calculator = NuclearDerivativeCalculatorFactory::create(calculator_name); @@ -136,12 +110,13 @@ NuclearDerivativeResult run_qdk_derivative_for_functional( calculator->settings().set("finite_difference_step", 1.0e-2); } - return calculator->run(testing::create_water_structure(), 0, 1, + return calculator->run(testing::create_lithium_hydride_structure(), 0, 1, std::string("sto-3g")); } void expect_qdk_analytic_gradient_matches_numeric( - const std::string& functional) { + const std::string& functional, + double gradient_tolerance = kDftAnalyticNumericGradientTolerance) { auto [numeric_energy, numeric_gradients, numeric_hessian, numeric_wavefunction] = run_qdk_derivative_for_functional("finite_difference", functional); @@ -152,6 +127,10 @@ void expect_qdk_analytic_gradient_matches_numeric( EXPECT_NEAR(analytic_energy, numeric_energy, 1.0e-8); ASSERT_NE(numeric_gradients, nullptr); ASSERT_NE(analytic_gradients, nullptr); + EXPECT_EQ(analytic_gradients->get_structure()->get_num_atoms(), 2); + EXPECT_EQ(analytic_gradients->get_values().size(), 6); + EXPECT_TRUE(analytic_gradients->get_values().allFinite()); + EXPECT_TRUE(numeric_gradients->get_values().allFinite()); EXPECT_FALSE(numeric_hessian.has_value()); EXPECT_FALSE(analytic_hessian.has_value()); ASSERT_TRUE(numeric_wavefunction.has_value()); @@ -161,7 +140,7 @@ void expect_qdk_analytic_gradient_matches_numeric( const auto& analytic_values = analytic_gradients->get_values(); ASSERT_EQ(analytic_values.size(), numeric_values.size()); for (Eigen::Index i = 0; i < analytic_values.size(); ++i) { - EXPECT_NEAR(analytic_values(i), numeric_values(i), 1.0e-3) + EXPECT_NEAR(analytic_values(i), numeric_values(i), gradient_tolerance) << "gradient component " << i; } } @@ -389,31 +368,18 @@ TEST(NuclearDerivativeCalculatorTest, } TEST(NuclearDerivativeCalculatorTest, - QdkAnalyticGradientRunsRealNonHybridDftForWater) { - expect_qdk_analytic_gradient_runs_for_functional("pbe"); -} - -TEST(NuclearDerivativeCalculatorTest, - QdkAnalyticGradientRunsRealHybridDftForWater) { - expect_qdk_analytic_gradient_runs_for_functional("b3lyp"); -} - -TEST(NuclearDerivativeCalculatorTest, - QdkAnalyticGradientRunsRealRangeSeparatedHybridDftForWater) { - expect_qdk_analytic_gradient_runs_for_functional("wB97x"); -} - -TEST(NuclearDerivativeCalculatorTest, - QdkAnalyticGradientMatchesNumericHybridGgaDftForWater) { + QdkAnalyticGradientMatchesNumericHybridGgaDftForLithiumHydride) { expect_qdk_analytic_gradient_matches_numeric("b3lyp"); } TEST(NuclearDerivativeCalculatorTest, - QdkAnalyticGradientMatchesNumericGgaDftForWater) { + QdkAnalyticGradientMatchesNumericGgaDftForLithiumHydride) { expect_qdk_analytic_gradient_matches_numeric("pbe"); } -TEST(NuclearDerivativeCalculatorTest, - QdkAnalyticGradientMatchesNumericRangeSeparatedHybridDftForWater) { - expect_qdk_analytic_gradient_matches_numeric("wB97x"); +TEST( + NuclearDerivativeCalculatorTest, + QdkAnalyticGradientMatchesNumericRangeSeparatedHybridDftForLithiumHydride) { + expect_qdk_analytic_gradient_matches_numeric( + "wB97x", kDftAnalyticNumericGradientTolerance); } diff --git a/cpp/tests/ut_common.hpp b/cpp/tests/ut_common.hpp index 546a9761c..2a2543315 100644 --- a/cpp/tests/ut_common.hpp +++ b/cpp/tests/ut_common.hpp @@ -234,6 +234,24 @@ inline std::shared_ptr create_li_structure() { return std::make_shared(coords, elements); } +/** + * @brief Creates a LiH structure + */ +inline std::shared_ptr create_lithium_hydride_structure() { + std::vector coords = { + {0.000000000, 0.000000000, 0.000000000}, + {0.000000000, 0.000000000, 1.595000000}}; + + for (auto& coord : coords) { + coord *= qdk::chemistry::constants::angstrom_to_bohr; + } + + std::vector elements = {qdk::chemistry::data::Element::Li, + qdk::chemistry::data::Element::H}; + + return std::make_shared(coords, elements); +} + /** * @brief Creates an O2 structure */ From 479b9f15524abd7bf8eacc3f8ffa88a9a334f637 Mon Sep 17 00:00:00 2001 From: Jan Unsleber Date: Sat, 6 Jun 2026 18:43:39 +0200 Subject: [PATCH 08/13] resolve comment --- .../src/eri/LIBINT2_DIRECT/libint2_direct.cpp | 30 +++++++++++++++++-- 1 file changed, 28 insertions(+), 2 deletions(-) diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/eri/LIBINT2_DIRECT/libint2_direct.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/eri/LIBINT2_DIRECT/libint2_direct.cpp index db04e5768..4ddcbded3 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/eri/LIBINT2_DIRECT/libint2_direct.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/eri/LIBINT2_DIRECT/libint2_direct.cpp @@ -683,6 +683,25 @@ class ERI { P_J = P_total.data(); } + const auto ln_max_engine_precision = std::log(max_engine_precision); + std::vector>> + shell_pair_data(nsh, std::vector>( + nsh, nullptr)); + for (size_t s1 = 0; s1 < nsh; ++s1) { + const auto& shell_pairs = splist_.at(s1); + const auto& shell_pair_data_for_s1 = spdata_.at(s1); + for (size_t pair_index = 0; pair_index < shell_pairs.size(); + ++pair_index) { + const size_t s2 = shell_pairs[pair_index]; + shell_pair_data[s1][s2] = shell_pair_data_for_s1[pair_index]; + if (s1 != s2) { + shell_pair_data[s2][s1] = std::make_shared<::libint2::ShellPair>( + obs_[s2], obs_[s1], ln_max_engine_precision, + ::libint2::ScreeningMethod::Original); + } + } + } + #ifdef _OPENMP const int nthreads = omp_get_max_threads(); #else @@ -729,12 +748,17 @@ class ERI { for (size_t s2 = 0; s2 < nsh; ++s2) { const auto bf2_st = shell2bf_[s2]; const auto n2 = obs_[s2].size(); + const auto* sp12_data = shell_pair_data[s1][s2].get(); for (size_t s3 = 0; s3 < nsh; ++s3) { const auto bf3_st = shell2bf_[s3]; const auto n3 = obs_[s3].size(); for (size_t s4 = 0; s4 < nsh; ++s4, ++s1234) { if (s1234 % nthreads != static_cast(thread_id)) continue; + if (!sp12_data) continue; + const auto* sp34_data = shell_pair_data[s3][s4].get(); + if (!sp34_data) continue; + if (K_schwarz_(s1, s2) * K_schwarz_(s3, s4) < eri_threshold_) continue; @@ -744,12 +768,14 @@ class ERI { if (need_coulomb) { engine.compute2<::libint2::Operator::coulomb, ::libint2::BraKet::xx_xx, 1>( - obs_[s1], obs_[s2], obs_[s3], obs_[s4]); + obs_[s1], obs_[s2], obs_[s3], obs_[s4], sp12_data, + sp34_data); } if (need_erf_exchange) { engine_erf->compute2<::libint2::Operator::erf_coulomb, ::libint2::BraKet::xx_xx, 1>( - obs_[s1], obs_[s2], obs_[s3], obs_[s4]); + obs_[s1], obs_[s2], obs_[s3], obs_[s4], sp12_data, + sp34_data); } const bool has_coulomb = need_coulomb && buf[0] != nullptr; From 2ce3fac45666cc1f9664597caafe11dc4888aee8 Mon Sep 17 00:00:00 2001 From: Jan Unsleber Date: Sat, 6 Jun 2026 23:05:01 +0200 Subject: [PATCH 09/13] fixes --- python/src/pybind11/data/nuclear_gradients.cpp | 8 ++++++-- python/src/pybind11/data/nuclear_hessian.cpp | 8 ++++++-- 2 files changed, 12 insertions(+), 4 deletions(-) diff --git a/python/src/pybind11/data/nuclear_gradients.cpp b/python/src/pybind11/data/nuclear_gradients.cpp index b99c8e5e9..99d09c7c6 100644 --- a/python/src/pybind11/data/nuclear_gradients.cpp +++ b/python/src/pybind11/data/nuclear_gradients.cpp @@ -90,8 +90,12 @@ Create nuclear gradients for a structure. R"(Return the atom-major gradient vector in Hartree/Bohr.)") .def("as_matrix", &NuclearGradients::as_matrix, R"(Return gradients as an ``(num_atoms, 3)`` matrix.)") - .def("get_num_atoms", &NuclearGradients::get_num_atoms, - R"(Return the number of atoms in the associated structure.)") + .def( + "get_num_atoms", + [](const NuclearGradients& gradients) { + return gradients.get_structure()->get_num_atoms(); + }, + R"(Return the number of atoms in the associated structure.)") .def("get_data_type_name", &NuclearGradients::get_data_type_name, R"(Return the serialized data type name.)") .def("get_summary", &NuclearGradients::get_summary, diff --git a/python/src/pybind11/data/nuclear_hessian.cpp b/python/src/pybind11/data/nuclear_hessian.cpp index 9dba729e9..1e0194281 100644 --- a/python/src/pybind11/data/nuclear_hessian.cpp +++ b/python/src/pybind11/data/nuclear_hessian.cpp @@ -85,8 +85,12 @@ Create a nuclear Hessian for a structure. .def("get_matrix", &NuclearHessian::get_matrix, py::return_value_policy::reference_internal, R"(Return the Hessian matrix in Hartree/Bohr^2.)") - .def("get_num_atoms", &NuclearHessian::get_num_atoms, - R"(Return the number of atoms in the associated structure.)") + .def( + "get_num_atoms", + [](const NuclearHessian& hessian) { + return hessian.get_structure()->get_num_atoms(); + }, + R"(Return the number of atoms in the associated structure.)") .def("get_data_type_name", &NuclearHessian::get_data_type_name, R"(Return the serialized data type name.)") .def("get_summary", &NuclearHessian::get_summary, From 9af1682d7fc059eb8b719d9289b9aa0212d0a906 Mon Sep 17 00:00:00 2001 From: Jan Unsleber Date: Sat, 6 Jun 2026 23:58:17 +0200 Subject: [PATCH 10/13] comments --- .../algorithms/nuclear_derivative.cpp | 58 ++++++++++++++++++- cpp/tests/test_nuclear_derivative.cpp | 40 ++++++++----- cpp/tests/ut_common.hpp | 2 +- .../src/pybind11/data/nuclear_gradients.cpp | 6 -- python/src/pybind11/data/nuclear_hessian.cpp | 6 -- python/tests/test_nuclear_derivative.py | 9 ++- 6 files changed, 90 insertions(+), 31 deletions(-) diff --git a/cpp/src/qdk/chemistry/algorithms/nuclear_derivative.cpp b/cpp/src/qdk/chemistry/algorithms/nuclear_derivative.cpp index 2a38b567c..fb74df7dd 100644 --- a/cpp/src/qdk/chemistry/algorithms/nuclear_derivative.cpp +++ b/cpp/src/qdk/chemistry/algorithms/nuclear_derivative.cpp @@ -13,6 +13,10 @@ #include #include +#ifdef QDK_CHEMISTRY_ENABLE_MPI +#include +#endif + #include "microsoft/scf.hpp" namespace qdk::chemistry::algorithms { @@ -64,6 +68,54 @@ std::shared_ptr copy_structure( return std::make_shared(*structure); } +std::optional broadcast_gradient_from_root( + const std::optional& gradient) { +#ifdef QDK_CHEMISTRY_ENABLE_MPI + int initialized = 0; + MPI_Initialized(&initialized); + if (!initialized) { + return gradient; + } + + int finalized = 0; + MPI_Finalized(&finalized); + if (finalized) { + return gradient; + } + + int world_size = 1; + MPI_Comm_size(MPI_COMM_WORLD, &world_size); + if (world_size <= 1) { + return gradient; + } + + int world_rank = 0; + MPI_Comm_rank(MPI_COMM_WORLD, &world_rank); + + int gradient_size = -1; + if (world_rank == 0 && gradient.has_value()) { + gradient_size = static_cast(gradient->size()); + } + MPI_Bcast(&gradient_size, 1, MPI_INT, 0, MPI_COMM_WORLD); + + if (gradient_size < 0) { + return std::nullopt; + } + + Eigen::VectorXd broadcast_gradient; + if (world_rank == 0) { + broadcast_gradient = *gradient; + } else { + broadcast_gradient.resize(gradient_size); + } + MPI_Bcast(broadcast_gradient.data(), gradient_size, MPI_DOUBLE, 0, + MPI_COMM_WORLD); + return broadcast_gradient; +#else + return gradient; +#endif +} + std::shared_ptr displace_structure( const std::shared_ptr& structure, Eigen::Index coordinate, double displacement) { @@ -505,12 +557,14 @@ NuclearDerivativeResult QdkNuclearDerivativeCalculator::_run_impl( auto scf_result = solver.run_with_analytic_gradient( structure, charge, spin_multiplicity, seed_to_scf_input(seed, true)); - if (!scf_result.nuclear_gradient.has_value()) { + auto nuclear_gradient = + broadcast_gradient_from_root(scf_result.nuclear_gradient); + if (!nuclear_gradient.has_value()) { throw std::runtime_error( "Internal SCF did not return the requested analytic nuclear gradient"); } auto gradients = std::make_shared( - copy_structure(structure), *scf_result.nuclear_gradient); + copy_structure(structure), *nuclear_gradient); return {scf_result.energy, gradients, std::nullopt, scf_result.wavefunction}; } diff --git a/cpp/tests/test_nuclear_derivative.cpp b/cpp/tests/test_nuclear_derivative.cpp index 082d1dadc..f9b39817d 100644 --- a/cpp/tests/test_nuclear_derivative.cpp +++ b/cpp/tests/test_nuclear_derivative.cpp @@ -14,6 +14,7 @@ #include #include #include +#include #include #include "ut_common.hpp" @@ -90,6 +91,23 @@ class RecordingLocalizer : public Localizer { } }; +template +class ScopedFactoryRegistration { + public: + explicit ScopedFactoryRegistration(std::string key) : key_(std::move(key)) { + Factory::unregister_instance(key_); + } + + ~ScopedFactoryRegistration() { Factory::unregister_instance(key_); } + + ScopedFactoryRegistration(const ScopedFactoryRegistration&) = delete; + ScopedFactoryRegistration& operator=(const ScopedFactoryRegistration&) = + delete; + + private: + std::string key_; +}; + void expect_same_structure(const std::shared_ptr& left, const std::shared_ptr& right) { ASSERT_NE(left, nullptr); @@ -110,7 +128,7 @@ NuclearDerivativeResult run_qdk_derivative_for_functional( calculator->settings().set("finite_difference_step", 1.0e-2); } - return calculator->run(testing::create_lithium_hydride_structure(), 0, 1, + return calculator->run(testing::create_lih_structure(), 0, 1, std::string("sto-3g")); } @@ -263,8 +281,9 @@ TEST(NuclearDerivativeCalculatorTest, MultiReferenceFiniteDifferenceReusesSeedActiveSpace) { recorded_active_spaces.clear(); recorded_inactive_spaces.clear(); - MultiConfigurationCalculatorFactory::unregister_instance( - "_test_nuclear_derivative_recording_mc"); + [[maybe_unused]] ScopedFactoryRegistration< + MultiConfigurationCalculatorFactory> + mc_guard("_test_nuclear_derivative_recording_mc"); MultiConfigurationCalculatorFactory::register_instance( []() -> MultiConfigurationCalculatorFactory::return_type { return std::make_unique(); @@ -300,9 +319,6 @@ TEST(NuclearDerivativeCalculatorTest, for (const auto& inactive_space : recorded_inactive_spaces) { EXPECT_TRUE(inactive_space.empty()); } - - MultiConfigurationCalculatorFactory::unregister_instance( - "_test_nuclear_derivative_recording_mc"); } TEST(NuclearDerivativeCalculatorTest, @@ -311,13 +327,14 @@ TEST(NuclearDerivativeCalculatorTest, recorded_inactive_spaces.clear(); recorded_localized_alpha_spaces.clear(); recorded_localized_beta_spaces.clear(); - MultiConfigurationCalculatorFactory::unregister_instance( - "_test_nuclear_derivative_recording_mc"); + [[maybe_unused]] ScopedFactoryRegistration< + MultiConfigurationCalculatorFactory> + mc_guard("_test_nuclear_derivative_recording_mc"); MultiConfigurationCalculatorFactory::register_instance( []() -> MultiConfigurationCalculatorFactory::return_type { return std::make_unique(); }); - LocalizerFactory::unregister_instance( + [[maybe_unused]] ScopedFactoryRegistration localizer_guard( "_test_nuclear_derivative_recording_localizer"); LocalizerFactory::register_instance([]() -> LocalizerFactory::return_type { return std::make_unique(); @@ -360,11 +377,6 @@ TEST(NuclearDerivativeCalculatorTest, for (const auto& localized_space : recorded_localized_beta_spaces) { EXPECT_EQ(localized_space, std::vector{0}); } - - LocalizerFactory::unregister_instance( - "_test_nuclear_derivative_recording_localizer"); - MultiConfigurationCalculatorFactory::unregister_instance( - "_test_nuclear_derivative_recording_mc"); } TEST(NuclearDerivativeCalculatorTest, diff --git a/cpp/tests/ut_common.hpp b/cpp/tests/ut_common.hpp index 2a2543315..d3a9a552e 100644 --- a/cpp/tests/ut_common.hpp +++ b/cpp/tests/ut_common.hpp @@ -237,7 +237,7 @@ inline std::shared_ptr create_li_structure() { /** * @brief Creates a LiH structure */ -inline std::shared_ptr create_lithium_hydride_structure() { +inline std::shared_ptr create_lih_structure() { std::vector coords = { {0.000000000, 0.000000000, 0.000000000}, {0.000000000, 0.000000000, 1.595000000}}; diff --git a/python/src/pybind11/data/nuclear_gradients.cpp b/python/src/pybind11/data/nuclear_gradients.cpp index 99d09c7c6..001315ecc 100644 --- a/python/src/pybind11/data/nuclear_gradients.cpp +++ b/python/src/pybind11/data/nuclear_gradients.cpp @@ -90,12 +90,6 @@ Create nuclear gradients for a structure. R"(Return the atom-major gradient vector in Hartree/Bohr.)") .def("as_matrix", &NuclearGradients::as_matrix, R"(Return gradients as an ``(num_atoms, 3)`` matrix.)") - .def( - "get_num_atoms", - [](const NuclearGradients& gradients) { - return gradients.get_structure()->get_num_atoms(); - }, - R"(Return the number of atoms in the associated structure.)") .def("get_data_type_name", &NuclearGradients::get_data_type_name, R"(Return the serialized data type name.)") .def("get_summary", &NuclearGradients::get_summary, diff --git a/python/src/pybind11/data/nuclear_hessian.cpp b/python/src/pybind11/data/nuclear_hessian.cpp index 1e0194281..fc1cb132a 100644 --- a/python/src/pybind11/data/nuclear_hessian.cpp +++ b/python/src/pybind11/data/nuclear_hessian.cpp @@ -85,12 +85,6 @@ Create a nuclear Hessian for a structure. .def("get_matrix", &NuclearHessian::get_matrix, py::return_value_policy::reference_internal, R"(Return the Hessian matrix in Hartree/Bohr^2.)") - .def( - "get_num_atoms", - [](const NuclearHessian& hessian) { - return hessian.get_structure()->get_num_atoms(); - }, - R"(Return the number of atoms in the associated structure.)") .def("get_data_type_name", &NuclearHessian::get_data_type_name, R"(Return the serialized data type name.)") .def("get_summary", &NuclearHessian::get_summary, diff --git a/python/tests/test_nuclear_derivative.py b/python/tests/test_nuclear_derivative.py index 104f51ecf..139afa5cf 100644 --- a/python/tests/test_nuclear_derivative.py +++ b/python/tests/test_nuclear_derivative.py @@ -11,16 +11,18 @@ def _h2_structure(): + """Create an H2 structure for derivative data tests.""" return data.Structure([[0.0, 0.0, 0.0], [0.0, 0.0, 1.4]], [1, 1]) def test_nuclear_gradients_roundtrip_json(): + """Round-trip nuclear gradients through JSON.""" structure = _h2_structure() values = np.arange(6, dtype=float) gradients = data.NuclearGradients(structure, values) - assert gradients.get_num_atoms() == 2 + assert gradients.get_structure().get_num_atoms() == 2 np.testing.assert_allclose(gradients.get_values(), values) assert gradients.as_matrix().shape == (2, 3) @@ -33,12 +35,13 @@ def test_nuclear_gradients_roundtrip_json(): def test_nuclear_hessian_roundtrip_json(): + """Round-trip a nuclear Hessian through JSON.""" structure = _h2_structure() matrix = np.eye(6) hessian = data.NuclearHessian(structure, matrix) - assert hessian.get_num_atoms() == 2 + assert hessian.get_structure().get_num_atoms() == 2 np.testing.assert_allclose(hessian.get_matrix(), matrix) loaded = data.NuclearHessian.from_json(hessian.to_json()) @@ -50,6 +53,7 @@ def test_nuclear_hessian_roundtrip_json(): def test_nuclear_derivative_factory_registered(): + """Create the default nuclear derivative calculator from the factory.""" calculator = algorithms.create("nuclear_derivative_calculator") assert isinstance(calculator, algorithms.NuclearDerivativeCalculator) @@ -57,6 +61,7 @@ def test_nuclear_derivative_factory_registered(): def test_qdk_nuclear_derivative_factory_registered(): + """Create the QDK nuclear derivative calculator from the factory.""" calculator = algorithms.create("nuclear_derivative_calculator", "qdk") assert isinstance(calculator, algorithms.QdkNuclearDerivativeCalculator) From bff84376ac79ffa84a09257af282136a33989bde Mon Sep 17 00:00:00 2001 From: Jan Unsleber Date: Mon, 8 Jun 2026 14:39:48 +0200 Subject: [PATCH 11/13] initial draft --- cpp/include/qdk/chemistry.hpp | 1 + .../algorithms/geometry_optimization.hpp | 128 ++++++++++++++ .../qdk/chemistry/algorithms/CMakeLists.txt | 1 + .../algorithms/algorithm_defaults.cpp | 2 + .../algorithms/geometry_optimization.cpp | 11 ++ python/CMakeLists.txt | 1 + python/pyproject.toml | 1 + .../algorithms/geometry_optimization.cpp | 99 +++++++++++ python/src/pybind11/module.cpp | 2 + python/src/qdk_chemistry/__init__.py | 4 + .../src/qdk_chemistry/algorithms/__init__.py | 2 + .../algorithms/geometry_optimization.py | 10 ++ .../src/qdk_chemistry/algorithms/registry.py | 2 + .../plugins/geometric/__init__.py | 36 ++++ .../plugins/geometric/geometry_optimizer.py | 165 ++++++++++++++++++ python/tests/test_algorithms_registry.py | 7 + python/tests/test_geometric_plugin.py | 53 ++++++ python/tests/test_nuclear_derivative.py | 8 + 18 files changed, 533 insertions(+) create mode 100644 cpp/include/qdk/chemistry/algorithms/geometry_optimization.hpp create mode 100644 cpp/src/qdk/chemistry/algorithms/geometry_optimization.cpp create mode 100644 python/src/pybind11/algorithms/geometry_optimization.cpp create mode 100644 python/src/qdk_chemistry/algorithms/geometry_optimization.py create mode 100644 python/src/qdk_chemistry/plugins/geometric/__init__.py create mode 100644 python/src/qdk_chemistry/plugins/geometric/geometry_optimizer.py create mode 100644 python/tests/test_geometric_plugin.py diff --git a/cpp/include/qdk/chemistry.hpp b/cpp/include/qdk/chemistry.hpp index 9392a6a98..b3fe3ec80 100644 --- a/cpp/include/qdk/chemistry.hpp +++ b/cpp/include/qdk/chemistry.hpp @@ -4,6 +4,7 @@ #include #include +#include #include #include #include diff --git a/cpp/include/qdk/chemistry/algorithms/geometry_optimization.hpp b/cpp/include/qdk/chemistry/algorithms/geometry_optimization.hpp new file mode 100644 index 000000000..31d7186b7 --- /dev/null +++ b/cpp/include/qdk/chemistry/algorithms/geometry_optimization.hpp @@ -0,0 +1,128 @@ +// Copyright (c) Microsoft Corporation. All rights reserved. +// Licensed under the MIT License. See LICENSE.txt in the project root for +// license information. + +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace qdk::chemistry::algorithms { + +/** + * @brief Energy, optimized structure, optional wavefunction, and optional + * Hessian from a geometry optimization run. + */ +using GeometryOptimizationResult = + std::tuple, + std::optional>, + std::optional>>; + +/** + * @class GeometryOptimizerSettings + * @brief Shared settings for geometry optimization algorithms. + */ +class GeometryOptimizerSettings : public data::Settings { + public: + /** + * @brief Construct geometry optimization settings with common defaults. + */ + GeometryOptimizerSettings() { + set_default("derivative_calculator", + data::AlgorithmRef("nuclear_derivative_calculator", + "finite_difference"), + "Nuclear derivative calculator used to evaluate energies and " + "gradients during optimization."); + set_default("transition_state", false, + "Whether to run transition-state optimization instead of " + "minimum-energy geometry optimization."); + set_default("max_iterations", static_cast(300), + "Maximum number of geometry optimization steps.", + data::BoundConstraint{1, 1000000}); + set_default("convergence_energy", 1.0e-6, + "Energy convergence threshold used by optimizers that expose " + "an energy convergence setting.", + data::BoundConstraint{0.0, 1.0}); + set_default("convergence_gradient", 3.0e-4, + "Gradient convergence threshold used by optimizers that expose " + "a gradient convergence setting.", + data::BoundConstraint{0.0, 1.0}); + set_default("convergence_displacement", 1.2e-3, + "Displacement convergence threshold used by optimizers that " + "expose a displacement convergence setting.", + data::BoundConstraint{0.0, 1.0}); + set_default("compute_hessian", false, + "Whether to compute a nuclear Hessian at the optimized " + "geometry before returning."); + } +}; + +/** + * @class GeometryOptimizer + * @brief Base class for geometry optimization algorithms. + * + * Implementations optimize molecular coordinates using the same input + * arguments as nuclear derivative calculators. Results include the optimized + * energy and structure plus optional wavefunction and Hessian values. + */ +class GeometryOptimizer + : public Algorithm, int, int, + NuclearDerivativeSeedType> { + public: + /** + * @brief Construct a geometry optimizer with shared settings. + */ + GeometryOptimizer() { + _settings = std::make_unique(); + } + virtual ~GeometryOptimizer() = default; + + using Algorithm::run; + + virtual std::string name() const = 0; + + /** + * @brief Return the factory type name for geometry optimizers. + */ + std::string type_name() const final { return "geometry_optimizer"; } + + protected: + /** + * @brief Implementation hook for derived geometry optimizers. + */ + virtual GeometryOptimizationResult _run_impl( + std::shared_ptr structure, int charge, + int spin_multiplicity, NuclearDerivativeSeedType seed) const = 0; +}; + +/** + * @brief Factory for geometry optimizer implementations. + */ +struct GeometryOptimizerFactory + : public AlgorithmFactory { + /** + * @brief Return the algorithm type name managed by this factory. + */ + static std::string algorithm_type_name() { return "geometry_optimizer"; } + + /** + * @brief Register built-in geometry optimizer implementations. + */ + static void register_default_instances(); + + /** + * @brief Return the default geometry optimizer implementation name. + */ + static std::string default_algorithm_name() { return "geometric"; } +}; + +} // namespace qdk::chemistry::algorithms diff --git a/cpp/src/qdk/chemistry/algorithms/CMakeLists.txt b/cpp/src/qdk/chemistry/algorithms/CMakeLists.txt index c06f2abae..dc026562f 100644 --- a/cpp/src/qdk/chemistry/algorithms/CMakeLists.txt +++ b/cpp/src/qdk/chemistry/algorithms/CMakeLists.txt @@ -7,6 +7,7 @@ target_sources(chemistry PRIVATE nuclear_derivative.cpp pmc.cpp dynamical_correlation_calculator.cpp + geometry_optimization.cpp scf.cpp stability.cpp ) diff --git a/cpp/src/qdk/chemistry/algorithms/algorithm_defaults.cpp b/cpp/src/qdk/chemistry/algorithms/algorithm_defaults.cpp index cbc009ac6..2f8e80d4d 100644 --- a/cpp/src/qdk/chemistry/algorithms/algorithm_defaults.cpp +++ b/cpp/src/qdk/chemistry/algorithms/algorithm_defaults.cpp @@ -5,6 +5,7 @@ #include #include #include +#include #include #include #include @@ -41,6 +42,7 @@ std::shared_ptr resolve_algorithm_defaults( REGISTER_FACTORY_SETTINGS_INIT(MultiConfigurationCalculatorFactory) REGISTER_FACTORY_SETTINGS_INIT(ProjectedMultiConfigurationCalculatorFactory) REGISTER_FACTORY_SETTINGS_INIT(DynamicalCorrelationCalculatorFactory) + REGISTER_FACTORY_SETTINGS_INIT(GeometryOptimizerFactory) REGISTER_FACTORY_SETTINGS_INIT(MultiConfigurationScfFactory) REGISTER_FACTORY_SETTINGS_INIT(NuclearDerivativeCalculatorFactory) REGISTER_FACTORY_SETTINGS_INIT(LocalizerFactory) diff --git a/cpp/src/qdk/chemistry/algorithms/geometry_optimization.cpp b/cpp/src/qdk/chemistry/algorithms/geometry_optimization.cpp new file mode 100644 index 000000000..d77d407ca --- /dev/null +++ b/cpp/src/qdk/chemistry/algorithms/geometry_optimization.cpp @@ -0,0 +1,11 @@ +// Copyright (c) Microsoft Corporation. All rights reserved. +// Licensed under the MIT License. See LICENSE.txt in the project root for +// license information. + +#include + +namespace qdk::chemistry::algorithms { + +void GeometryOptimizerFactory::register_default_instances() {} + +} // namespace qdk::chemistry::algorithms diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index d28b861e4..f0df2e37f 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -142,6 +142,7 @@ pybind11_add_module(_core src/pybind11/algorithms/hamiltonian.cpp src/pybind11/algorithms/scf.cpp src/pybind11/algorithms/nuclear_derivative.cpp + src/pybind11/algorithms/geometry_optimization.cpp src/pybind11/algorithms/active_space.cpp src/pybind11/algorithms/dynamical_correlation_calculator.cpp src/pybind11/algorithms/davidson_solver.cpp diff --git a/python/pyproject.toml b/python/pyproject.toml index 6e9a66ba1..32d9bca44 100644 --- a/python/pyproject.toml +++ b/python/pyproject.toml @@ -87,6 +87,7 @@ openfermion-extras = [ "openfermion>=1.0.0; python_version < '3.14'" ] plugins = [ + "geometric>=1.0", "pyscf>=2.9.0,<2.12.1" ] qiskit-extras = [ diff --git a/python/src/pybind11/algorithms/geometry_optimization.cpp b/python/src/pybind11/algorithms/geometry_optimization.cpp new file mode 100644 index 000000000..dd6cad6d9 --- /dev/null +++ b/python/src/pybind11/algorithms/geometry_optimization.cpp @@ -0,0 +1,99 @@ +// Copyright (c) Microsoft Corporation. All rights reserved. +// Licensed under the MIT License. See LICENSE.txt in the project root for +// license information. + +#include +#include + +#include + +#include "factory_bindings.hpp" + +namespace py = pybind11; +using namespace qdk::chemistry::algorithms; +using namespace qdk::chemistry::data; + +class GeometryOptimizerBase : public GeometryOptimizer, + public pybind11::trampoline_self_life_support { + public: + std::string name() const override { + PYBIND11_OVERRIDE_PURE(std::string, GeometryOptimizer, name); + } + + std::vector aliases() const override { + PYBIND11_OVERRIDE(std::vector, GeometryOptimizer, aliases); + } + + void replace_settings( + std::unique_ptr new_settings) { + this->_settings = std::move(new_settings); + } + + protected: + GeometryOptimizationResult _run_impl( + std::shared_ptr structure, int charge, int spin_multiplicity, + NuclearDerivativeSeedType seed) const override { + PYBIND11_OVERRIDE_PURE(GeometryOptimizationResult, GeometryOptimizer, + _run_impl, structure, charge, spin_multiplicity, + seed); + } +}; + +void bind_geometry_optimization(py::module& m) { + py::class_ + optimizer(m, "GeometryOptimizer", + R"( + Base class for geometry optimization algorithms. + + Optimizers take the same arguments as nuclear derivative calculators and + return the optimized energy, optimized structure, optional wavefunction, + and optional Hessian. + )"); + + optimizer.def(py::init<>(), R"(Create a geometry optimizer.)"); + optimizer.def( + "run", + [](const GeometryOptimizer& self, std::shared_ptr structure, + int charge, int spin_multiplicity, NuclearDerivativeSeedType seed) { + return self.run(structure, charge, spin_multiplicity, seed); + }, + py::arg("structure"), py::arg("charge"), py::arg("spin_multiplicity"), + py::arg("seed_or_basis"), + R"( +Optimize a molecular structure. + +Args: + structure: Initial molecular structure to optimize. + charge: Total molecular charge. + spin_multiplicity: Spin multiplicity of the molecular system. + seed_or_basis: Basis name, basis set, orbitals, or wavefunction seed. + +Returns: + tuple: ``(energy, structure, wavefunction, hessian)``. +)"); + optimizer.def("settings", &GeometryOptimizer::settings, + py::return_value_policy::reference_internal, + R"(Return the optimizer settings.)"); + optimizer.def_property( + "_settings", + [](GeometryOptimizerBase& self) -> Settings& { return self.settings(); }, + [](GeometryOptimizerBase& self, + std::unique_ptr new_settings) { + self.replace_settings(std::move(new_settings)); + }, + py::return_value_policy::reference_internal, + R"(Internal settings replacement hook for Python subclasses.)"); + optimizer.def("name", &GeometryOptimizer::name, + R"(Return the implementation name.)"); + optimizer.def("type_name", &GeometryOptimizer::type_name, + R"(Return the algorithm type name.)"); + optimizer.def("__repr__", [](const GeometryOptimizer& self) { + return ""; + }); + + qdk::chemistry::python::bind_algorithm_factory< + GeometryOptimizerFactory, GeometryOptimizer, GeometryOptimizerBase>( + m, "GeometryOptimizerFactory"); + qdk::chemistry::python::bind_create_nested(optimizer); +} diff --git a/python/src/pybind11/module.cpp b/python/src/pybind11/module.cpp index 6228e8e9b..f99d8f57b 100644 --- a/python/src/pybind11/module.cpp +++ b/python/src/pybind11/module.cpp @@ -29,6 +29,7 @@ void bind_mcscf(py::module& m); void bind_hamiltonian_constructor(py::module& m); void bind_scf(py::module& m); void bind_nuclear_derivative(py::module& m); +void bind_geometry_optimization(py::module& m); void bind_active_space(py::module& m); void bind_constants(py::module& m); void bind_pmc(py::module& m); @@ -85,6 +86,7 @@ PYBIND11_MODULE(_core, m) { bind_hamiltonian_constructor(algorithms); bind_scf(algorithms); bind_nuclear_derivative(algorithms); + bind_geometry_optimization(algorithms); bind_active_space(algorithms); bind_dynamical_correlation_calculator(algorithms); bind_pmc(algorithms); diff --git a/python/src/qdk_chemistry/__init__.py b/python/src/qdk_chemistry/__init__.py index 3d448916f..e4e19d27b 100644 --- a/python/src/qdk_chemistry/__init__.py +++ b/python/src/qdk_chemistry/__init__.py @@ -134,6 +134,10 @@ def _import_plugins() -> None: import qdk_chemistry.plugins.networkx as networkx_plugin # noqa: PLC0415 networkx_plugin.load() + with contextlib.suppress(ImportError): + import qdk_chemistry.plugins.geometric as geometric_plugin # noqa: PLC0415 + + geometric_plugin.load() def _is_placeholder_stub(stub_file: Path) -> bool: diff --git a/python/src/qdk_chemistry/algorithms/__init__.py b/python/src/qdk_chemistry/algorithms/__init__.py index dc90d3beb..87d96d8da 100644 --- a/python/src/qdk_chemistry/algorithms/__init__.py +++ b/python/src/qdk_chemistry/algorithms/__init__.py @@ -26,6 +26,7 @@ from qdk_chemistry.algorithms.dynamical_correlation_calculator import DynamicalCorrelationCalculator from qdk_chemistry.algorithms.energy_estimator.energy_estimator import EnergyEstimator from qdk_chemistry.algorithms.energy_estimator.qdk import QdkEnergyEstimator +from qdk_chemistry.algorithms.geometry_optimization import GeometryOptimizer from qdk_chemistry.algorithms.hamiltonian_constructor import ( HamiltonianConstructor, QdkHamiltonianConstructor, @@ -70,6 +71,7 @@ "DynamicalCorrelationCalculator", "EnergyEstimator", "FiniteDifferenceNuclearDerivativeCalculator", + "GeometryOptimizer", "HamiltonianConstructor", "HamiltonianUnitaryBuilder", "MultiConfigurationCalculator", diff --git a/python/src/qdk_chemistry/algorithms/geometry_optimization.py b/python/src/qdk_chemistry/algorithms/geometry_optimization.py new file mode 100644 index 000000000..c2a585010 --- /dev/null +++ b/python/src/qdk_chemistry/algorithms/geometry_optimization.py @@ -0,0 +1,10 @@ +"""Public entry point for geometry optimizer algorithms.""" + +# -------------------------------------------------------------------------------------------- +# Copyright (c) Microsoft Corporation. All rights reserved. +# Licensed under the MIT License. See LICENSE.txt in the project root for license information. +# -------------------------------------------------------------------------------------------- + +from qdk_chemistry._core._algorithms import GeometryOptimizer + +__all__ = ["GeometryOptimizer"] diff --git a/python/src/qdk_chemistry/algorithms/registry.py b/python/src/qdk_chemistry/algorithms/registry.py index a03d38790..675896d83 100644 --- a/python/src/qdk_chemistry/algorithms/registry.py +++ b/python/src/qdk_chemistry/algorithms/registry.py @@ -473,6 +473,7 @@ def _register_cpp_factories(): from qdk_chemistry._core._algorithms import ( # noqa: PLC0415 ActiveSpaceSelectorFactory, DynamicalCorrelationCalculatorFactory, + GeometryOptimizerFactory, HamiltonianConstructorFactory, LocalizerFactory, MultiConfigurationCalculatorFactory, @@ -491,6 +492,7 @@ def _register_cpp_factories(): register_factory(NuclearDerivativeCalculatorFactory) register_factory(ProjectedMultiConfigurationCalculatorFactory) register_factory(DynamicalCorrelationCalculatorFactory) + register_factory(GeometryOptimizerFactory) register_factory(ScfSolverFactory) register_factory(StabilityCheckerFactory) diff --git a/python/src/qdk_chemistry/plugins/geometric/__init__.py b/python/src/qdk_chemistry/plugins/geometric/__init__.py new file mode 100644 index 000000000..876d62ce1 --- /dev/null +++ b/python/src/qdk_chemistry/plugins/geometric/__init__.py @@ -0,0 +1,36 @@ +"""geomeTRIC plugin for QDK/Chemistry geometry optimization.""" + +# -------------------------------------------------------------------------------------------- +# Copyright (c) Microsoft Corporation. All rights reserved. +# Licensed under the MIT License. See LICENSE.txt in the project root for license information. +# -------------------------------------------------------------------------------------------- + +import importlib.util + +from qdk_chemistry.utils import Logger + +_loaded = False +QDK_CHEMISTRY_HAS_GEOMETRIC = False + + +def load(): + """Load the geomeTRIC plugin into QDK/Chemistry.""" + Logger.trace_entering() + global _loaded, QDK_CHEMISTRY_HAS_GEOMETRIC # noqa: PLW0603 + if _loaded: + return + _loaded = True + + if importlib.util.find_spec("geometric") is not None: + QDK_CHEMISTRY_HAS_GEOMETRIC = True + _register_algorithms() + + +def _register_algorithms(): + """Register geomeTRIC-backed optimizer algorithms.""" + Logger.trace_entering() + from qdk_chemistry.algorithms import register # noqa: PLC0415 + from qdk_chemistry.plugins.geometric.geometry_optimizer import GeometricOptimizer # noqa: PLC0415 + + register(lambda: GeometricOptimizer()) + Logger.debug(f"geomeTRIC plugin loaded: [{GeometricOptimizer().type_name()}: {GeometricOptimizer().name()}].") diff --git a/python/src/qdk_chemistry/plugins/geometric/geometry_optimizer.py b/python/src/qdk_chemistry/plugins/geometric/geometry_optimizer.py new file mode 100644 index 000000000..4bd1c7a4f --- /dev/null +++ b/python/src/qdk_chemistry/plugins/geometric/geometry_optimizer.py @@ -0,0 +1,165 @@ +"""geomeTRIC-backed geometry optimizer.""" + +# -------------------------------------------------------------------------------------------- +# Copyright (c) Microsoft Corporation. All rights reserved. +# Licensed under the MIT License. See LICENSE.txt in the project root for license information. +# -------------------------------------------------------------------------------------------- + +from __future__ import annotations + +from tempfile import TemporaryDirectory +from typing import Any + +import numpy as np +from geometric.engine import Engine + +from qdk_chemistry.algorithms import GeometryOptimizer +from qdk_chemistry.data import AlgorithmRef, NuclearHessian, Settings, Structure +from qdk_chemistry.utils import Logger + +__all__ = ["GeometricOptimizer", "GeometricOptimizerSettings"] + + +class GeometricOptimizerSettings(Settings): + """Settings for the geomeTRIC geometry optimizer.""" + + def __init__(self): + """Initialize geomeTRIC optimizer defaults.""" + super().__init__() + self._set_default( + "derivative_calculator", + "algorithm_ref", + AlgorithmRef("nuclear_derivative_calculator", "finite_difference"), + "Nuclear derivative calculator used to evaluate energies and gradients.", + ) + self._set_default( + "transition_state", "bool", False, "Run transition-state optimization instead of minimization." + ) + self._set_default("max_iterations", "int", 300, "Maximum number of geometry optimization steps.") + self._set_default("convergence_energy", "double", 1.0e-6, "Energy convergence threshold.") + self._set_default("convergence_gradient", "double", 3.0e-4, "Gradient convergence threshold.") + self._set_default("convergence_displacement", "double", 1.2e-3, "Displacement convergence threshold.") + self._set_default("compute_hessian", "bool", False, "Compute a Hessian at the optimized geometry.") + self._set_default("print_level", "int", 0, "geomeTRIC output verbosity level.") + + +class _QdkDerivativeEngine(Engine): + """geomeTRIC engine that evaluates QDK/Chemistry nuclear derivatives.""" + + def __init__( + self, + structure: Structure, + charge: int, + spin_multiplicity: int, + seed_or_basis: Any, + derivative_calculator: Any, + ): + super().__init__() + self._structure = structure + self._charge = charge + self._spin_multiplicity = spin_multiplicity + self._seed_or_basis = seed_or_basis + self._derivative_calculator = derivative_calculator + self._last_energy = None + self._last_structure = structure + self._last_wavefunction = None + + def structure_from_coordinates(self, coordinates: np.ndarray) -> Structure: + """Create a QDK/Chemistry structure with updated coordinates.""" + matrix = np.asarray(coordinates, dtype=float).reshape((-1, 3)) + return Structure( + matrix, self._structure.get_elements(), self._structure.get_masses(), self._structure.get_nuclear_charges() + ) + + def last_coordinates(self) -> np.ndarray: + """Return the most recently evaluated coordinates.""" + return np.asarray(self._last_structure.get_coordinates(), dtype=float) + + def calc_new(self, coordinates: np.ndarray, dirname: str) -> dict[str, np.ndarray | float]: # noqa: ARG002 + """Evaluate energy and gradients for geomeTRIC.""" + structure = self.structure_from_coordinates(coordinates) + energy, gradients, _hessian, wavefunction = self._derivative_calculator.run( + structure, self._charge, self._spin_multiplicity, self._seed_or_basis + ) + self._last_energy = energy + self._last_structure = structure + self._last_wavefunction = wavefunction + return {"energy": energy, "gradient": np.asarray(gradients.get_values(), dtype=float)} + + +class GeometricOptimizer(GeometryOptimizer): + """Geometry optimizer implemented with the geomeTRIC Python library.""" + + def __init__(self): + """Initialize the geomeTRIC optimizer.""" + Logger.trace_entering() + super().__init__() + self._settings = GeometricOptimizerSettings() + + def name(self) -> str: + """Return the implementation name.""" + return "geometric" + + def aliases(self) -> list[str]: + """Return accepted factory aliases.""" + return ["geometric", "geomeTRIC"] + + def _run_impl( + self, structure: Structure, charge: int, spin_multiplicity: int, seed_or_basis: Any + ) -> tuple[float, Structure, Any | None, NuclearHessian | None]: + """Optimize a molecular structure using geomeTRIC.""" + Logger.trace_entering() + from geometric.molecule import Molecule # noqa: PLC0415 + from geometric.optimize import run_optimizer # noqa: PLC0415 + + derivative_calculator = self._create_nested("derivative_calculator") + derivative_calculator.settings().set("compute_hessian", False) + engine = _QdkDerivativeEngine(structure, charge, spin_multiplicity, seed_or_basis, derivative_calculator) + + molecule = Molecule() + molecule.elem = structure.get_atomic_symbols() + molecule.xyzs = [np.asarray(structure.get_coordinates(), dtype=float)] + engine.M = molecule + + params = { + "customengine": engine, + "input": None, + "maxiter": self._settings["max_iterations"], + "convergence_energy": self._settings["convergence_energy"], + "convergence_grms": self._settings["convergence_gradient"], + "convergence_drms": self._settings["convergence_displacement"], + "transition": self._settings["transition_state"], + "verbose": self._settings["print_level"], + } + + with TemporaryDirectory(prefix="qdk-chemistry-geometric-") as tmpdir: + result = run_optimizer(**params, dirname=tmpdir) + + optimized_coordinates = _extract_coordinates(result, engine) + optimized_structure = engine.structure_from_coordinates(optimized_coordinates) + + final_calculator = self._create_nested("derivative_calculator") + final_calculator.settings().set("compute_hessian", self._settings["compute_hessian"]) + final_energy, _gradients, hessian, wavefunction = final_calculator.run( + optimized_structure, charge, spin_multiplicity, seed_or_basis + ) + if not self._settings["compute_hessian"]: + hessian = None + + return final_energy, optimized_structure, wavefunction, hessian + + +def _extract_coordinates(result: Any, engine: _QdkDerivativeEngine) -> np.ndarray: + """Extract final coordinates from geomeTRIC's return value.""" + if isinstance(result, np.ndarray): + return result + if hasattr(result, "xyzs") and result.xyzs: + return np.asarray(result.xyzs[-1], dtype=float) + if isinstance(result, dict): + for key in ("coords", "coordinates", "xyz", "xyzs"): + if key in result: + value = result[key] + if key == "xyzs" and value: + value = value[-1] + return np.asarray(value, dtype=float) + return engine.last_coordinates() diff --git a/python/tests/test_algorithms_registry.py b/python/tests/test_algorithms_registry.py index 799383751..8b9574f92 100644 --- a/python/tests/test_algorithms_registry.py +++ b/python/tests/test_algorithms_registry.py @@ -37,6 +37,7 @@ def test_show_default_returns_dict_when_no_type_specified(self): expected_types = [ "active_space_selector", "dynamical_correlation_calculator", + "geometry_optimizer", "hamiltonian_constructor", "orbital_localizer", "multi_configuration_calculator", @@ -125,6 +126,8 @@ def test_show_default_consistent_with_available(self): for algorithm_type, default_name in defaults.items(): # Each default should be in the available list for that type assert algorithm_type in available, f"Algorithm type '{algorithm_type}' not in available algorithms" + if algorithm_type == "geometry_optimizer" and not available[algorithm_type]: + continue # geomeTRIC-backed optimizer is optional if default_name == "pyscf" and not PYSCF_AVAILABLE: continue # Skip check if pyscf is not available assert default_name in available[algorithm_type], ( @@ -137,6 +140,8 @@ def test_default_algorithm_can_be_created(self): defaults = registry.show_default() for algorithm_type, default_name in defaults.items(): + if default_name not in registry.available(algorithm_type): + continue # Optional plugin-backed default is not installed if default_name == "pyscf" and not PYSCF_AVAILABLE: continue # Skip check if pyscf is not available # Should be able to create the default algorithm @@ -253,6 +258,8 @@ def test_available_all_types_have_algorithms(self): all_algorithms = registry.available() for algorithm_type, algorithms in all_algorithms.items(): assert isinstance(algorithms, list), f"Expected list for {algorithm_type}" + if algorithm_type == "geometry_optimizer" and not algorithms: + continue # geomeTRIC-backed optimizer is optional assert len(algorithms) > 0, f"No algorithms available for {algorithm_type}" def test_available_algorithms_can_be_created(self): diff --git a/python/tests/test_geometric_plugin.py b/python/tests/test_geometric_plugin.py new file mode 100644 index 000000000..0b16a9e04 --- /dev/null +++ b/python/tests/test_geometric_plugin.py @@ -0,0 +1,53 @@ +"""Tests for the optional geomeTRIC geometry optimizer plugin.""" + +# -------------------------------------------------------------------------------------------- +# Copyright (c) Microsoft Corporation. All rights reserved. +# Licensed under the MIT License. See LICENSE.txt in the project root for license information. +# -------------------------------------------------------------------------------------------- + +import numpy as np +import pytest + +from qdk_chemistry import algorithms, data +from qdk_chemistry.plugins.geometric import QDK_CHEMISTRY_HAS_GEOMETRIC + +pytestmark = pytest.mark.skipif(not QDK_CHEMISTRY_HAS_GEOMETRIC, reason="geomeTRIC not available") + + +def _h2_structure(): + return data.Structure([[0.0, 0.0, 0.0], [0.0, 0.0, 1.4]], [1, 1]) + + +def test_geometric_plugin_registration(): + """The geomeTRIC plugin registers a geometry optimizer when installed.""" + assert "geometric" in algorithms.available("geometry_optimizer") + optimizer = algorithms.create("geometry_optimizer", "geometric") + assert isinstance(optimizer, algorithms.GeometryOptimizer) + assert optimizer.name() == "geometric" + + +def test_geometric_optimizer_settings(): + """The geomeTRIC optimizer exposes shared optimizer settings.""" + optimizer = algorithms.create("geometry_optimizer", "geometric") + settings = optimizer.settings() + + assert settings.get("transition_state") is False + assert settings.get("compute_hessian") is False + assert settings.get("max_iterations") == 300 + + +def test_geometric_optimizer_smoke_run(): + """Run a small geometry optimization through the optional geomeTRIC backend.""" + optimizer = algorithms.create("geometry_optimizer", "geometric") + optimizer.settings().set("max_iterations", 2) + derivative_ref = data.AlgorithmRef("nuclear_derivative_calculator", "finite_difference") + derivative_ref.get_settings().set("finite_difference_step", 1.0e-2) + optimizer.settings().set("derivative_calculator", derivative_ref) + + energy, structure, wavefunction, hessian = optimizer.run(_h2_structure(), 0, 1, "sto-3g") + + assert np.isfinite(energy) + assert structure.get_num_atoms() == 2 + assert structure.get_coordinates().shape == (2, 3) + assert wavefunction is not None + assert hessian is None diff --git a/python/tests/test_nuclear_derivative.py b/python/tests/test_nuclear_derivative.py index 139afa5cf..c8608ac2e 100644 --- a/python/tests/test_nuclear_derivative.py +++ b/python/tests/test_nuclear_derivative.py @@ -66,3 +66,11 @@ def test_qdk_nuclear_derivative_factory_registered(): assert isinstance(calculator, algorithms.QdkNuclearDerivativeCalculator) assert calculator.name() == "qdk" + + +def test_geometry_optimizer_factory_registered(): + """The geometry optimizer algorithm type is registered in the core registry.""" + assert algorithms.show_default("geometry_optimizer") == "geometric" + available = algorithms.available("geometry_optimizer") + assert isinstance(available, list) + assert not available or "geometric" in available From 32694c2767290916a2d7a8102c8d1740b3592092 Mon Sep 17 00:00:00 2001 From: Jan Unsleber Date: Mon, 8 Jun 2026 15:46:18 +0200 Subject: [PATCH 12/13] add tsopt --- .../plugins/geometric/__init__.py | 16 +++- .../plugins/geometric/geometry_optimizer.py | 83 +++++++++++++------ python/tests/test_geometric_plugin.py | 30 +++++-- 3 files changed, 94 insertions(+), 35 deletions(-) diff --git a/python/src/qdk_chemistry/plugins/geometric/__init__.py b/python/src/qdk_chemistry/plugins/geometric/__init__.py index 876d62ce1..bcc7a5149 100644 --- a/python/src/qdk_chemistry/plugins/geometric/__init__.py +++ b/python/src/qdk_chemistry/plugins/geometric/__init__.py @@ -30,7 +30,15 @@ def _register_algorithms(): """Register geomeTRIC-backed optimizer algorithms.""" Logger.trace_entering() from qdk_chemistry.algorithms import register # noqa: PLC0415 - from qdk_chemistry.plugins.geometric.geometry_optimizer import GeometricOptimizer # noqa: PLC0415 - - register(lambda: GeometricOptimizer()) - Logger.debug(f"geomeTRIC plugin loaded: [{GeometricOptimizer().type_name()}: {GeometricOptimizer().name()}].") + from qdk_chemistry.plugins.geometric.geometry_optimizer import ( # noqa: PLC0415 + GEOMETRIC_OPTIMIZER_ALGORITHMS, + GeometricOptimizer, + ) + + for algorithm in GEOMETRIC_OPTIMIZER_ALGORITHMS: + register(lambda algorithm=algorithm: GeometricOptimizer(algorithm=algorithm)) + register(lambda algorithm=algorithm: GeometricOptimizer(algorithm=algorithm, transition_state=True)) + + Logger.debug( + f"geomeTRIC plugin loaded: [{GeometricOptimizer().type_name()}: {', '.join(GEOMETRIC_OPTIMIZER_ALGORITHMS)}]." + ) diff --git a/python/src/qdk_chemistry/plugins/geometric/geometry_optimizer.py b/python/src/qdk_chemistry/plugins/geometric/geometry_optimizer.py index 4bd1c7a4f..3f880c2c5 100644 --- a/python/src/qdk_chemistry/plugins/geometric/geometry_optimizer.py +++ b/python/src/qdk_chemistry/plugins/geometric/geometry_optimizer.py @@ -17,13 +17,27 @@ from qdk_chemistry.data import AlgorithmRef, NuclearHessian, Settings, Structure from qdk_chemistry.utils import Logger -__all__ = ["GeometricOptimizer", "GeometricOptimizerSettings"] +__all__ = [ + "GEOMETRIC_OPTIMIZER_ALGORITHMS", + "GeometricOptimizer", + "GeometricOptimizerSettings", +] + + +GEOMETRIC_OPTIMIZER_ALGORITHMS = { + "tric": "tric", + "tric_p": "tric-p", + "dlc": "dlc", + "hdlc": "hdlc", + "prim": "prim", + "cartesian": "cart", +} class GeometricOptimizerSettings(Settings): """Settings for the geomeTRIC geometry optimizer.""" - def __init__(self): + def __init__(self, *, transition_state: bool = False, coordinate_system: str = "tric"): """Initialize geomeTRIC optimizer defaults.""" super().__init__() self._set_default( @@ -33,7 +47,14 @@ def __init__(self): "Nuclear derivative calculator used to evaluate energies and gradients.", ) self._set_default( - "transition_state", "bool", False, "Run transition-state optimization instead of minimization." + "transition_state", "bool", transition_state, "Run transition-state optimization instead of minimization." + ) + self._set_default( + "coordinate_system", + "string", + coordinate_system, + "geomeTRIC coordinate-system optimizer algorithm.", + limit=list(GEOMETRIC_OPTIMIZER_ALGORITHMS.values()), ) self._set_default("max_iterations", "int", 300, "Maximum number of geometry optimization steps.") self._set_default("convergence_energy", "double", 1.0e-6, "Energy convergence threshold.") @@ -53,8 +74,9 @@ def __init__( spin_multiplicity: int, seed_or_basis: Any, derivative_calculator: Any, + molecule: Any, ): - super().__init__() + super().__init__(molecule) self._structure = structure self._charge = charge self._spin_multiplicity = spin_multiplicity @@ -90,19 +112,28 @@ def calc_new(self, coordinates: np.ndarray, dirname: str) -> dict[str, np.ndarra class GeometricOptimizer(GeometryOptimizer): """Geometry optimizer implemented with the geomeTRIC Python library.""" - def __init__(self): + def __init__(self, *, algorithm: str = "tric", transition_state: bool = False): """Initialize the geomeTRIC optimizer.""" Logger.trace_entering() super().__init__() - self._settings = GeometricOptimizerSettings() + if algorithm not in GEOMETRIC_OPTIMIZER_ALGORITHMS: + raise ValueError(f"Unknown geomeTRIC optimizer algorithm: {algorithm}") + self._algorithm = algorithm + self._coordinate_system = GEOMETRIC_OPTIMIZER_ALGORITHMS[algorithm] + self._transition_state = transition_state + self._settings = GeometricOptimizerSettings( + transition_state=transition_state, + coordinate_system=self._coordinate_system, + ) def name(self) -> str: """Return the implementation name.""" - return "geometric" + mode = "tsopt" if self._transition_state else "geoopt" + return f"geometric_{mode}_{self._algorithm}" def aliases(self) -> list[str]: """Return accepted factory aliases.""" - return ["geometric", "geomeTRIC"] + return [self.name()] def _run_impl( self, structure: Structure, charge: int, spin_multiplicity: int, seed_or_basis: Any @@ -112,28 +143,21 @@ def _run_impl( from geometric.molecule import Molecule # noqa: PLC0415 from geometric.optimize import run_optimizer # noqa: PLC0415 - derivative_calculator = self._create_nested("derivative_calculator") - derivative_calculator.settings().set("compute_hessian", False) - engine = _QdkDerivativeEngine(structure, charge, spin_multiplicity, seed_or_basis, derivative_calculator) - molecule = Molecule() molecule.elem = structure.get_atomic_symbols() molecule.xyzs = [np.asarray(structure.get_coordinates(), dtype=float)] - engine.M = molecule - params = { - "customengine": engine, - "input": None, - "maxiter": self._settings["max_iterations"], - "convergence_energy": self._settings["convergence_energy"], - "convergence_grms": self._settings["convergence_gradient"], - "convergence_drms": self._settings["convergence_displacement"], - "transition": self._settings["transition_state"], - "verbose": self._settings["print_level"], - } + derivative_calculator = self._create_nested("derivative_calculator") + derivative_calculator.settings().set("compute_hessian", False) + engine = _QdkDerivativeEngine( + structure, charge, spin_multiplicity, seed_or_basis, derivative_calculator, molecule + ) + + params = self._geometric_options() + params.update({"customengine": engine, "input": None}) with TemporaryDirectory(prefix="qdk-chemistry-geometric-") as tmpdir: - result = run_optimizer(**params, dirname=tmpdir) + result = run_optimizer(**params, prefix=f"{tmpdir}/qdk-chemistry", dirname=tmpdir) optimized_coordinates = _extract_coordinates(result, engine) optimized_structure = engine.structure_from_coordinates(optimized_coordinates) @@ -148,6 +172,17 @@ def _run_impl( return final_energy, optimized_structure, wavefunction, hessian + def _geometric_options(self) -> dict[str, Any]: + return { + "transition": self._settings["transition_state"], + "coordsys": self._settings["coordinate_system"], + "maxiter": self._settings["max_iterations"], + "convergence_energy": self._settings["convergence_energy"], + "convergence_grms": self._settings["convergence_gradient"], + "convergence_drms": self._settings["convergence_displacement"], + "verbose": self._settings["print_level"], + } + def _extract_coordinates(result: Any, engine: _QdkDerivativeEngine) -> np.ndarray: """Extract final coordinates from geomeTRIC's return value.""" diff --git a/python/tests/test_geometric_plugin.py b/python/tests/test_geometric_plugin.py index 0b16a9e04..4b4e7abd8 100644 --- a/python/tests/test_geometric_plugin.py +++ b/python/tests/test_geometric_plugin.py @@ -10,6 +10,7 @@ from qdk_chemistry import algorithms, data from qdk_chemistry.plugins.geometric import QDK_CHEMISTRY_HAS_GEOMETRIC +from qdk_chemistry.plugins.geometric.geometry_optimizer import GEOMETRIC_OPTIMIZER_ALGORITHMS pytestmark = pytest.mark.skipif(not QDK_CHEMISTRY_HAS_GEOMETRIC, reason="geomeTRIC not available") @@ -20,28 +21,43 @@ def _h2_structure(): def test_geometric_plugin_registration(): """The geomeTRIC plugin registers a geometry optimizer when installed.""" - assert "geometric" in algorithms.available("geometry_optimizer") - optimizer = algorithms.create("geometry_optimizer", "geometric") + available = algorithms.available("geometry_optimizer") + for algorithm in GEOMETRIC_OPTIMIZER_ALGORITHMS: + assert f"geometric_geoopt_{algorithm}" in available + assert f"geometric_tsopt_{algorithm}" in available + + optimizer = algorithms.create("geometry_optimizer", "geometric_geoopt_tric") assert isinstance(optimizer, algorithms.GeometryOptimizer) - assert optimizer.name() == "geometric" + assert optimizer.name() == "geometric_geoopt_tric" + assert optimizer.settings().get("coordinate_system") == "tric" + + ts_optimizer = algorithms.create("geometry_optimizer", "geometric_tsopt_dlc") + assert isinstance(ts_optimizer, algorithms.GeometryOptimizer) + assert ts_optimizer.name() == "geometric_tsopt_dlc" + assert ts_optimizer.settings().get("coordinate_system") == "dlc" + assert ts_optimizer.settings().get("transition_state") is True def test_geometric_optimizer_settings(): """The geomeTRIC optimizer exposes shared optimizer settings.""" - optimizer = algorithms.create("geometry_optimizer", "geometric") + optimizer = algorithms.create("geometry_optimizer", "geometric_geoopt_cartesian") settings = optimizer.settings() assert settings.get("transition_state") is False + assert settings.get("coordinate_system") == "cart" assert settings.get("compute_hessian") is False assert settings.get("max_iterations") == 300 + assert settings.get("convergence_energy") == pytest.approx(1.0e-6) + assert settings.get("convergence_gradient") == pytest.approx(3.0e-4) + assert settings.get("convergence_displacement") == pytest.approx(1.2e-3) def test_geometric_optimizer_smoke_run(): """Run a small geometry optimization through the optional geomeTRIC backend.""" - optimizer = algorithms.create("geometry_optimizer", "geometric") - optimizer.settings().set("max_iterations", 2) + optimizer = algorithms.create("geometry_optimizer", "geometric_geoopt_tric") + optimizer.settings().set("max_iterations", 20) derivative_ref = data.AlgorithmRef("nuclear_derivative_calculator", "finite_difference") - derivative_ref.get_settings().set("finite_difference_step", 1.0e-2) + derivative_ref.set("finite_difference_step", 1.0e-2) optimizer.settings().set("derivative_calculator", derivative_ref) energy, structure, wavefunction, hessian = optimizer.run(_h2_structure(), 0, 1, "sto-3g") From c983a8d4f537ff913656199203f929d5db9b7c17 Mon Sep 17 00:00:00 2001 From: Jan Unsleber Date: Mon, 8 Jun 2026 17:26:10 +0200 Subject: [PATCH 13/13] fixes --- .../_static/examples/python/custom_plugin.py | 24 +++++++------------ .../plugins/geometric/__init__.py | 1 + .../plugins/geometric/geometry_optimizer.py | 5 +++- python/tests/test_geometric_plugin.py | 6 +++++ 4 files changed, 20 insertions(+), 16 deletions(-) diff --git a/docs/source/_static/examples/python/custom_plugin.py b/docs/source/_static/examples/python/custom_plugin.py index 7e961cd0f..54a0a8006 100644 --- a/docs/source/_static/examples/python/custom_plugin.py +++ b/docs/source/_static/examples/python/custom_plugin.py @@ -147,14 +147,7 @@ def __init__(self): ################################################################################ # start-cell-geometry-base-class -from qdk_chemistry.algorithms.base import Algorithm - - -class GeometryOptimizer(Algorithm): - """Abstract base class for geometry optimization algorithms.""" - - def type_name(self) -> str: - return "geometry_optimizer" +from qdk_chemistry.algorithms import GeometryOptimizer # end-cell-geometry-base-class @@ -193,14 +186,16 @@ def __init__(self): def name(self) -> str: return "bfgs" - def _run_impl(self, structure: Structure) -> Structure: + def _run_impl( + self, structure: Structure, charge: int, spin_multiplicity: int, seed + ): # max_steps = self.settings().get("max_steps") # threshold = self.settings().get("convergence_threshold") # BFGS optimization implementation # Placeholder for optimized structure optimized_structure = Structure() - return optimized_structure + return 0.0, optimized_structure, None, None # end-cell-geometry-implementations @@ -219,11 +214,13 @@ def __init__(self): def name(self) -> str: return "steepest_descent" - def _run_impl(self, structure: Structure) -> Structure: + def _run_impl( + self, structure: Structure, charge: int, spin_multiplicity: int, seed + ): # Steepest descent implementation # Placeholder for optimized structure optimized_structure = Structure() - return optimized_structure + return 0.0, optimized_structure, None, None # end-cell-steepest-descent @@ -233,9 +230,6 @@ def _run_impl(self, structure: Structure) -> Structure: # start-cell-geometry-registration import qdk_chemistry.algorithms as algorithms -# Register the factory -algorithms.registry.register_factory(GeometryOptimizerFactory()) - # Register implementations algorithms.register(lambda: BfgsOptimizer()) algorithms.register(lambda: SteepestDescentOptimizer()) diff --git a/python/src/qdk_chemistry/plugins/geometric/__init__.py b/python/src/qdk_chemistry/plugins/geometric/__init__.py index bcc7a5149..4055f37cb 100644 --- a/python/src/qdk_chemistry/plugins/geometric/__init__.py +++ b/python/src/qdk_chemistry/plugins/geometric/__init__.py @@ -35,6 +35,7 @@ def _register_algorithms(): GeometricOptimizer, ) + register(lambda: GeometricOptimizer(algorithm="tric", name="geometric")) for algorithm in GEOMETRIC_OPTIMIZER_ALGORITHMS: register(lambda algorithm=algorithm: GeometricOptimizer(algorithm=algorithm)) register(lambda algorithm=algorithm: GeometricOptimizer(algorithm=algorithm, transition_state=True)) diff --git a/python/src/qdk_chemistry/plugins/geometric/geometry_optimizer.py b/python/src/qdk_chemistry/plugins/geometric/geometry_optimizer.py index 3f880c2c5..c573832ec 100644 --- a/python/src/qdk_chemistry/plugins/geometric/geometry_optimizer.py +++ b/python/src/qdk_chemistry/plugins/geometric/geometry_optimizer.py @@ -112,7 +112,7 @@ def calc_new(self, coordinates: np.ndarray, dirname: str) -> dict[str, np.ndarra class GeometricOptimizer(GeometryOptimizer): """Geometry optimizer implemented with the geomeTRIC Python library.""" - def __init__(self, *, algorithm: str = "tric", transition_state: bool = False): + def __init__(self, *, algorithm: str = "tric", transition_state: bool = False, name: str | None = None): """Initialize the geomeTRIC optimizer.""" Logger.trace_entering() super().__init__() @@ -121,6 +121,7 @@ def __init__(self, *, algorithm: str = "tric", transition_state: bool = False): self._algorithm = algorithm self._coordinate_system = GEOMETRIC_OPTIMIZER_ALGORITHMS[algorithm] self._transition_state = transition_state + self._name = name self._settings = GeometricOptimizerSettings( transition_state=transition_state, coordinate_system=self._coordinate_system, @@ -128,6 +129,8 @@ def __init__(self, *, algorithm: str = "tric", transition_state: bool = False): def name(self) -> str: """Return the implementation name.""" + if self._name is not None: + return self._name mode = "tsopt" if self._transition_state else "geoopt" return f"geometric_{mode}_{self._algorithm}" diff --git a/python/tests/test_geometric_plugin.py b/python/tests/test_geometric_plugin.py index 4b4e7abd8..adee97d82 100644 --- a/python/tests/test_geometric_plugin.py +++ b/python/tests/test_geometric_plugin.py @@ -22,10 +22,16 @@ def _h2_structure(): def test_geometric_plugin_registration(): """The geomeTRIC plugin registers a geometry optimizer when installed.""" available = algorithms.available("geometry_optimizer") + assert "geometric" in available for algorithm in GEOMETRIC_OPTIMIZER_ALGORITHMS: assert f"geometric_geoopt_{algorithm}" in available assert f"geometric_tsopt_{algorithm}" in available + default_optimizer = algorithms.create("geometry_optimizer", "geometric") + assert isinstance(default_optimizer, algorithms.GeometryOptimizer) + assert default_optimizer.name() == "geometric" + assert default_optimizer.settings().get("coordinate_system") == "tric" + optimizer = algorithms.create("geometry_optimizer", "geometric_geoopt_tric") assert isinstance(optimizer, algorithms.GeometryOptimizer) assert optimizer.name() == "geometric_geoopt_tric"