Add SW-PT2 downfolding

cpp/include/qdk/chemistry/algorithms/effective_hamiltonian.hpp

@@ -22,7 +22,9 @@ namespace qdk::chemistry::algorithms {
  * specified target P-space.
  *
  * @c p_indices holds absolute molecular-orbital indices, drawn from the same
- * index universe as @c data::Orbitals::active_indices().
+ * index universe as @c data::Orbitals::active_indices(). They select the target
+ * space out of the input Hamiltonian's active window, independently of the
+ * reference wavefunction's own active space.
  *
  * The returned Hamiltonian is expressed over P and must satisfy:
  * - its orbitals have @c active_indices() equal to @c p_indices;
@@ -160,7 +162,7 @@ struct EffectiveHamiltonianConstructorFactory
     return "effective_hamiltonian_constructor";
   }
   static void register_default_instances();
-  static std::string default_algorithm_name() { return ""; }
+  static std::string default_algorithm_name() { return "qdk_swpt2"; }
 };
 
 }  // namespace qdk::chemistry::algorithms

cpp/src/qdk/chemistry/algorithms/CMakeLists.txt

@@ -11,8 +11,8 @@ target_sources(chemistry PRIVATE
     population_analysis.cpp
     qdk_nuclear_derivative.cpp
     dynamical_correlation_calculator.cpp
-    geometry_optimization.cpp
     effective_hamiltonian.cpp
+    geometry_optimization.cpp
     scf.cpp
     stability.cpp
 )

cpp/src/qdk/chemistry/algorithms/effective_hamiltonian.cpp

@@ -8,9 +8,12 @@
 #include <qdk/chemistry/data/basis_set.hpp>
 #include <qdk/chemistry/data/orbitals.hpp>
 #include <qdk/chemistry/data/symmetry/spin_channel_indices.hpp>
+#include <qdk/chemistry/utils/logger.hpp>
 #include <stdexcept>
 #include <string>
 
+#include "microsoft/effective_hamiltonian/swpt2.hpp"
+
 namespace qdk::chemistry::algorithms {
 
 namespace {
@@ -146,6 +149,15 @@ void EffectiveHamiltonianConstructor::_validate_inputs(
   }
 }
 
-void EffectiveHamiltonianConstructorFactory::register_default_instances() {}
+std::unique_ptr<EffectiveHamiltonianConstructor> make_swpt2_constructor() {
+  QDK_LOG_TRACE_ENTERING();
+  return std::make_unique<microsoft::SchriefferWolffPT2Constructor>();
+}
+
+void EffectiveHamiltonianConstructorFactory::register_default_instances() {
+  QDK_LOG_TRACE_ENTERING();
+  EffectiveHamiltonianConstructorFactory::register_instance(
+      &make_swpt2_constructor);
+}
 
 }  // namespace qdk::chemistry::algorithms

cpp/src/qdk/chemistry/algorithms/microsoft/CMakeLists.txt

@@ -20,6 +20,8 @@ target_sources(chemistry PRIVATE
     scf.cpp
     stabilized_scf.cpp
     stability.cpp
+    effective_hamiltonian/swpt2_kernel.cpp
+    effective_hamiltonian/swpt2.cpp
     utils.cpp
 )
 

cpp/src/qdk/chemistry/algorithms/microsoft/effective_hamiltonian/swpt2.hpp

