// Γ-only restricted Hartree-Fock for a periodic system, molecular-limit
// regime.
//
// Scope note. This driver solves a single-k (Γ) closed-shell SCF with the
// real-space density taken as Γ-only — i.e., P_μν(g) is nonzero only for
// g = 0. That's physically meaningful in the molecular-limit regime
// (unit cell large enough that cross-cell density matrix elements are
// numerically zero), where it must reproduce molecular RHF to machine
// precision.
//
// Bulk calculations with genuinely metallic or covalent inter-cell
// bonding need multi-k sampling of P(g) — that arrives in Phase 12c.
//
// Implementation reuses the existing DIIS machinery: at Γ the Fock /
// overlap / density matrices are real nbf × nbf, so Pulay's FDS - SDF
// error vector and the linear-combination extrapolation work unchanged.

#pragma once

#include <Eigen/Dense>
#include <cstddef>
#include <vector>

#include "basis.hpp"
#include "dft_plus_u.hpp"  // HubbardSiteCxx
#include "ediis.hpp"    // SCFAccelerator
#include "guess.hpp"
#include "lattice_sum.hpp"
#include "periodic.hpp"
#include "rhf.hpp"      // SCFIteration

namespace vibeqc {

struct PeriodicRHFOptions {
    // SCF convergence controls (same shape as molecular RHFOptions).
    int max_iter = 100;
    double conv_tol_energy = 1.0e-8;
    double conv_tol_grad = 1.0e-6;
    double damping = 0.5;
    // Dynamic damping (Zerner-Hehenberger 1979). See RHFOptions::dynamic_damping.
    bool dynamic_damping = false;
    double dynamic_damping_min = 0.0;
    double dynamic_damping_max = 0.95;
    // Fock matrix mixing; see RHFOptions::fock_mixing. This mirrors
    // CRYSTAL's FMIXING keyword and is independent of density damping.
    double fock_mixing = 0.0;

    bool use_diis = true;
    int diis_start_iter = 2;
    std::size_t diis_subspace_size = 8;

    // SCF Fock-extrapolation accelerator (DIIS / KDIIS / EDIIS /
    // EDIIS_DIIS). See RHFOptions::scf_accelerator and
    // cpp/include/vibeqc/ediis.hpp for the family. Default DIIS
    // (back-compat). Ignored when use_diis = false.
    SCFAccelerator scf_accelerator = SCFAccelerator::DIIS;
    // EDIIS → DIIS switch threshold for SCFAccelerator::EDIIS_DIIS.
    // See RHFOptions::ediis_diis_switch_threshold.
    double ediis_diis_switch_threshold = 1e-1;

    // v0.9.x: default flipped to AUTO. The unified GuessEngine
    // (cpp/src/guess.cpp) resolves AUTO → SAD for any periodic system
    // — fixes the historic NaCl/MgO bombing class of failure where
    // Hcore's neglect of electron repulsion sends DIIS into
    // +30 000 / −16 000 Ha territory before recovery. The v0.6.2
    // calibration freeze that kept Hcore in place is lifted; tests
    // that were calibrated to Hcore iteration counts have been refreshed.
    InitialGuess initial_guess = InitialGuess::AUTO;

    // Phase C1a — Saunders-Hillier level shift (Hartree). Adds
    // ``b · (S - ½ S D S)`` to F before diagonalisation, which
    // raises virtual orbital eigenvalues by ``b`` while leaving
    // occupied orbitals unchanged. The shift is "inert at the
    // converged density" — it doesn't change the SCF fixed point,
    // only the iteration dynamics. Useful when DIIS oscillates
    // between near-degenerate occupied / virtual swaps on small-
    // HOMO–LUMO-gap insulators or weakly metallic cells. Default
    // 0.0 (no shift). Typical values: 0.1 – 0.5 Hartree.
    double level_shift = 0.0;

    // CRYSTAL-style level-shift warm-up. ``-1`` means auto: if
    // level_shift is nonzero, use up to five shifted startup cycles
    // and leave at least one unshifted tail cycle. ``0`` keeps the
    // historic persistent level-shift behavior. Positive values request
    // an explicit warm-up length.
    int level_shift_warmup_cycles = -1;

