exchange_commutator.h

Go to the documentation of this file.

//
// exchange_commutator.h — library-style API for evaluating [K̂, f₁₂] |ij⟩
//
// Phase A: thin wrappers around existing CCPotentials::apply_Kfxy /
// apply_KffK_low_rank / apply_KffK_low_rank_direct, plus the scalar
// split-α diagnostic extracted from benchmark_exchange_commutator.cc.
// All entry points return a KffKResult with uniform wall-time / memory /
// rank instrumentation; all diagnostics run through a single diagnose()
// function operating on a harmonic-basis projection of the 6D result.
//
 
#ifndef MADNESS_CHEM_EXCHANGE_COMMUTATOR_H
#define MADNESS_CHEM_EXCHANGE_COMMUTATOR_H
 
#include <madness/chem/CCStructures.h>
#include <madness/chem/ccpairfunction.h>
#include <madness/chem/lowrankfunction.h>
#include <madness/chem/CCParameters.h>
#include <madness/chem/operator_diagnostics.h>
 
namespace madness {
 
struct ExchangeCommutator {
 
    ExchangeCommutator() = default;
 
    /// Uniform result envelope for every algorithm path.  Each "piece" may
    /// be empty when the algorithm only produces the combined commutator
    /// (e.g. lrf-split-alpha constructs KffK directly).
    struct KffKResult {
        std::vector<CCPairFunction<double, 6>> Kf;     ///<  K̂ f₁₂ | i j ⟩
        std::vector<CCPairFunction<double, 6>> fK;     ///<  f₁₂ K̂ | i j ⟩
        std::vector<CCPairFunction<double, 6>> KffK;   ///<  [K̂, f₁₂] | i j ⟩ = Kf − fK (signs baked in)
        double      t_wall = 0.0;                      ///< wall time (s)
        double      mem_gb = 0.0;                      ///< combined size of pair factors (GByte)
        long        rank   = 0;                        ///< combined LRF rank (0 for 6D)
        std::string algo;                              ///< "6d" / "lrf" / "lrf-direct" / "lrf-split-alpha"
    };
 
    /// Options for the split-α LRF k-commutator variant
    /// (see lrf_k_commutator_analysis.md and lrf_k_commutator_results.md).
    struct SplitAlphaOptions {
        double alpha_star = 1.0e4;   ///< GFit partition threshold; α_μ > α* is discarded
        double lo         = 1.0e-6;
        double hi         = 10.0;
        double eps_gfit   = 1.0e-6;
    };
 
    /// Options for the three-range LRF k-commutator variant.
    /// Partition the Coulomb GFit into:
    ///   * diffuse  (α_μ < alpha_lo): wide Gaussians, large-ve LRF, low rank
    ///   * medium   (alpha_lo ≤ α_μ ≤ alpha_hi): tighter LRF
    ///   * tight    (α_μ > alpha_hi): discard (δ-cancellation in commutator)
    /// See lrf_three_range_gfit.md for the derivation.
    struct ThreeRangeOptions {
        double alpha_lo  = 1.0;      ///< diffuse / medium boundary
        double alpha_hi  = 1.0e4;    ///< medium / tight boundary  (≡ alpha_star)
        double lo        = 1.0e-6;
        double hi        = 10.0;
        double eps_gfit  = 1.0e-6;
        /// Per-range LRF construction parameters; both default-constructed.
        /// Caller is expected to set radius / volume_element / tol / lmax /
        /// f12type per range before invocation.
        LowRankFunctionParameters lrfparam_diffuse;
        LowRankFunctionParameters lrfparam_medium;
        /// If true, build the medium-range LRF via project_from_operator
        /// (Taylor / direct construction) instead of random-Y project().
        /// Mandatory for the medium slab in isolation: random-Y stalls at
        /// large error without diffuse anchor functions.
        bool medium_use_taylor = true;
        /// Max Taylor order for project_from_operator (passed straight
        /// through; default matches lowrankfunction.h's own default).
        int  max_taylor_order  = 15;
        /// If true, the medium slab is treated by a *6D* apply with a
        /// partial-Coulomb operator (only the medium Gaussians) instead
        /// of any LRF.  The diffuse slab still uses the LRF path.
        /// Mutually exclusive with `medium_use_taylor`; takes precedence.
        bool medium_use_6d     = false;
    };
 
 
    /// compute the LRF of the exchange operator: K(r,r') = \sum_k k(r)g(r,r') R2(r') k(r')
    ///
    /// apply K to a function f via: Kf(r) = inner(K,f,p2,p1);
    /// note K is not self-adjoint due to the ncf factor
    static LowRankFunction<double,6> compute_lrf_exchange_operator(
            World& world,
            const Info& info,
            const SplitAlphaOptions& opt,
            const LowRankFunctionParameters& lrfparam);
 