@@ -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.
+
+#pragma once
+#include <cmath>
+#include <limits>
+#include <memory>
+#include <qdk/chemistry/algorithms/effective_hamiltonian.hpp>
+#include <qdk/chemistry/data/settings.hpp>
+#include <string>
+#include <vector>
+
+namespace qdk::chemistry::algorithms::microsoft {
+
+/**
+ * @brief Settings for the second-order Schrieffer-Wolff downfold.
+ *
+ * `regularizer_sigma2` damps the energy denominators; zero selects a guarded
+ * bare pseudoinverse that omits coupled channels below the fixed 1e-8 Eh
+ * denominator floor. The current denominator operator is a semicanonical,
+ * spin-free generalized Fock.
+ */
+class SchriefferWolffPT2Settings : public qdk::chemistry::data::Settings {
+ public:
+  SchriefferWolffPT2Settings() {
+    set_default(
+        "regularizer_sigma2", 1.0,
+        "Sigma parameter of the sigma^2 denominator regularizer, "
+        "1/D -> (1-exp(-sigma*D^2))/D, in units of Eh^-2; equivalently the "
+        "DSRG flow parameter. Larger values regularize less; set to 0 to "
+        "disable. The unregularized calculation uses a guarded pseudoinverse "
+        "that maps coupled denominators smaller than 1e-8 Eh to zero and "
+        "logs a warning. This borrows the DSRG damping form; it does not turn "
+        "the downfold into a full DSRG calculation.",
+        data::BoundConstraint<double>{0.0, std::numeric_limits<double>::max()});
+    set_default("semicanonicalize", true,
+                "Diagonalize the generalized Fock independently within the "
+                "inactive, active, and virtual blocks before forming Fock "
+                "denominators.");
+    set_default(
+        "fold_above_two_body", true,
+        "Approximate the three-body terms generated by the transformation "
+        "with pair contractions from the reference 1-RDM instead of "
+        "discarding them. This is a Gaussian-reference approximation that "
+        "neglects higher density cumulants, and it makes the downfold "
+        "more expensive, increasingly so as the kept space grows.");
+    set_default(
+        "max_folded_occupation_deviation", 0.5,
+        "Largest allowed deviation from an integer reference occupation (0 or "
+        "2) for an orbital folded into the external space. Folded occupations "
+        "are rounded to the nearest of 0 or 2; the total electron count is "
+        "preserved because the active space receives whatever the folded "
+        "orbitals do not take. Must be below 1, so a singly occupied orbital "
+        "is never folded on an arbitrary rounding.",
+        data::BoundConstraint<double>{0.0, std::nextafter(1.0, 0.0)});
+  }
+  ~SchriefferWolffPT2Settings() override = default;
+};
+
+/**
+ * @brief Second-order Schrieffer-Wolff (Van Vleck) effective-Hamiltonian
+ * downfold with semicanonical generalized-Fock orbital-energy denominators.
+ *
+ * Computes `H_eff = H_BD + 1/2 [S, H_OD]`, emitting <= 2-body, folding the
+ * window's external space Q onto its kept space P. S solves
+ * `[F0, S] = H_OD` for a diagonal generalized-Fock F0; a regularizer setting
+ * replaces the bare inverse denominators to damp intruder-state channels.
+ *
+ * The commutator also generates three-body terms, which a Hamiltonian cannot
+ * hold. Rather than discarding them, `fold_above_two_body` contracts them with
+ * the reference 1-RDM and keeps what falls to two-body. For a Gaussian
+ * reference these are the ordinary Wick contractions, but the residual
+ * three-body operator is still discarded. For a correlated reference the
+ * approximation additionally neglects the two-body and higher density
+ * cumulants.
+ *
+ * The implementation assumes a common restricted MO basis, however it was
+ * produced: RHF, ROHF, CAS, localized and natural orbitals all qualify. Every
+ * singly occupied orbital must lie in the reference active space. Noncanonical
+ * orbitals are semicanonicalized independently within inactive, active, and
+ * virtual blocks. Intruder diagnostics are logged, with an additional warning
+ * for large raw amplitudes.
+ *
+ * The emitted two-body block is only 4-fold symmetric: hermiticity and electron
+ * exchange survive, the bra swap `(pq|rs) = (qp|rs)` does not. Consumers must
+ * read the full dense norb^4 block, not the canonical 8-fold-unique elements.
+ * See `to_spatial_chemist` in `swpt2_kernel.hpp`.
+ *
+ * The window W is the input Hamiltonian's own active space, so that Hamiltonian
+ * must be built over the whole window; a P-only active-space Hamiltonian has
+ * already folded away the couplings this transforms. The reference supplies
+ * only the density over W (its active 1-RDM, with inactive orbitals doubly
+ * occupied and the rest empty), never the partition: P comes from `p_indices`
+ * alone and need not coincide with the reference active space. Occupations of
+ * folded orbitals are rounded, bounded by `max_folded_occupation_deviation`.
+ *
+ * See `swpt2_kernel.hpp` for the operator and tensor conventions.
+ */
+class SchriefferWolffPT2Constructor
+    : public qdk::chemistry::algorithms::EffectiveHamiltonianConstructor {
+ public:
+  SchriefferWolffPT2Constructor() {
+    _settings = std::make_unique<SchriefferWolffPT2Settings>();
+  }
+  ~SchriefferWolffPT2Constructor() override = default;
+
+  std::string name() const final { return "qdk_swpt2"; }
+
+ protected:
+  std::shared_ptr<data::Hamiltonian> _run_impl(
+      std::shared_ptr<data::Wavefunction> reference,
+      std::shared_ptr<data::Hamiltonian> hamiltonian,
+      std::shared_ptr<const data::SymmetryBlockedIndexSet> p_indices)
+      const override;
+};
+
+}  // namespace qdk::chemistry::algorithms::microsoft