diff --git a/cpp/include/qdk/chemistry.hpp b/cpp/include/qdk/chemistry.hpp
index c253148c7..9392a6a98 100644
--- a/cpp/include/qdk/chemistry.hpp
+++ b/cpp/include/qdk/chemistry.hpp
@@ -8,11 +8,14 @@
 #include <qdk/chemistry/algorithms/localization.hpp>
 #include <qdk/chemistry/algorithms/mc.hpp>
 #include <qdk/chemistry/algorithms/mcscf.hpp>
+#include <qdk/chemistry/algorithms/nuclear_derivative.hpp>
 #include <qdk/chemistry/algorithms/pmc.hpp>
 #include <qdk/chemistry/algorithms/scf.hpp>
 #include <qdk/chemistry/algorithms/stability.hpp>
 #include <qdk/chemistry/data/hamiltonian.hpp>
 #include <qdk/chemistry/data/lattice_graph.hpp>
+#include <qdk/chemistry/data/nuclear_gradients.hpp>
+#include <qdk/chemistry/data/nuclear_hessian.hpp>
 #include <qdk/chemistry/data/orbitals.hpp>
 #include <qdk/chemistry/data/stability_result.hpp>
 #include <qdk/chemistry/data/structure.hpp>
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 <limits>
+#include <memory>
+#include <optional>
+#include <qdk/chemistry/algorithms/algorithm.hpp>
+#include <qdk/chemistry/algorithms/scf.hpp>
+#include <qdk/chemistry/data/nuclear_gradients.hpp>
+#include <qdk/chemistry/data/nuclear_hessian.hpp>
+#include <qdk/chemistry/data/settings.hpp>
+#include <qdk/chemistry/data/wavefunction.hpp>
+#include <string>
+#include <tuple>
+#include <variant>
+#include <vector>
+
+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<data::Orbitals>,
+                 std::shared_ptr<data::BasisSet>,
+                 std::shared_ptr<data::Wavefunction>, std::string>;
+
+/**
+ * @brief Energy, gradients, optional Hessian, and optional wavefunction.
+ */
+using NuclearDerivativeResult =
+    std::tuple<double, std::shared_ptr<data::NuclearGradients>,
+               std::optional<std::shared_ptr<data::NuclearHessian>>,
+               std::optional<std::shared_ptr<data::Wavefunction>>>;
+
+/**
+ * @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<int64_t>(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<int64_t>{0, std::numeric_limits<int64_t>::max()});
+    set_default(
+        "n_active_beta_electrons", static_cast<int64_t>(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<int64_t>{0, std::numeric_limits<int64_t>::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<double>{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<NuclearDerivativeCalculator, NuclearDerivativeResult,
+                       std::shared_ptr<data::Structure>, int, int,
+                       NuclearDerivativeSeedType> {
+ public:
+  /**
+   * @brief Construct a calculator with nuclear derivative settings.
+   */
+  NuclearDerivativeCalculator() {
+    _settings = std::make_unique<NuclearDerivativeSettings>();
+  }
+  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<data::Structure> structure, int charge,
+      int spin_multiplicity, NuclearDerivativeSeedType seed) const = 0;
+};
+
+/**
+ * @brief Factory for nuclear derivative calculator implementations.
+ */
+struct NuclearDerivativeCalculatorFactory
+    : public AlgorithmFactory<NuclearDerivativeCalculator,
+                              NuclearDerivativeCalculatorFactory> {
+  /**
+   * @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<FiniteDifferenceNuclearDerivativeSettings>();
+  }
+
+  /**
+   * @brief Return the implementation name.
+   */
+  std::string name() const final { return "finite_difference"; }
+
+  /**
+   * @brief Return accepted factory aliases for this implementation.
+   */
+  std::vector<std::string> aliases() const final {
+    return {"finite_difference", "numeric"};
+  }
+
+ protected:
+  /**
+   * @brief Compute finite-difference nuclear derivatives.
+   */
+  NuclearDerivativeResult _run_impl(
+      std::shared_ptr<data::Structure> 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<std::string> aliases() const final {
+    return {"qdk", "analytical_gradient"};
+  }
+
+ protected:
+  /**
+   * @brief Compute analytic nuclear gradients with the internal SCF backend.
+   */
+  NuclearDerivativeResult _run_impl(
+      std::shared_ptr<data::Structure> 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..2c3b1507d
--- /dev/null
+++ b/cpp/include/qdk/chemistry/data/nuclear_gradients.hpp
@@ -0,0 +1,139 @@
+// 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 <Eigen/Dense>
+#include <memory>
+#include <nlohmann/json.hpp>
+#include <qdk/chemistry/data/data_class.hpp>
+#include <qdk/chemistry/data/structure.hpp>
+#include <qdk/chemistry/utils/string_utils.hpp>
+#include <string>
+
+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<NuclearGradients> {
+ 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> structure,
+                   const Eigen::VectorXd& gradient_values);
+
+  /**
+   * @brief Get the molecular structure associated with these gradients.
+   */
+  const std::shared_ptr<Structure> 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<NuclearGradients> from_file(
+      const std::string& filename, const std::string& type);
+
+  /**
+   * @brief Load gradients from a JSON file.
+   */
+  static std::shared_ptr<NuclearGradients> from_json_file(
+      const std::string& filename);
+
+  /**
+   * @brief Load gradients from JSON.
+   */
+  static std::shared_ptr<NuclearGradients> from_json(const nlohmann::json& j);
+
+  /**
+   * @brief Load gradients from an HDF5 file.
+   */
+  static std::shared_ptr<NuclearGradients> from_hdf5_file(
+      const std::string& filename);
+
+  /**
+   * @brief Load gradients from an HDF5 group.
+   */
+  static std::shared_ptr<NuclearGradients> 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<NuclearGradients> _from_json_file(
+      const std::string& filename);
+  static std::shared_ptr<NuclearGradients> _from_hdf5_file(
+      const std::string& filename);
+  void hash_update(qdk::chemistry::utils::HashContext& ctx) const override;
+
+  std::shared_ptr<Structure> structure_;
+  Eigen::VectorXd values_;
+};
+
+static_assert(DataClassCompliant<NuclearGradients>,
+              "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..524e82429
--- /dev/null
+++ b/cpp/include/qdk/chemistry/data/nuclear_hessian.hpp
@@ -0,0 +1,133 @@
+// 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 <Eigen/Dense>
+#include <memory>
+#include <nlohmann/json.hpp>
+#include <qdk/chemistry/data/data_class.hpp>
+#include <qdk/chemistry/data/structure.hpp>
+#include <qdk/chemistry/utils/string_utils.hpp>
+#include <string>
+
+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<NuclearHessian> {
+ 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> structure,
+                 const Eigen::MatrixXd& hessian_matrix);
+
+  /**
+   * @brief Get the molecular structure associated with this Hessian.
+   */
+  const std::shared_ptr<Structure> 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<NuclearHessian> from_file(const std::string& filename,
+                                                   const std::string& type);
+
+  /**
+   * @brief Load the Hessian from a JSON file.
+   */
+  static std::shared_ptr<NuclearHessian> from_json_file(
+      const std::string& filename);
+
+  /**
+   * @brief Load the Hessian from JSON.
+   */
+  static std::shared_ptr<NuclearHessian> from_json(const nlohmann::json& j);
+
+  /**
+   * @brief Load the Hessian from an HDF5 file.
+   */
+  static std::shared_ptr<NuclearHessian> from_hdf5_file(
+      const std::string& filename);
+
+  /**
+   * @brief Load the Hessian from an HDF5 group.
+   */
+  static std::shared_ptr<NuclearHessian> 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<NuclearHessian> _from_json_file(
+      const std::string& filename);
+  static std::shared_ptr<NuclearHessian> _from_hdf5_file(
+      const std::string& filename);
+  void hash_update(qdk::chemistry::utils::HashContext& ctx) const override;
+
+  std::shared_ptr<Structure> structure_;
+  Eigen::MatrixXd matrix_;
+};
+
+static_assert(DataClassCompliant<NuclearHessian>,
+              "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 1ecf51e09..79470d3a0 100644
--- a/cpp/include/qdk/chemistry/data/settings.hpp
+++ b/cpp/include/qdk/chemistry/data/settings.hpp
@@ -17,6 +17,7 @@
 #include <optional>
 #include <qdk/chemistry/data/data_class.hpp>
 #include <qdk/chemistry/utils/string_utils.hpp>
+#include <set>
 #include <sstream>
 #include <stdexcept>
 #include <string>
@@ -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<std::string, std::string> descriptions_;
   std::map<std::string, Constraint> limits_;
   std::map<std::string, bool> documented_;
+  std::set<std::string> algorithm_ref_type_change_allowed_;
 
   void hash_update(qdk::chemistry::utils::HashContext& ctx) const override;
 
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..c06f2abae 100644
--- a/cpp/src/qdk/chemistry/algorithms/CMakeLists.txt
+++ b/cpp/src/qdk/chemistry/algorithms/CMakeLists.txt
@@ -4,6 +4,7 @@ target_sources(chemistry PRIVATE
     hamiltonian.cpp
     localization.cpp
     mc.cpp
+    nuclear_derivative.cpp
     pmc.cpp
     dynamical_correlation_calculator.cpp
     scf.cpp
diff --git a/cpp/src/qdk/chemistry/algorithms/algorithm_defaults.cpp b/cpp/src/qdk/chemistry/algorithms/algorithm_defaults.cpp
index b77316bbf..cbc009ac6 100644
--- a/cpp/src/qdk/chemistry/algorithms/algorithm_defaults.cpp
+++ b/cpp/src/qdk/chemistry/algorithms/algorithm_defaults.cpp
@@ -9,6 +9,7 @@
 #include <qdk/chemistry/algorithms/localization.hpp>
 #include <qdk/chemistry/algorithms/mc.hpp>
 #include <qdk/chemistry/algorithms/mcscf.hpp>
+#include <qdk/chemistry/algorithms/nuclear_derivative.hpp>
 #include <qdk/chemistry/algorithms/pmc.hpp>
 #include <qdk/chemistry/algorithms/scf.hpp>
 #include <qdk/chemistry/algorithms/stability.hpp>
@@ -41,6 +42,7 @@ std::shared_ptr<data::Settings> resolve_algorithm_defaults(
   REGISTER_FACTORY_SETTINGS_INIT(ProjectedMultiConfigurationCalculatorFactory)
   REGISTER_FACTORY_SETTINGS_INIT(DynamicalCorrelationCalculatorFactory)
   REGISTER_FACTORY_SETTINGS_INIT(MultiConfigurationScfFactory)
+  REGISTER_FACTORY_SETTINGS_INIT(NuclearDerivativeCalculatorFactory)
   REGISTER_FACTORY_SETTINGS_INIT(LocalizerFactory)
   REGISTER_FACTORY_SETTINGS_INIT(StabilityCheckerFactory)
 
diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/scf.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/scf.cpp
index 952c570fc..64ef38786 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<int, int> calculate_electron_counts(int nuclear_charge, int charge,
   return {n_alpha, n_beta};
 }
 
-std::pair<double, std::shared_ptr<data::Wavefunction>> ScfSolver::_run_impl(
+ScfCalculationResult ScfSolver::_run_with_options(
     std::shared_ptr<data::Structure> 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<double, std::shared_ptr<data::Wavefunction>> ScfSolver::_run_impl(
   // Create SCFConfig
   auto ms_scf_config = std::make_unique<qcs::SCFConfig>();
   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(),
@@ -495,7 +495,38 @@ std::pair<double, std::shared_ptr<data::Wavefunction>> ScfSolver::_run_impl(
   // Return total energy
   double total_energy = context.result.scf_total_energy;
 
-  return std::make_pair(total_energy, std::make_shared<data::Wavefunction>(
-                                          std::move(wavefunction)));
+  std::optional<Eigen::VectorXd> 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<const Eigen::VectorXd>(
+        raw_gradient.data(), static_cast<Eigen::Index>(raw_gradient.size()));
+  }
+
+  auto wavefunction_ptr =
+      std::make_shared<data::Wavefunction>(std::move(wavefunction));
+
+  return {total_energy, wavefunction_ptr, gradient};
+}
+
+ScfCalculationResult ScfSolver::run_with_analytic_gradient(
+    std::shared_ptr<data::Structure> 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<double, std::shared_ptr<data::Wavefunction>> ScfSolver::_run_impl(
+    std::shared_ptr<data::Structure> 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 <Eigen/Dense>
 #include <qdk/chemistry/algorithms/scf.hpp>
 
 namespace qdk::chemistry::algorithms::microsoft {
 
+/**
+ * @brief Internal SCF result including optional analytic nuclear gradients.
+ */
+struct ScfCalculationResult {
+  double energy;
+  std::shared_ptr<data::Wavefunction> wavefunction;
+  std::optional<Eigen::VectorXd> 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<data::Structure> 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<double, std::shared_ptr<data::Wavefunction>> _run_impl(
       std::shared_ptr<data::Structure> structure, int charge,
       int spin_multiplicity, BasisOrGuessType basis_or_guess) const override;
+
+ private:
+  ScfCalculationResult _run_with_options(
+      std::shared_ptr<data::Structure> 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<size_t>
       shell2bf_;  ///< Mapping from shell index to first atomic orbital
+  std::vector<int> 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<int>(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<double> dJ_buffer;
+    std::vector<double> 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<double> 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<std::vector<std::shared_ptr<::libint2::ShellPair>>>
+        shell_pair_data(nsh, std::vector<std::shared_ptr<::libint2::ShellPair>>(
+                                 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<size_t>(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<size_t>(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..4b08e1657
--- /dev/null
+++ b/cpp/src/qdk/chemistry/algorithms/nuclear_derivative.cpp
@@ -0,0 +1,580 @@
+// Copyright (c) Microsoft Corporation. All rights reserved.
+// Licensed under the MIT License. See LICENSE.txt in the project root for
+// license information.
+
+#include <qdk/chemistry/algorithms/active_space.hpp>
+#include <qdk/chemistry/algorithms/hamiltonian.hpp>
+#include <qdk/chemistry/algorithms/localization.hpp>
+#include <qdk/chemistry/algorithms/mc.hpp>
+#include <qdk/chemistry/algorithms/mcscf.hpp>
+#include <qdk/chemistry/algorithms/nuclear_derivative.hpp>
+#include <qdk/chemistry/algorithms/scf.hpp>
+#include <qdk/chemistry/data/wavefunction_containers/state_vector.hpp>
+#include <qdk/chemistry/utils/logger.hpp>
+#include <stdexcept>
+
+#ifdef QDK_CHEMISTRY_ENABLE_MPI
+#include <mpi.h>
+#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<int>(previous_level_) < static_cast<int>(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<std::shared_ptr<data::Wavefunction>> wavefunction;
+};
+
+struct ReferenceOrbitals {
+  std::shared_ptr<data::Orbitals> orbitals;
+  std::shared_ptr<data::Wavefunction> wavefunction;
+};
+
+std::shared_ptr<data::Structure> copy_structure(
+    const std::shared_ptr<data::Structure>& structure) {
+  if (!structure) {
+    throw std::invalid_argument("Structure must not be null");
+  }
+  return std::make_shared<data::Structure>(*structure);
+}
+
+std::optional<Eigen::VectorXd> broadcast_gradient_from_root(
+    const std::optional<Eigen::VectorXd>& 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<int>(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<data::Structure> displace_structure(
+    const std::shared_ptr<data::Structure>& 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<data::Structure>(
+      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<std::string>(seed)) {
+    return std::get<std::string>(seed);
+  }
+  if (std::holds_alternative<std::shared_ptr<data::BasisSet>>(seed)) {
+    auto basis = std::get<std::shared_ptr<data::BasisSet>>(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<data::Orbitals> orbitals;
+  if (std::holds_alternative<std::shared_ptr<data::Orbitals>>(seed)) {
+    orbitals = std::get<std::shared_ptr<data::Orbitals>>(seed);
+  } else if (std::holds_alternative<std::shared_ptr<data::Wavefunction>>(
+                 seed)) {
+    auto wavefunction = std::get<std::shared_ptr<data::Wavefunction>>(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>
+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<data::Orbitals> seed_to_orbitals(
+    const NuclearDerivativeSeedType& seed) {
+  std::shared_ptr<data::Orbitals> orbitals;
+  if (std::holds_alternative<std::shared_ptr<data::Orbitals>>(seed)) {
+    orbitals = std::get<std::shared_ptr<data::Orbitals>>(seed);
+    if (!orbitals) {
+      throw std::invalid_argument("Orbital seed must not be null");
+    }
+  } else if (std::holds_alternative<std::shared_ptr<data::Wavefunction>>(
+                 seed)) {
+    auto wavefunction = std::get<std::shared_ptr<data::Wavefunction>>(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<data::Wavefunction> seed_to_wavefunction(
+    const NuclearDerivativeSeedType& seed) {
+  if (!std::holds_alternative<std::shared_ptr<data::Wavefunction>>(seed)) {
+    return nullptr;
+  }
+  auto wavefunction = std::get<std::shared_ptr<data::Wavefunction>>(seed);
+  if (!wavefunction) {
+    throw std::invalid_argument("Wavefunction seed must not be null");
+  }
+  return wavefunction;
+}
+
+std::shared_ptr<data::Orbitals> copy_active_space_metadata(
+    const std::shared_ptr<data::Orbitals>& orbitals,
+    const std::shared_ptr<data::Orbitals>& 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<Eigen::MatrixXd> ao_overlap;
+  if (orbitals->has_overlap_matrix()) {
+    ao_overlap = orbitals->get_overlap_matrix();
+  }
+  std::shared_ptr<data::BasisSet> basis_set = orbitals->get_basis_set();
+
+  if (orbitals->is_restricted()) {
+    std::optional<Eigen::VectorXd> energies;
+    if (orbitals->has_energies()) {
+      energies = orbitals->get_energies().first;
+    }
+    return std::make_shared<data::Orbitals>(
+        orbitals->get_coefficients().first, energies, ao_overlap, basis_set,
+        std::make_tuple(
+            std::vector<size_t>(active_a.begin(), active_a.end()),
+            std::vector<size_t>(inactive_a.begin(), inactive_a.end())));
+  }
+
+  std::optional<Eigen::VectorXd> energies_a;
+  std::optional<Eigen::VectorXd> 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<data::Orbitals>(
+      orbitals->get_coefficients().first, orbitals->get_coefficients().second,
+      energies_a, energies_b, ao_overlap, basis_set,
+      std::make_tuple(
+          std::vector<size_t>(active_a.begin(), active_a.end()),
+          std::vector<size_t>(active_b.begin(), active_b.end()),
+          std::vector<size_t>(inactive_a.begin(), inactive_a.end()),
+          std::vector<size_t>(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<data::Wavefunction> wavefunction_from_orbitals(
+    std::shared_ptr<data::Orbitals> 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::from_spin_half_string(
+      hartree_fock_determinant_string(active_a.size(), n_active_alpha,
+                                      n_active_beta));
+  auto container =
+      std::make_unique<data::StateVectorContainer>(determinant, orbitals);
+  return std::make_shared<data::Wavefunction>(std::move(container));
+}
+
+ReferenceOrbitals localize_reference_orbitals(const data::Settings& settings,
+                                              ReferenceOrbitals reference) {
+  if (!settings.get<bool>("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<LocalizerFactory>(
+      settings.get<data::AlgorithmRef>("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<data::Structure> 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<ScfSolverFactory>(
+      settings.get<data::AlgorithmRef>("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<bool>("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>
+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<int64_t>(key);
+  if (value <= 0) {
+    throw std::invalid_argument(
+        key + " must be set to a positive value for this energy calculator");
+  }
+  return static_cast<unsigned int>(value);
+}
+
+EnergyEvaluation evaluate_energy(const data::Settings& settings,
+                                 std::shared_ptr<data::Structure> structure,
+                                 int charge, int spin_multiplicity,
+                                 const NuclearDerivativeSeedType& seed,
+                                 bool allow_wavefunction_seed) {
+  const auto ref = settings.get<data::AlgorithmRef>("energy_calculator");
+  const auto& algorithm_type = ref.get_algorithm_type();
+
+  if (algorithm_type == ScfSolverFactory::algorithm_type_name()) {
+    auto solver = create_from_ref<ScfSolverFactory>(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<MultiConfigurationScfFactory>(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<HamiltonianConstructorFactory>(
+            settings.get<data::AlgorithmRef>("hamiltonian_constructor"));
+    auto hamiltonian = hamiltonian_constructor->run(reference.orbitals);
+    auto calculator = create_from_ref<MultiConfigurationCalculatorFactory>(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<data::Structure> 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<double>("finite_difference_step");
+  const bool compute_hessian = _settings->get<bool>("compute_hessian");
+  const auto dimension =
+      static_cast<Eigen::Index>(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<double> plus_energies(static_cast<size_t>(dimension));
+  std::vector<double> minus_energies(static_cast<size_t>(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<std::shared_ptr<data::NuclearHessian>> 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<bool>("symmetrize_hessian")) {
+      hessian_matrix = 0.5 * (hessian_matrix + hessian_matrix.transpose());
+    }
+    hessian = std::make_shared<data::NuclearHessian>(copy_structure(structure),
+                                                     hessian_matrix);
+  }
+
+  return {central.energy,
+          std::make_shared<data::NuclearGradients>(copy_structure(structure),
+                                                   gradients),
+          hessian, central.wavefunction};
+}
+
+std::unique_ptr<NuclearDerivativeCalculator>
+make_finite_difference_nuclear_derivative_calculator() {
+  return std::make_unique<FiniteDifferenceNuclearDerivativeCalculator>();
+}
+
+std::unique_ptr<NuclearDerivativeCalculator>
+make_qdk_nuclear_derivative_calculator() {
+  return std::make_unique<QdkNuclearDerivativeCalculator>();
+}
+
+NuclearDerivativeResult QdkNuclearDerivativeCalculator::_run_impl(
+    std::shared_ptr<data::Structure> 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<bool>("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<data::AlgorithmRef>("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<data::NuclearGradients>(
+      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 5624422a6..1c6d6b67d 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..44dad1dd6
--- /dev/null
+++ b/cpp/src/qdk/chemistry/data/nuclear_gradients.cpp
@@ -0,0 +1,226 @@
+// Copyright (c) Microsoft Corporation. All rights reserved.
+// Licensed under the MIT License. See LICENSE.txt in the project root for
+// license information.
+
+#include <H5Cpp.h>
+
+#include <fstream>
+#include <qdk/chemistry/data/nuclear_gradients.hpp>
+#include <qdk/chemistry/utils/logger.hpp>
+#include <sstream>
+#include <stdexcept>
+
+#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> 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<Structure> 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<Eigen::Index>(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> 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> 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> NuclearGradients::from_json(
+    const nlohmann::json& j) {
+  if (j.contains("serialization_version")) {
+    validate_serialization_version(
+        SERIALIZATION_VERSION, j["serialization_version"].get<std::string>());
+  }
+  if (j.contains("type") &&
+      j["type"].get<std::string>() != "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<NuclearGradients>(
+      Structure::from_json(j["structure"]), json_to_vector(j["values"]));
+}
+
+std::shared_ptr<NuclearGradients> 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> 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<NuclearGradients>(
+      Structure::from_hdf5(structure_group),
+      load_vector_from_group(group, "values"));
+}
+
+bool NuclearGradients::_is_valid() const {
+  return structure_ != nullptr &&
+         values_.size() ==
+             static_cast<Eigen::Index>(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> 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> 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);
+}
+
+void NuclearGradients::hash_update(
+    qdk::chemistry::utils::HashContext& ctx) const {
+  hash_value(ctx, get_data_type_name());
+  hash_value(ctx, get_structure()->content_hash());
+  hash_value(ctx, values_);
+}
+
+}  // 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..f56cbd260
--- /dev/null
+++ b/cpp/src/qdk/chemistry/data/nuclear_hessian.cpp
@@ -0,0 +1,215 @@
+// Copyright (c) Microsoft Corporation. All rights reserved.
+// Licensed under the MIT License. See LICENSE.txt in the project root for
+// license information.
+
+#include <H5Cpp.h>
+
+#include <fstream>
+#include <qdk/chemistry/data/nuclear_hessian.hpp>
+#include <qdk/chemistry/utils/logger.hpp>
+#include <sstream>
+#include <stdexcept>
+
+#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> 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<Structure> 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> 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> 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> NuclearHessian::from_json(
+    const nlohmann::json& j) {
+  if (j.contains("serialization_version")) {
+    validate_serialization_version(
+        SERIALIZATION_VERSION, j["serialization_version"].get<std::string>());
+  }
+  if (j.contains("type") && j["type"].get<std::string>() != "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<NuclearHessian>(Structure::from_json(j["structure"]),
+                                          json_to_matrix(j["matrix"]));
+}
+
+std::shared_ptr<NuclearHessian> 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> 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<NuclearHessian>(
+      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<Eigen::Index>(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> 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> 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);
+}
+
+void NuclearHessian::hash_update(
+    qdk::chemistry::utils::HashContext& ctx) const {
+  hash_value(ctx, get_data_type_name());
+  hash_value(ctx, get_structure()->content_hash());
+  hash_value(ctx, matrix_);
+}
+
+}  // namespace qdk::chemistry::data
diff --git a/cpp/src/qdk/chemistry/data/settings.cpp b/cpp/src/qdk/chemistry/data/settings.cpp
index 568d3dc32..773682b07 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<AlgorithmRef>(value)) {
     const auto& existing = std::get<AlgorithmRef>(settings_[key]);
     const auto& incoming = std::get<AlgorithmRef>(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<spdlog::logger> logger;
+  std::lock_guard<std::mutex> lock(g_level_mutex);
+  set_global_level_locked(spdlog_level);
+}
 
-  // Update both our tracked level and spdlog's global level
-  {
-    std::lock_guard<std::mutex> 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<std::mutex> 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..4f2b40c12
--- /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 <gtest/gtest.h>
+
+#include <cmath>
+#include <memory>
+#include <qdk/chemistry/algorithms/localization.hpp>
+#include <qdk/chemistry/algorithms/mc.hpp>
+#include <qdk/chemistry/algorithms/nuclear_derivative.hpp>
+#include <qdk/chemistry/algorithms/scf.hpp>
+#include <qdk/chemistry/data/nuclear_gradients.hpp>
+#include <qdk/chemistry/data/nuclear_hessian.hpp>
+#include <qdk/chemistry/data/wavefunction_containers/state_vector.hpp>
+#include <stdexcept>
+#include <utility>
+#include <vector>
+
+#include "ut_common.hpp"
+
+using namespace qdk::chemistry::algorithms;
+using namespace qdk::chemistry::data;
+
+namespace {
+
+std::vector<std::vector<size_t>> recorded_active_spaces;
+std::vector<std::vector<size_t>> recorded_inactive_spaces;
+std::vector<std::vector<size_t>> recorded_localized_alpha_spaces;
+std::vector<std::vector<size_t>> 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<double, std::shared_ptr<Wavefunction>> _run_impl(
+      std::shared_ptr<Hamiltonian> 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::from_spin_half_string(determinant_string(
+        orbitals->get_active_space_indices().first.size(), n_alpha, n_beta));
+    auto container =
+        std::make_unique<StateVectorContainer>(determinant, orbitals);
+    return {0.0, std::make_shared<Wavefunction>(std::move(container))};
+  }
+};
+
+class RecordingLocalizer : public Localizer {
+ public:
+  std::string name() const override {
+    return "_test_nuclear_derivative_recording_localizer";
+  }
+
+ protected:
+  std::shared_ptr<Wavefunction> _run_impl(
+      std::shared_ptr<Wavefunction> wavefunction,
+      const std::vector<size_t>& loc_indices_a,
+      const std::vector<size_t>& 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 <typename Factory>
+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<Structure>& left,
+                           const std::shared_ptr<Structure>& 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<RecordingMultiConfigurationCalculator>();
+      });
+
+  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<size_t>{0}, std::vector<size_t>{});
+
+  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<size_t>{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<RecordingMultiConfigurationCalculator>();
+      });
+  [[maybe_unused]] ScopedFactoryRegistration<LocalizerFactory> localizer_guard(
+      "_test_nuclear_derivative_recording_localizer");
+  LocalizerFactory::register_instance([]() -> LocalizerFactory::return_type {
+    return std::make_unique<RecordingLocalizer>();
+  });
+
+  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<size_t>{0}, std::vector<size_t>{});
+
+  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<size_t>{0});
+  }
+  for (const auto& localized_space : recorded_localized_beta_spaces) {
+    EXPECT_EQ(localized_space, std::vector<size_t>{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 ce32823fe..963446945 100644
--- a/cpp/tests/ut_common.hpp
+++ b/cpp/tests/ut_common.hpp
@@ -268,6 +268,24 @@ inline std::shared_ptr<Structure> create_li_structure() {
   return std::make_shared<Structure>(coords, elements);
 }
 
+/**
+ * @brief Creates a LiH structure
+ */
+inline std::shared_ptr<Structure> create_lih_structure() {
+  std::vector<Eigen::Vector3d> 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<Element> elements = {qdk::chemistry::data::Element::Li,
+                                   qdk::chemistry::data::Element::H};
+
+  return std::make_shared<Structure>(coords, elements);
+}
+
 /**
  * @brief Creates an O2 structure
  */
diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt
index 3ac21adcd..a73061800 100644
--- a/python/CMakeLists.txt
+++ b/python/CMakeLists.txt
@@ -133,6 +133,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
@@ -144,6 +146,7 @@ pybind11_add_module(_core
     src/pybind11/algorithms/mcscf.cpp
     src/pybind11/algorithms/hamiltonian.cpp
     src/pybind11/algorithms/scf.cpp
+    src/pybind11/algorithms/nuclear_derivative.cpp
     src/pybind11/algorithms/active_space.cpp
     src/pybind11/algorithms/dynamical_correlation_calculator.cpp
     src/pybind11/algorithms/davidson_solver.cpp
diff --git a/python/src/pybind11/algorithms/nuclear_derivative.cpp b/python/src/pybind11/algorithms/nuclear_derivative.cpp
new file mode 100644
index 000000000..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 <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+
+#include <qdk/chemistry.hpp>
+
+#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<std::string> aliases() const override {
+    PYBIND11_OVERRIDE(std::vector<std::string>, NuclearDerivativeCalculator,
+                      aliases);
+  }
+
+  void replace_settings(
+      std::unique_ptr<qdk::chemistry::data::Settings> new_settings) {
+    this->_settings = std::move(new_settings);
+  }
+
+ protected:
+  NuclearDerivativeResult _run_impl(
+      std::shared_ptr<Structure> 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_<NuclearDerivativeCalculator, NuclearDerivativeCalculatorBase,
+             py::smart_holder>
+      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> 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<qdk::chemistry::data::Settings> 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.algorithms.NuclearDerivativeCalculator name='" +
+           self.name() + "'>";
+  });
+
+  qdk::chemistry::python::bind_algorithm_factory<
+      NuclearDerivativeCalculatorFactory, NuclearDerivativeCalculator,
+      NuclearDerivativeCalculatorBase>(m, "NuclearDerivativeCalculatorFactory");
+  qdk::chemistry::python::bind_create_nested(calculator);
+
+  py::class_<FiniteDifferenceNuclearDerivativeCalculator,
+             NuclearDerivativeCalculator, py::smart_holder>(
+      m, "FiniteDifferenceNuclearDerivativeCalculator",
+      R"(Numeric nuclear derivative calculator using central finite differences.)")
+      .def(py::init<>(),
+           R"(Create a finite-difference derivative calculator.)");
+
+  py::class_<QdkNuclearDerivativeCalculator, NuclearDerivativeCalculator,
+             py::smart_holder>(
+      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 <pybind11/eigen.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+
+#include <nlohmann/json.hpp>
+#include <qdk/chemistry/data/nuclear_gradients.hpp>
+#include <qdk/chemistry/utils/string_utils.hpp>
+
+#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<NuclearGradients> 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<NuclearGradients> 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<NuclearGradients> 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_<NuclearGradients, DataClass, py::smart_holder>(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<std::shared_ptr<Structure>, 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 <pybind11/eigen.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+
+#include <nlohmann/json.hpp>
+#include <qdk/chemistry/data/nuclear_hessian.hpp>
+#include <qdk/chemistry/utils/string_utils.hpp>
+
+#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<NuclearHessian> 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<NuclearHessian> 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<NuclearHessian> 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_<NuclearHessian, DataClass, py::smart_holder>(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<std::shared_ptr<Structure>, 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 fd8bfa7c2..5b5d2b1e6 100644
--- a/python/src/pybind11/module.cpp
+++ b/python/src/pybind11/module.cpp
@@ -23,6 +23,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);
@@ -31,6 +33,7 @@ void bind_mc(py::module& m);
 void bind_mcscf(py::module& m);
 void bind_hamiltonian_constructor(py::module& m);
 void bind_scf(py::module& m);
+void bind_nuclear_derivative(py::module& m);
 void bind_active_space(py::module& m);
 void bind_constants(py::module& m);
 void bind_pmc(py::module& m);
@@ -83,6 +86,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
@@ -93,6 +98,7 @@ PYBIND11_MODULE(_core, m) {
   bind_mcscf(algorithms);
   bind_hamiltonian_constructor(algorithms);
   bind_scf(algorithms);
+  bind_nuclear_derivative(algorithms);
   bind_active_space(algorithms);
   bind_dynamical_correlation_calculator(algorithms);
   bind_pmc(algorithms);
diff --git a/python/src/qdk_chemistry/algorithms/__init__.py b/python/src/qdk_chemistry/algorithms/__init__.py
index 59e1444f5..ae3d25dc3 100644
--- a/python/src/qdk_chemistry/algorithms/__init__.py
+++ b/python/src/qdk_chemistry/algorithms/__init__.py
@@ -38,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,
@@ -65,11 +70,13 @@
     "ControlledCircuitMapper",
     "DynamicalCorrelationCalculator",
     "EnergyEstimator",
+    "FiniteDifferenceNuclearDerivativeCalculator",
     "HadamardTest",
     "HamiltonianConstructor",
     "HamiltonianUnitaryBuilder",
     "MultiConfigurationCalculator",
     "MultiConfigurationScf",
+    "NuclearDerivativeCalculator",
     "OrbitalLocalizer",
     "PhaseEstimation",
     "ProjectedMultiConfigurationCalculator",
@@ -81,6 +88,7 @@
     "QdkMacisAsci",
     "QdkMacisCas",
     "QdkMacisPmc",
+    "QdkNuclearDerivativeCalculator",
     "QdkOccupationActiveSpaceSelector",
     "QdkPipekMezeyLocalizer",
     "QdkQubitMapper",
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 7df312a91..ef4c457c1 100644
--- a/python/src/qdk_chemistry/algorithms/registry.py
+++ b/python/src/qdk_chemistry/algorithms/registry.py
@@ -647,6 +647,7 @@ def _register_cpp_factories():
         LocalizerFactory,
         MultiConfigurationCalculatorFactory,
         MultiConfigurationScfFactory,
+        NuclearDerivativeCalculatorFactory,
         ProjectedMultiConfigurationCalculatorFactory,
         ScfSolverFactory,
         StabilityCheckerFactory,
@@ -657,6 +658,7 @@ def _register_cpp_factories():
     register_factory(LocalizerFactory)
     register_factory(MultiConfigurationCalculatorFactory)
     register_factory(MultiConfigurationScfFactory)
+    register_factory(NuclearDerivativeCalculatorFactory)
     register_factory(ProjectedMultiConfigurationCalculatorFactory)
     register_factory(DynamicalCorrelationCalculatorFactory)
     register_factory(ScfSolverFactory)
diff --git a/python/src/qdk_chemistry/data/__init__.py b/python/src/qdk_chemistry/data/__init__.py
index e8acacd9a..0f17857f2 100644
--- a/python/src/qdk_chemistry/data/__init__.py
+++ b/python/src/qdk_chemistry/data/__init__.py
@@ -28,6 +28,8 @@
 - :class:`MeasurementData`: Measurement bitstring data and metadata for QubitHamiltonian objects.
 - :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:`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.
@@ -83,6 +85,8 @@
     LatticeGraph,
     MajoranaMapping,
     ModelOrbitals,
+    NuclearGradients,
+    NuclearHessian,
     Orbitals,
     OrbitalType,
     PauliOperator,
@@ -157,6 +161,8 @@
     "MajoranaMapping",
     "MeasurementData",
     "ModelOrbitals",
+    "NuclearGradients",
+    "NuclearHessian",
     "OrbitalType",
     "Orbitals",
     "PauliOperator",
diff --git a/python/tests/test_nuclear_derivative.py b/python/tests/test_nuclear_derivative.py
new file mode 100644
index 000000000..139afa5cf
--- /dev/null
+++ b/python/tests/test_nuclear_derivative.py
@@ -0,0 +1,68 @@
+"""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"