    // ---------------------------------------------------------------------
    // Top-level entry points — each produces a KffKResult wrapping one of
    // the four existing algorithm paths.
    // ---------------------------------------------------------------------
 
    /// 6D algorithm: Kf via apply_Kfxy, fK via make_f_xy_macrotask(Kx,...).
    /// Mirrors the body of test_kcomm_6d in the original benchmark file.
    static KffKResult apply_KffK_6d(
            World& world,
            const CCFunction<double, 3>& phi_i,
            const CCFunction<double, 3>& phi_j,
            const Info& info,
            const real_convolution_6d* Gscreen = nullptr);
 
    // ---------------------------------------------------------------------
    // Split-α k-commutator.  Scalar-only variant for Phase A — 6D pair
    // assembly is deferred.
    // ---------------------------------------------------------------------
 
    /// Assemble the 6D exchange commutator [K̂, f] |ij⟩ using the split-α
    /// LRF k-commutator algorithm.  Sums over occupied orbitals in
    /// info.mo_bra/info.mo_ket.  Currently assumes phi_i == phi_j
    /// (symmetric case) and exploits swap_particles for the K̂₂ piece.
    ///
    /// Output pair layout:
    ///   result.pair[0] = Σ_k Σ_ρ  g_kρ(1) · (phi_j · A_kρ)(2)       (separable)
    ///   result.pair[1] = Σ_k Σ_ρ  B_kρ · (g_kρ·phi_i)(1) · f(1,2) · phi_j(2)
    ///                    with coefficient −1
    ///   plus their particle-2 mirrors when i == j.
    static KffKResult apply_KffK_lowrank_split_alpha(
            World& world,
            const CCFunction<double, 3>& phi_i,
            const CCFunction<double, 3>& phi_j,
            const Info& info,
            const LowRankFunction<double,6>& exchange_op,
            const LowRankFunctionParameters& lrfparam,
            const SplitAlphaOptions& opt);
 
    /// Three-range LRF k-commutator: builds two LRFs (diffuse + medium) with
    /// independent LowRankFunctionParameters, concatenates their g/h vectors
    /// and assembles Kf / fK / KffK with the same algebra as
    /// apply_KffK_lowrank_split_alpha.  Tight terms (α_μ > opt.alpha_hi) are
    /// discarded.  Returns KffKResult with algo = "lrf-three-range" and rank
    /// = rank_diffuse + rank_medium (printed individually in t_wall summary).
    static KffKResult apply_KffK_lowrank_three_range(
            World& world,
            const CCFunction<double, 3>& phi_i,
            const CCFunction<double, 3>& phi_j,
            const Info& info,
            const ThreeRangeOptions& opt);
 
    // ---------------------------------------------------------------------
    // [K̂,f] diagnostic — plain ⟨ab| projection vs analytic 3D reference
    // ---------------------------------------------------------------------
 
    /// Diagnose ⟨ab | [K̂,f₁₂] | ij⟩ (no G, no Q₁₂): the 6D pieces are
    /// projected onto the observer basis and compared against the analytic
    /// Kf/fK matrix elements from kf_provider/fk_provider evaluated over the
    /// plain ⟨ab| bra.  Unassigned CCPairFunctions are skipped.
    ///
    /// NOTE: diagnose_KffK, diagnose_GKffK and diagnose_GQKffK share the same
    /// signature; parameters not needed by a particular diagnostic are
    /// ignored (here: energy — pass any value).
    /// Returns a DiagnosticMatrix with entries "Kf", "fK", "KffK" (= Kf − fK).
    /// ref=3D Exchange+f12 analytic formula, result=6D projection, error=||ref-result||.
    DiagnosticMatrix<> diagnose_KffK(
            World& world,
            const std::vector<CCPairFunction<double,6>>& Kf_cc,
            const std::vector<CCPairFunction<double,6>>& fK_cc,
            const real_function_3d& phi_i,
            const real_function_3d& phi_j,
            const real_function_3d& Kphi_i,
            const real_function_3d& Kphi_j,
            const Info& info,
            double energy,
            const std::vector<real_function_3d>& aobasis,
            bool orthonormalize_basis = true) const;
 
