diff --git a/cpp/include/qdk/chemistry.hpp b/cpp/include/qdk/chemistry.hpp index c253148c7..b3fe3ec80 100644 --- a/cpp/include/qdk/chemistry.hpp +++ b/cpp/include/qdk/chemistry.hpp @@ -4,15 +4,19 @@ #include #include +#include #include #include #include #include +#include #include #include #include #include #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/include/qdk/chemistry/algorithms/nuclear_derivative.hpp b/cpp/include/qdk/chemistry/algorithms/nuclear_derivative.hpp new file mode 100644 index 000000000..9238c6b16 --- /dev/null +++ b/cpp/include/qdk/chemistry/algorithms/nuclear_derivative.hpp @@ -0,0 +1,258 @@ +// 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-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, + 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 shared nuclear derivative defaults. + */ + NuclearDerivativeSettings() { + 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"), + "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."); + 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), + "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), + "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()}); + } +}; + +/** + * @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. + * + * 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 Construct a calculator with finite-difference settings. + */ + FiniteDifferenceNuclearDerivativeCalculator() { + _settings = std::make_unique(); + } + + /** + * @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 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..a6b376eb8 --- /dev/null +++ b/cpp/include/qdk/chemistry/data/nuclear_gradients.hpp @@ -0,0 +1,138 @@ +// 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 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& gradient_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 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 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..963901d93 --- /dev/null +++ b/cpp/include/qdk/chemistry/data/nuclear_hessian.hpp @@ -0,0 +1,132 @@ +// 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 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& hessian_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 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 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/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/CMakeLists.txt b/cpp/src/qdk/chemistry/algorithms/CMakeLists.txt index f28aa2d61..dc026562f 100644 --- a/cpp/src/qdk/chemistry/algorithms/CMakeLists.txt +++ b/cpp/src/qdk/chemistry/algorithms/CMakeLists.txt @@ -4,8 +4,10 @@ target_sources(chemistry PRIVATE hamiltonian.cpp localization.cpp mc.cpp + 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 b77316bbf..2f8e80d4d 100644 --- a/cpp/src/qdk/chemistry/algorithms/algorithm_defaults.cpp +++ b/cpp/src/qdk/chemistry/algorithms/algorithm_defaults.cpp @@ -5,10 +5,12 @@ #include #include #include +#include #include #include #include #include +#include #include #include #include @@ -40,7 +42,9 @@ 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) REGISTER_FACTORY_SETTINGS_INIT(StabilityCheckerFactory) 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/cpp/src/qdk/chemistry/algorithms/microsoft/scf.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf.cpp index 0c5ed92de..4e2d34ed3 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(), @@ -505,7 +505,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 && 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) { + 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..416b159a9 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,20 @@ 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 + * when available. + * + * 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, + int spin_multiplicity, BasisOrGuessType basis_or_guess) const; + protected: /** * @brief Perform an SCF calculation on the given molecular structure @@ -154,6 +178,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/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..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 @@ -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_ @@ -477,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 * @@ -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,216 @@ 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(); + } + + 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 + 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(); + 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; + + 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], 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], sp12_data, + sp34_data); + } + + 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.5 * dJ_coord; + if (need_K) { + dK_out[coord] -= 0.25 * dK_coord; + } + } + } + } + } + } + } } /** 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..fb74df7dd --- /dev/null +++ b/cpp/src/qdk/chemistry/algorithms/nuclear_derivative.cpp @@ -0,0 +1,579 @@ +// 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 +#include +#include + +#ifdef QDK_CHEMISTRY_ENABLE_MPI +#include +#endif + +#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); + scoped_level_ = minimum_level; + changed_ = true; + } + } + + ~ScopedLogLevel() { + if (changed_) { + utils::Logger::restore_global_level_if_unchanged(scoped_level_, + previous_level_); + } + } + + ScopedLogLevel(const ScopedLogLevel&) = delete; + ScopedLogLevel& operator=(const ScopedLogLevel&) = delete; + + private: + utils::LogLevel previous_level_; + utils::LogLevel scoped_level_ = utils::LogLevel::off; + bool changed_ = false; +}; + +struct EnergyEvaluation { + double energy = 0.0; + 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) { + throw std::invalid_argument("Structure must not be null"); + } + 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) { + 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); + +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 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"); + } + + 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(hartree_fock_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()); + 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 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.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 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.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); + 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 = + 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 (!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 " + "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)); + 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), *nuclear_gradient); + + return {scf_result.energy, gradients, 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 diff --git a/cpp/src/qdk/chemistry/data/CMakeLists.txt b/cpp/src/qdk/chemistry/data/CMakeLists.txt index 0269e03a6..4422b9abb 100644 --- a/cpp/src/qdk/chemistry/data/CMakeLists.txt +++ b/cpp/src/qdk/chemistry/data/CMakeLists.txt @@ -15,6 +15,8 @@ target_sources(chemistry PRIVATE majorana_mapping.cpp majorana_mapping_factories.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..6314fb649 --- /dev/null +++ b/cpp/src/qdk/chemistry/data/nuclear_gradients.cpp @@ -0,0 +1,219 @@ +// 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& gradient_values) + : structure_(std::move(structure)), values_(gradient_values) { + if (!_is_valid()) { + throw std::invalid_argument( + "NuclearGradients requires 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 { + 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); + matrix(atom, 2) = values_(3 * atom + 2); + } + return matrix; +} + +std::string NuclearGradients::get_summary() const { + std::ostringstream oss; + oss << "NuclearGradients(" << get_structure()->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 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..61227e742 --- /dev/null +++ b/cpp/src/qdk/chemistry/data/nuclear_hessian.cpp @@ -0,0 +1,208 @@ +// 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& 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"); + } +} + +const std::shared_ptr NuclearHessian::get_structure() const { + if (!structure_) { + throw std::runtime_error("No structure is associated with this Hessian"); + } + return structure_; +} + +std::string NuclearHessian::get_summary() const { + std::ostringstream oss; + oss << "NuclearHessian(" << get_structure()->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 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/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 new file mode 100644 index 000000000..f9b39817d --- /dev/null +++ b/cpp/tests/test_nuclear_derivative.cpp @@ -0,0 +1,397 @@ +// 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 +#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::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) { + 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))}; + } +}; + +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; + } +}; + +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); + 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()); +} + +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); + if (calculator->settings().has("finite_difference_step")) { + calculator->settings().set("finite_difference_step", 1.0e-2); + } + + return calculator->run(testing::create_lih_structure(), 0, 1, + std::string("sto-3g")); +} + +void expect_qdk_analytic_gradient_matches_numeric( + 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); + 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_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()); + 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), gradient_tolerance) + << "gradient component " << i; + } +} + +} // namespace + +TEST(NuclearDerivativeSettingsTest, SettingsHaveUserFacingDescriptions) { + NuclearDerivativeSettings settings; + for (const auto& key : + {"energy_calculator", "orbital_solver", "hamiltonian_constructor", + "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 = 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_structure()->get_num_atoms(), 3); + EXPECT_TRUE(gradients.get_values().isApprox(values)); + EXPECT_EQ(gradients.as_matrix().rows(), 3); + 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 = testing::create_water_structure(); + Eigen::MatrixXd matrix = Eigen::MatrixXd::Identity(9, 9); + + NuclearHessian hessian(structure, matrix); + + EXPECT_EQ(hessian.get_structure()->get_num_atoms(), 3); + 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, 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( + 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-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_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)); + + ASSERT_TRUE(wavefunction.has_value()); + 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(); + recorded_inactive_spaces.clear(); + [[maybe_unused]] ScopedFactoryRegistration< + MultiConfigurationCalculatorFactory> + mc_guard("_test_nuclear_derivative_recording_mc"); + MultiConfigurationCalculatorFactory::register_instance( + []() -> MultiConfigurationCalculatorFactory::return_type { + return std::make_unique(); + }); + + auto structure = testing::create_water_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()); + } +} + +TEST(NuclearDerivativeCalculatorTest, + MultiReferenceFiniteDifferenceLocalizesReferenceOrbitals) { + recorded_active_spaces.clear(); + recorded_inactive_spaces.clear(); + recorded_localized_alpha_spaces.clear(); + recorded_localized_beta_spaces.clear(); + [[maybe_unused]] ScopedFactoryRegistration< + MultiConfigurationCalculatorFactory> + mc_guard("_test_nuclear_derivative_recording_mc"); + MultiConfigurationCalculatorFactory::register_instance( + []() -> MultiConfigurationCalculatorFactory::return_type { + return std::make_unique(); + }); + [[maybe_unused]] ScopedFactoryRegistration localizer_guard( + "_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}); + } +} + +TEST(NuclearDerivativeCalculatorTest, + QdkAnalyticGradientMatchesNumericHybridGgaDftForLithiumHydride) { + expect_qdk_analytic_gradient_matches_numeric("b3lyp"); +} + +TEST(NuclearDerivativeCalculatorTest, + QdkAnalyticGradientMatchesNumericGgaDftForLithiumHydride) { + expect_qdk_analytic_gradient_matches_numeric("pbe"); +} + +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..d3a9a552e 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_lih_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 */ 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/CMakeLists.txt b/python/CMakeLists.txt index cf75bb488..9966a0411 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -132,6 +132,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/majorana_mapping.cpp src/pybind11/data/serialization.cpp @@ -143,6 +145,8 @@ 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/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/algorithms/nuclear_derivative.cpp b/python/src/pybind11/algorithms/nuclear_derivative.cpp new file mode 100644 index 000000000..d972a883a --- /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.)"); +} diff --git a/python/src/pybind11/data/nuclear_gradients.cpp b/python/src/pybind11/data/nuclear_gradients.cpp new file mode 100644 index 000000000..001315ecc --- /dev/null +++ b/python/src/pybind11/data/nuclear_gradients.cpp @@ -0,0 +1,123 @@ +// 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_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(); + }); +} diff --git a/python/src/pybind11/data/nuclear_hessian.cpp b/python/src/pybind11/data/nuclear_hessian.cpp new file mode 100644 index 000000000..fc1cb132a --- /dev/null +++ b/python/src/pybind11/data/nuclear_hessian.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 +#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_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(); }); +} diff --git a/python/src/pybind11/module.cpp b/python/src/pybind11/module.cpp index a50f5687a..9d2804649 100644 --- a/python/src/pybind11/module.cpp +++ b/python/src/pybind11/module.cpp @@ -22,6 +22,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); @@ -30,6 +32,8 @@ 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_geometry_optimization(py::module& m); void bind_active_space(py::module& m); void bind_constants(py::module& m); void bind_pmc(py::module& m); @@ -81,6 +85,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_majorana_mapping( data); // Depends on SparsePauliWord from pauli_operator @@ -91,6 +97,8 @@ PYBIND11_MODULE(_core, m) { bind_mcscf(algorithms); 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 91423195a..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, @@ -37,6 +38,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, @@ -64,10 +70,13 @@ "ControlledCircuitMapper", "DynamicalCorrelationCalculator", "EnergyEstimator", + "FiniteDifferenceNuclearDerivativeCalculator", + "GeometryOptimizer", "HamiltonianConstructor", "HamiltonianUnitaryBuilder", "MultiConfigurationCalculator", "MultiConfigurationScf", + "NuclearDerivativeCalculator", "OrbitalLocalizer", "PhaseEstimation", "ProjectedMultiConfigurationCalculator", @@ -79,6 +88,7 @@ "QdkMacisAsci", "QdkMacisCas", "QdkMacisPmc", + "QdkNuclearDerivativeCalculator", "QdkOccupationActiveSpaceSelector", "QdkPipekMezeyLocalizer", "QdkQubitMapper", 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/nuclear_derivative.py b/python/src/qdk_chemistry/algorithms/nuclear_derivative.py new file mode 100644 index 000000000..1dfbf3d34 --- /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 ( + FiniteDifferenceNuclearDerivativeCalculator, + NuclearDerivativeCalculator, + QdkNuclearDerivativeCalculator, +) + +__all__ = [ + "FiniteDifferenceNuclearDerivativeCalculator", + "NuclearDerivativeCalculator", + "QdkNuclearDerivativeCalculator", +] diff --git a/python/src/qdk_chemistry/algorithms/registry.py b/python/src/qdk_chemistry/algorithms/registry.py index afa591928..675896d83 100644 --- a/python/src/qdk_chemistry/algorithms/registry.py +++ b/python/src/qdk_chemistry/algorithms/registry.py @@ -473,10 +473,12 @@ def _register_cpp_factories(): from qdk_chemistry._core._algorithms import ( # noqa: PLC0415 ActiveSpaceSelectorFactory, DynamicalCorrelationCalculatorFactory, + GeometryOptimizerFactory, HamiltonianConstructorFactory, LocalizerFactory, MultiConfigurationCalculatorFactory, MultiConfigurationScfFactory, + NuclearDerivativeCalculatorFactory, ProjectedMultiConfigurationCalculatorFactory, ScfSolverFactory, StabilityCheckerFactory, @@ -487,8 +489,10 @@ def _register_cpp_factories(): register_factory(LocalizerFactory) register_factory(MultiConfigurationCalculatorFactory) register_factory(MultiConfigurationScfFactory) + 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/data/__init__.py b/python/src/qdk_chemistry/data/__init__.py index e60415053..e7dcceefc 100644 --- a/python/src/qdk_chemistry/data/__init__.py +++ b/python/src/qdk_chemistry/data/__init__.py @@ -30,6 +30,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. @@ -87,6 +89,8 @@ MajoranaMapping, ModelOrbitals, MP2Container, + NuclearGradients, + NuclearHessian, Orbitals, OrbitalType, PauliOperator, @@ -165,6 +169,8 @@ "MajoranaMapping", "MeasurementData", "ModelOrbitals", + "NuclearGradients", + "NuclearHessian", "OrbitalType", "Orbitals", "PauliOperator", 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..4055f37cb --- /dev/null +++ b/python/src/qdk_chemistry/plugins/geometric/__init__.py @@ -0,0 +1,45 @@ +"""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 ( # noqa: PLC0415 + GEOMETRIC_OPTIMIZER_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)) + + 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 new file mode 100644 index 000000000..c573832ec --- /dev/null +++ b/python/src/qdk_chemistry/plugins/geometric/geometry_optimizer.py @@ -0,0 +1,203 @@ +"""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__ = [ + "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, *, transition_state: bool = False, coordinate_system: str = "tric"): + """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", 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.") + 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, + molecule: Any, + ): + super().__init__(molecule) + 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, *, algorithm: str = "tric", transition_state: bool = False, name: str | None = None): + """Initialize the geomeTRIC optimizer.""" + Logger.trace_entering() + super().__init__() + 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._name = name + self._settings = GeometricOptimizerSettings( + transition_state=transition_state, + coordinate_system=self._coordinate_system, + ) + + 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}" + + def aliases(self) -> list[str]: + """Return accepted factory aliases.""" + return [self.name()] + + 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 + + molecule = Molecule() + molecule.elem = structure.get_atomic_symbols() + molecule.xyzs = [np.asarray(structure.get_coordinates(), dtype=float)] + + 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, prefix=f"{tmpdir}/qdk-chemistry", 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 _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.""" + 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..adee97d82 --- /dev/null +++ b/python/tests/test_geometric_plugin.py @@ -0,0 +1,75 @@ +"""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 +from qdk_chemistry.plugins.geometric.geometry_optimizer import GEOMETRIC_OPTIMIZER_ALGORITHMS + +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.""" + 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" + 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_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_geoopt_tric") + optimizer.settings().set("max_iterations", 20) + derivative_ref = data.AlgorithmRef("nuclear_derivative_calculator", "finite_difference") + 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") + + 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 new file mode 100644 index 000000000..c8608ac2e --- /dev/null +++ b/python/tests/test_nuclear_derivative.py @@ -0,0 +1,76 @@ +"""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(): + """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_structure().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(): + """Round-trip a nuclear Hessian through JSON.""" + structure = _h2_structure() + matrix = np.eye(6) + + hessian = data.NuclearHessian(structure, matrix) + + assert hessian.get_structure().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(): + """Create the default nuclear derivative calculator from the factory.""" + calculator = algorithms.create("nuclear_derivative_calculator") + + assert isinstance(calculator, algorithms.NuclearDerivativeCalculator) + assert calculator.name() == "finite_difference" + + +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) + 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