    // Phase C1b — Fermi-Dirac smearing temperature (Hartree).
    // ``T > 0`` replaces hard Aufbau (occupations ∈ {0, 2}) with a
    // smooth Fermi function ``n_i = 2 / (1 + exp((ε_i - μ) / T))``
    // and bisects for ``μ`` to satisfy the per-cell electron-count
    // constraint. Required for metals and small-gap insulators where
    // Aufbau gives oscillating fixed points. Adds an electronic-
    // entropy contribution to the free energy ``A = E - T S``;
    // converged calculations report both the internal energy and the
    // free energy. Default 0.0 (no smearing → standard Aufbau).
    // Typical values: 0.001–0.02 Hartree (~315–6300 K).
    double smearing_temperature = 0.0;

    // Phase C1c — second-order ("quadratic") SCF fallback.
    //
    // When standard SCF (damping + DIIS + level shift) fails to
    // converge — typically on small-gap insulators where DIIS
    // oscillates between near-degenerate occ/vir swaps — switch from
    // "diagonalize F" to a Newton step in MO space:
    //
    //   κ_{ai}  =  -F_{ai}^{MO} / (ε_a - ε_i + quadratic_fallback_shift)
    //   C_new   =  C_prev · exp(κ)            (skew-symmetric κ)
    //   D_new   =  2 · C_new_occ · C_new_occ^T
    //
    // The step is preconditioned by orbital-energy differences (the
    // diagonal of the orbital-rotation Hessian) and trust-region
    // capped at ``quadratic_fallback_max_step`` to keep the Taylor
    // exponential well-conditioned.
    //
    // Activation. ``quadratic_fallback_iter > 0`` enables the
    // fallback after iteration ``quadratic_fallback_iter``.
    // ``= 0`` (default) disables — SCF behaves exactly as before.
    // DIIS and level shift are skipped during the quadratic phase
    // (the Newton step is its own update mechanism, and mixing
    // with DIIS extrapolation undoes the trust-region cap).
    //
    // When this helps. Insulators with HOMO-LUMO gap < ~0.1 eV;
    // open-shell systems with near-degenerate ground / excited
    // states; high-symmetry crystals where many MO pairs are close.
    int quadratic_fallback_iter = 0;
    double quadratic_fallback_shift = 0.1;
    double quadratic_fallback_max_step = 0.1;

    // Lattice-sum controls applied to one-electron and J/K builds alike.
    LatticeSumOptions lattice_opts = {};

    // DFT+U (Dudarev). Γ-only periodic +U uses the *same* real-density
    // kernel as the molecular path: at Γ the density is a single real
    // nbf x nbf matrix, the AO overlap is a single real matrix, and
    // ao_group_indices(basis) returns the same shape because the basis
    // is built from PeriodicSystem.unit_cell_molecule(). Empty by
    // default. Multi-k periodic +U lands in Increment 4c — see
    // HANDOVER_DFT_PLUS_U.md.
    std::vector<HubbardSiteCxx> dft_plus_u_sites;
    std::vector<std::vector<int>> dft_plus_u_ao_groups;
};

struct PeriodicRHFResult {
    double energy = 0.0;           // Total Hartree–Fock energy per unit cell
    double e_electronic = 0.0;
    double e_nuclear = 0.0;        // Sum over unit-cell pairs within nuclear
                                   // cutoff (matches V lattice sum)
    // Dudarev DFT+U contribution per unit cell. 0 when +U is not active.
    double e_dft_plus_u = 0.0;
    int n_iter = 0;
    bool converged = false;

    // Γ-point molecular orbitals.
    Eigen::VectorXd mo_energies;   // ascending (Hartree)
    Eigen::MatrixXd mo_coeffs;     // real at Γ

    // Final Γ-aggregated matrices.
    Eigen::MatrixXd density;
    Eigen::MatrixXd fock;
    Eigen::MatrixXd overlap;

    std::vector<SCFIteration> scf_trace;
};

// Nuclear-repulsion energy per unit cell, with the reference cell's nuclei
// interacting with all images within the lattice-sum cutoff. A half-sum
// handles the pair double-counting correctly: each pair of unit-cell atoms
// is counted once across all cell shifts (including g = 0, which reduces
// to the intra-cell molecular repulsion).
double nuclear_repulsion_per_cell(const PeriodicSystem& system,
                                  const LatticeSumOptions& opts);

PeriodicRHFResult run_rhf_periodic_gamma(const PeriodicSystem& system,
                                         const BasisSet& basis,
                                         const PeriodicRHFOptions& options = {});

}  // namespace vibeqc