    // ---------------------------------------------------------------------
    // G·[K̂,f] diagnostics — Schwinger-quadrature 3D reference vs 6D projection
    // ---------------------------------------------------------------------
 
    /// Diagnose ⟨ab | G [K̂,f₁₂] | ij⟩ by comparing a 6D projection against
    /// a 3D Schwinger eigentime quadrature.
    /// Returns a DiagnosticMatrix with entries "GKf", "GfK", "GKffK".
    /// ref=3D Schwinger quadrature, result=6D projection <ab|G·piece|ij>, error=||ref-result||.
    ///
    /// The 6D references GKf_cc and GfK_cc are the G-applied commutator pieces
    /// as CCPairFunction vectors.  Accepting CCPairFunction (instead of a bare
    /// real_function_6d) allows any algorithm variant to be tested:
    ///
    ///   * 6D path (apply_KffK_6d): Kf and fK are pure-6D; caller applies G6d
    ///     via apply(G6d, kffk.Kf[0].get_function()) and wraps in a vector.
    ///   * LRF path (apply_KffK_lowrank_split_alpha, three_range): Kf may be
    ///     decomposed.  Caller converts to pure (sum outer products), applies G6d,
    ///     and passes the result wrapped in a CCPairFunction vector.
    ///   * Mixed/multi-piece results: sum all pieces (each entry is summed inside).
    ///
    /// The 6D reference is projected using the same partial_inner / matrix_inner
    /// machinery as diagnose(), handling both pure and decomposed CCPairFunctions.
    ///
    /// The 3D Schwinger formulas use only K̂, K̂†, and the Slater f₁₂:
    ///
    ///   G Kf:  ⟨ã b̃ | K̂₁ f₁₂ | φᵢ φⱼ⟩ = inner((K̂†ã)·φᵢ, f₁₂(b̃·φⱼ))  + K̂₂ swap
    ///   G fK:  ⟨ã b̃ | f₁₂ K̂₁ | φᵢ φⱼ⟩ = inner(ã·K̂φᵢ,    f₁₂(b̃·φⱼ))  + K̂₂ swap
    ///
    /// where ã = g_{αₙ}*a and b̃ = g_{αₙ}*b are Gaussian-smoothed AO basis
    /// functions and the sum over quadrature nodes n (with 6D weight
    /// w_n = c_n^{bsh}·(αₙ/π)^{3/2}) approximates G.
    ///
    /// NOTE: diagnose_GKffK and diagnose_GQKffK share the same signature.
    /// @param GKf_cc   G·K̂f₁₂|ij⟩ pieces (caller applies G6d externally)
    /// @param GfK_cc   G·f₁₂K̂|ij⟩ pieces
    /// @param Kphi_i   K̂φᵢ (precomputed once by caller)
    /// @param Kphi_j   K̂φⱼ (precomputed once by caller)
    /// @param info     molecular/orbital info carrying mo_ket, mo_bra, parameters
    /// @param energy   ε_i + ε_j  (must be negative)
    /// @param aobasis  observer basis
    /// @param orthonormalize_basis  Löwdin-orthonormalize aobasis (false: keep raw)
    DiagnosticMatrix<> diagnose_GKffK(
            World& world,
            const std::vector<CCPairFunction<double,6>>& GKf_cc,
            const std::vector<CCPairFunction<double,6>>& GfK_cc,
            const real_function_3d& phi_i,
            const real_function_3d& phi_j,
            const real_function_3d& Kphi_i,
            const real_function_3d& Kphi_j,
            const Info& info,
            double energy,
            const std::vector<real_function_3d>& aobasis,
            bool orthonormalize_basis = true) const;
 
    // ---------------------------------------------------------------------
    // G·Q₁₂·[K̂,f] diagnostic — 3D-only reference, Q₁₂ expanded on the ket
    // ---------------------------------------------------------------------
 
    /// Diagnose ⟨ab | G Q₁₂ [K̂,f₁₂] | ij⟩, Q₁₂ = (1-O₁)(1-O₂), O = Σ|k_ket⟩⟨k_bra|.
    ///
    /// Caller supplies G·Q₁₂·K̂f₁₂|ij⟩ and G·Q₁₂·f₁₂K̂|ij⟩ (neither G nor Q is
    /// applied internally).  The reference is computed via 3D functions only:
    /// G is moved to the bra via the Schwinger/BSH fit and Q₁₂ is expanded on
    /// the ket side as scalar contractions of analytic Kf/fK matrix elements —
    /// see operator_diagnostics.md for the full derivation.  Entries: "GQKf",
    /// "GQfK", "GQKffK" (= GQKf − GQfK).
    ///
    /// The observer basis is an explicit parameter (unlike diagnose_GKffK):
    /// pass info.ao for the accuracy check, or the raw pair-bra functions
    /// {mo_bra[i], mo_bra[j]} with orthonormalize_basis=false so that the
    /// (0,1) element ⟨ī j̄|G Q₁₂ [K̂,f]|ij⟩ is the exchange-commutator
    /// contribution to the MP2 pair energy.
    ///
    /// @param GQKf_cc  G Q₁₂ K̂f₁₂|ij⟩ pieces (caller applies Q12 and G externally)
    /// @param GQfK_cc  G Q₁₂ f₁₂K̂|ij⟩ pieces
    /// @param Kphi_i   K̂φᵢ (precomputed once by caller)
    /// @param Kphi_j   K̂φⱼ (precomputed once by caller)
    /// @param info     carries mo_ket, mo_bra (Q₁₂ and K̂ spaces) and parameters
    /// @param energy   ε_i + ε_j  (must be negative)
    /// @param aobasis  observer basis
    /// @param orthonormalize_basis  Löwdin-orthonormalize aobasis (false: keep raw)
    DiagnosticMatrix<> diagnose_GQKffK(
            World& world,
            const std::vector<CCPairFunction<double,6>>& GQKf_cc,
            const std::vector<CCPairFunction<double,6>>& GQfK_cc,
            const real_function_3d& phi_i,
            const real_function_3d& phi_j,
            const real_function_3d& Kphi_i,
            const real_function_3d& Kphi_j,
            const Info& info,
            double energy,
            const std::vector<real_function_3d>& aobasis,
            bool orthonormalize_basis = true) const;
 
    /// element provider: M(x,y) = Σ_p ⟨p1[x](1) p2[y](2)|K̂_p f₁₂|φᵢφⱼ⟩
    ///   K̂₁: ⟨(K̂†x)·φᵢ | f₁₂⋆(y·φⱼ)⟩      K̂₂: ⟨f₁₂⋆(x·φᵢ) | (K̂†y)·φⱼ⟩
    /// K̂† (Exchange with bra=mo_ket, ket=mo_bra) and the f₁₂ convolution are
    /// built once from `info` and captured in the closure.
    static DiagnosticMatrix<>::ElementProvider kf_provider(
            World& world,
            const real_function_3d& phi_i,
            const real_function_3d& phi_j,
            const Info& info,
            double eps = -1.0);
 
    /// element provider: M(x,y) = Σ_p ⟨p1[x](1) p2[y](2)|f₁₂ K̂_p|φᵢφⱼ⟩
    ///   K̂₁: ⟨x·(K̂φᵢ) | f₁₂⋆(y·φⱼ)⟩       K̂₂: ⟨x·φᵢ | f₁₂⋆(y·(K̂φⱼ))⟩
    /// The f₁₂ convolution is built once from `info` and captured in the closure.
    static DiagnosticMatrix<>::ElementProvider fk_provider(
            World& world,
            const real_function_3d& phi_i,
            const real_function_3d& phi_j,
            const real_function_3d& Kphi_i,
            const real_function_3d& Kphi_j,
            const Info& info,
            double eps = -1.0);
};
 
} // namespace madness
 
#endif // MADNESS_CHEM_EXCHANGE_COMMUTATOR_H