module tdhf_mrsf_ekt_mod

  implicit none

  character(len=*), parameter :: module_name = "tdhf_mrsf_ekt_mod"

contains

  subroutine tdhf_mrsf_ekt_ip_C(c_handle) bind(C, name="tdhf_mrsf_ekt_ip")
    use c_interop, only: oqp_handle_t, oqp_handle_get_info
    use types, only: information
    type(oqp_handle_t) :: c_handle
    type(information), pointer :: inf
    inf => oqp_handle_get_info(c_handle)
    call tdhf_mrsf_ekt(inf, .false.)
  end subroutine tdhf_mrsf_ekt_ip_C

  subroutine tdhf_mrsf_ekt_ea_C(c_handle) bind(C, name="tdhf_mrsf_ekt_ea")
    use c_interop, only: oqp_handle_t, oqp_handle_get_info
    use types, only: information
    type(oqp_handle_t) :: c_handle
    type(information), pointer :: inf
    inf => oqp_handle_get_info(c_handle)
    call tdhf_mrsf_ekt(inf, .true.)
  end subroutine tdhf_mrsf_ekt_ea_C

  subroutine tdhf_mrsf_ekt(infos, electron_affinity)
    use io_constants, only: iw
    use oqp_tagarray_driver
    use types, only: information, energy_results
    use basis_tools, only: basis_set
    use messages, only: show_message, with_abort
    use util, only: measure_time
    use precision, only: dp
    use mathlib, only: unpack_matrix, pack_matrix, orthogonal_transform_sym
    use eigen, only: diag_symm_full
    use printing, only: print_module_info
    use tdhf_mrsf_energy_mod, only: tdhf_mrsf_energy
    use tdhf_mrsf_z_vector_mod, only: tdhf_mrsf_z_vector
    use dft, only: dft_initialize
    use mod_dft_molgrid, only: dft_grid_t
    use scf_addons, only: calc_fock, scf_energy_t

    implicit none

    character(len=*), parameter :: subroutine_name = "tdhf_mrsf_ekt"
    type(information), target, intent(inout) :: infos
    logical, intent(in) :: electron_affinity

    type(basis_set), pointer :: basis
    integer :: nbf, nbf2, nroot, ok, i, nfocks, nschwz
    integer :: nactive
    integer, allocatable :: dom_mo_idx(:), dom_no_idx(:)
    ! EKT natural-orbital deflation threshold, matching GAMESS EKT-MRSF
    ! (thresh = 1.0d-2, sfgrad.src:3298).
    real(kind=dp), parameter :: ekt_occ_tol = 1.0e-2_dp
    real(kind=dp), parameter :: hartree_to_ev = 27.211386245988_dp
    real(kind=dp), allocatable :: fock_pack(:)
    type(dft_grid_t), target :: ea_molgrid
    type(scf_energy_t) :: ea_energy
    type(energy_results) :: saved_mol_energy
    real(kind=dp), allocatable :: density_alpha_ao(:), density_beta_ao(:)
    real(kind=dp), allocatable :: density_alpha_mo(:,:), density_beta_mo(:,:)
    real(kind=dp), allocatable :: density_ip_mo(:,:), density_ea_mo(:,:)
    real(kind=dp), allocatable :: smat(:)
    real(kind=dp), allocatable :: density_mo(:,:), fock_mo(:,:), lagrangian_mo(:,:)
    real(kind=dp), allocatable :: ea_density_ao(:,:), ea_fock_ao(:,:), saved_fock_ao(:,:)
    real(kind=dp), allocatable :: ekt_metric(:,:), ekt_operator(:,:)
    real(kind=dp), allocatable :: eig(:)
    real(kind=dp), allocatable :: orbitals(:,:), strengths(:), metric_norms(:)
    real(kind=dp), allocatable :: dom_mo_coeff(:), dom_no_coeff(:), dom_no_occ(:)
    real(kind=dp), contiguous, pointer :: fock_a(:), fock_b(:), dmat_a(:), dmat_b(:), mo_a(:,:), td_p(:,:), wao(:)
    real(kind=dp), contiguous, pointer :: density_store(:,:), lagrangian_store(:,:)
    real(kind=dp), contiguous, pointer :: fock_store(:,:), orbital_store(:,:), eig_store(:), strength_store(:)
    character(len=*), parameter :: tags_required(7) = (/ character(len=80) :: &
      OQP_FOCK_A, OQP_DM_A, OQP_DM_B, OQP_VEC_MO_A, OQP_td_p, OQP_WAO, OQP_td_bvec_mo /)
    character(len=*), parameter :: tags_alloc(6) = (/ character(len=80) :: &
      OQP_mrsf_ekt_density_mo, OQP_mrsf_ekt_lagrangian_mo, OQP_mrsf_ekt_fock_mo, &
      OQP_mrsf_ekt_orbitals_mo, OQP_mrsf_ekt_eigenvalues, OQP_mrsf_ekt_strengths /)

    ! Run the parent MRSF calculation and Z-vector first so the selected
    ! neutral state's relaxed density P and energy-weighted density /
    ! Lagrangian W are available.  The EKT step then builds the IP/EA
    ! generalized eigenproblem in the MRSF molecular-orbital basis.
    call tdhf_mrsf_energy(infos)
    call tdhf_mrsf_z_vector(infos)

    open(unit=iw, file=infos%log_filename, position="append")
    if (electron_affinity) then
      call print_module_info('MRSF_EKT_EA','Computing MRSF-EKT electron affinities')
    else
      call print_module_info('MRSF_EKT_IP','Computing MRSF-EKT ionization potentials')
    end if

    basis => infos%basis
    nbf = basis%nbf
    nbf2 = nbf*(nbf+1)/2
    ! Print/store all retained EKT Dyson roots available after the natural-
    ! orbital deflation, matching the GAMESS behavior.  The TDDFT NSTATE input
    ! still controls the parent MRSF calculation; it is not an EKT print cap.
    nroot = nbf
    select case (infos%control%scftype)
    case (1)
      nfocks = 1
    case default
      nfocks = 2
    end select

    call data_has_tags(infos%dat, tags_required, module_name, subroutine_name, WITH_ABORT)
    call tagarray_get_data(infos%dat, OQP_FOCK_A, fock_a)
    call tagarray_get_data(infos%dat, OQP_DM_A, dmat_a)
    call tagarray_get_data(infos%dat, OQP_DM_B, dmat_b)
    call tagarray_get_data(infos%dat, OQP_VEC_MO_A, mo_a)
    call tagarray_get_data(infos%dat, OQP_td_p, td_p)
    call tagarray_get_data(infos%dat, OQP_WAO, wao)

    allocate(fock_pack(nbf2), &
             ea_density_ao(nbf2,nfocks), ea_fock_ao(nbf2,nfocks), saved_fock_ao(nbf2,nfocks), &
             density_alpha_ao(nbf2), density_beta_ao(nbf2), &
             density_alpha_mo(nbf,nbf), density_beta_mo(nbf,nbf), &
             density_ip_mo(nbf,nbf), density_ea_mo(nbf,nbf), &
             density_mo(nbf,nbf), fock_mo(nbf,nbf), &
             lagrangian_mo(nbf,nbf), ekt_metric(nbf,nbf), ekt_operator(nbf,nbf), &
             orbitals(nbf,nroot), strengths(nroot), smat(nbf2), source=0.0_dp, stat=ok)
    if (ok /= 0) call show_message('Cannot allocate memory', WITH_ABORT)

    ! MRSF-EKT uses the relaxed target-state one-particle density P and the
    ! relaxed energy-weighted density / Lagrangian W generated by the MRSF
    ! Z-vector.  OQP_td_p is the MRSF Z-vector relaxed *difference* density in
    ! the AO basis, spin-blocked (alpha,beta), in the same contravariant
    ! C*M*C^T convention as dmat and already 0.5-scaled in tdhf_mrsf_z_vector.
    ! Hence dmat_sigma + td_p(:,sigma) is the total relaxed one-particle
    ! density of the target (S0) response state for each spin.
    density_alpha_ao = dmat_a + td_p(:,1)
    density_beta_ao = dmat_b + td_p(:,2)

    ! AO -> MO transform of the one-particle DENSITY.  A density is
    ! contravariant: its natural-orbital (MO) representation is obtained with
    ! the metric-covariant form  D_MO = C^T S P S C, NOT C^T P C.  The SCF
    ! MOs satisfy C^T S C = I, so plain C^T P C mis-scales by the overlap and
    ! gives Tr(D_MO) != N_elec.  Operators (Fock, Lagrangian) keep the plain
    ! C^T A C transform below.
    call get_overlap_matrix(infos, nbf, smat)
    call density_ao_to_mo(nbf, nbf2, density_alpha_ao, smat, mo_a, density_alpha_mo)
    call density_ao_to_mo(nbf, nbf2, density_beta_ao,  smat, mo_a, density_beta_mo)
    density_ip_mo = density_alpha_mo
    density_ea_mo = density_beta_mo

    ! Acceptance checks (printed): AO electron counts Tr(P_sigma * S) and the
    ! MO-density traces must both equal nelec_sigma, and MO occupations must
    ! be physically bounded.
    call ekt_density_trace_checks(infos, nbf, nbf2, density_alpha_ao, &
         density_beta_ao, smat, density_alpha_mo, density_beta_mo, iw)

    call orthogonal_transform_sym(nbf, nbf, fock_a, mo_a, nbf, fock_pack)
    call unpack_matrix(fock_pack, fock_mo, nbf, 'U')

    ! ---------------------------------------------------------------------
    ! EKT-MRSF state-specific Lagrangian W~, reproduced exactly from the
    ! GAMESS-US reference implementation (sfgrad.src, the IF(MREKT) block;
    ! Pomogaev et al., JPCL 2021, eq 5):
    !
    !   W~_MO = -1/2 ( WMO + Lambda_MO )
    !
    ! where, in the MO basis and using the density metric C^T S M S C
    ! (a Lagrangian density is contravariant, like D):
    !   * Lambda = energy-weighted density Sum_k eps_k C_k C_k^T  (= eijden /
    !     GAMESS DAF record 36).  OpenQP's eijden halves the packed diagonal
    !     (a gradient-consumer convention, grd1.src:624); the GAMESS EKT path
    !     reads the raw record, so the half-diagonal is undone here.
    !   * WMO = the MRSF response Lagrangian.  OpenQP stores
    !     wao = 1/4 C WMO C^T (z-vector applies two halves), while GAMESS uses
    !     the bare WMO (records 419/429 = 1/2 C WMO C^T).  Hence
    !     WMO = 2 C^T S wao S C.
    ! ---------------------------------------------------------------------
    block
      use grd1, only: eijden
      real(kind=dp), allocatable :: lam_ao(:), lam_mo(:,:), wmo_mo(:,:)
      integer :: gok, kk, ijd
      allocate(lam_ao(nbf2), lam_mo(nbf,nbf), wmo_mo(nbf,nbf), &
               source=0.0_dp, stat=gok)
      if (gok /= 0) call show_message('Cannot allocate memory', WITH_ABORT)
      call eijden(lam_ao, nbf, infos)
      ijd = 0
      do kk = 1, nbf
        ijd = ijd + kk                       ! packed-upper diagonal index
        lam_ao(ijd) = 2.0_dp*lam_ao(ijd)     ! undo eijden's half-diagonal
      end do
      call density_ao_to_mo(nbf, nbf2, lam_ao, smat, mo_a, lam_mo)  ! C^T S Lambda S C
      call density_ao_to_mo(nbf, nbf2, wao,    smat, mo_a, wmo_mo)  ! C^T S wao S C
      wmo_mo = 2.0_dp*wmo_mo                  ! bare GAMESS WMO
      lagrangian_mo = -0.5_dp*(wmo_mo + lam_mo)
      lagrangian_mo = 0.5_dp*(lagrangian_mo + transpose(lagrangian_mo))
      deallocate(lam_ao, lam_mo, wmo_mo)
    end block

    if (electron_affinity) then
      ! EKT-EA: (F - W) * x = (I - P) * x * lambda.
      ! GAMESS rebuilds the Fock matrix for EA from the relaxed MRSF density
      ! (sfgrad.src:3269-3280, BUILDFOCK(..., MRDEA)) before forming F-W.
      ! Reproduce that path here instead of using the reference SCF Fock.
      ! GAMESS MRDEA uses the spin-averaged relaxed MRSF density
      ! P = 1/2*(Da+Db) for the EA complement metric, not beta alone.
      density_mo = 0.5_dp*(density_ip_mo + density_ea_mo)
      ekt_metric = -density_mo
      do i = 1, nbf
        ekt_metric(i,i) = ekt_metric(i,i) + 1.0_dp
      end do

      ! GAMESS MRDEA builds a closed-shell Fock from the spin-summed
      ! relaxed density: LDVAL = AO(Pavg), DSCAL(2), BUILDFOCK; internally
      ! its DFT path halves that total density for alpha=beta.
      ea_density_ao(:,1) = 0.5_dp*(density_alpha_ao + density_beta_ao)
      if (nfocks > 1) ea_density_ao(:,2) = ea_density_ao(:,1)
      saved_fock_ao(:,1) = fock_a
      if (nfocks > 1) then
        call tagarray_get_data(infos%dat, OQP_FOCK_B, fock_b)
        saved_fock_ao(:,2) = fock_b
      end if
      saved_mol_energy = infos%mol_energy

      nschwz = 0
      if (infos%control%hamilton >= 20) call dft_initialize(infos, basis, ea_molgrid)
      call calc_fock(basis, infos, ea_molgrid, ea_fock_ao, ea_energy, &
                     mo_a_in=mo_a, dens_in=ea_density_ao, nschwz=nschwz)

      ! calc_fock updates the stored SCF Fock and infos%mol_energy%energy as
      ! side effects; restore them because this EA-specific relaxed-density
      ! Fock is only an EKT operator ingredient, not a new SCF reference.
      fock_a = saved_fock_ao(:,1)
      if (nfocks > 1) fock_b = saved_fock_ao(:,2)
      infos%mol_energy = saved_mol_energy

      call orthogonal_transform_sym(nbf, nbf, ea_fock_ao(:,1), mo_a, nbf, fock_pack)
      call unpack_matrix(fock_pack, fock_mo, nbf, 'U')
      write(iw,'(/,2x,"EKT-EA: rebuilt relaxed-density Fock for EA operator (GAMESS MRDEA path)")')
      ekt_operator = fock_mo - lagrangian_mo
    else
      ! EKT: W * x = P * x * lambda.  Metric is the spin-summed relaxed
      ! density P = 1/2 (Da_MO + Db_MO) (GAMESS sfgrad.src:3115); for H2O S0
      ! this gives 5 occupied natural orbitals at occupation ~1.
      density_mo = 0.5_dp*(density_ip_mo + density_ea_mo)
      ekt_metric = density_mo
      ekt_operator = lagrangian_mo
    end if

    ! Natural-orbital occupation deflation + EKT generalized solve.
    ! The metric D_MO is non-diagonal, so a diag(D) > tol active-space
    ! selection is mathematically invalid.  Diagonalize the metric to obtain
    ! natural orbitals/occupations, keep occ > occ_tol, project W and D into
    ! the retained NO subspace, solve the generalized eigenproblem there,
    ! metric-normalize, and back-transform the Dyson orbitals.  Pole strengths
    ! come out <= 1 by construction.
    allocate(eig(nroot), metric_norms(nroot), dom_mo_coeff(nroot), dom_no_coeff(nroot), &
             dom_no_occ(nroot), source=0.0_dp, stat=ok)
    if (ok /= 0) call show_message('Cannot allocate memory', WITH_ABORT)
    allocate(dom_mo_idx(nroot), dom_no_idx(nroot), source=0, stat=ok)
    if (ok /= 0) call show_message('Cannot allocate memory', WITH_ABORT)
    ! Deflation threshold matches GAMESS EKT (thresh = 1.0d-2, sfgrad.src:3298):
    ! keep only natural orbitals with occupation > 1e-2, removing low-occupation
    ! null/relaxation channels.
    call solve_ekt_no_deflation(nbf, nroot, ekt_operator, ekt_metric, &
         ekt_occ_tol, eig, orbitals, strengths, metric_norms, dom_mo_idx, &
         dom_mo_coeff, dom_no_idx, dom_no_coeff, dom_no_occ, nactive, iw)

    call infos%dat%alloc_or_die(OQP_mrsf_ekt_density_mo, (/ nbf, nbf /), density_store, &
        description='MRSF-EKT metric density P in MO basis')
    call infos%dat%alloc_or_die(OQP_mrsf_ekt_lagrangian_mo, (/ nbf, nbf /), lagrangian_store, &
        description='MRSF-EKT Lagrangian W in MO basis')
    call infos%dat%alloc_or_die(OQP_mrsf_ekt_fock_mo, (/ nbf, nbf /), fock_store, &
        description='MRSF-EKT Fock matrix in MO basis')
    call infos%dat%alloc_or_die(OQP_mrsf_ekt_orbitals_mo, (/ nbf, nroot /), orbital_store, &
        description='MRSF-EKT Dyson-like orbital coefficients in MO basis')
    call infos%dat%alloc_or_die(OQP_mrsf_ekt_eigenvalues, (/ nroot /), eig_store, &
        description='MRSF-EKT generalized-eigenproblem eigenvalues in Hartree')
    call infos%dat%alloc_or_die(OQP_mrsf_ekt_strengths, (/ nroot /), strength_store, &
        description='MRSF-EKT pole strengths')

    density_store = density_mo
    lagrangian_store = lagrangian_mo
    fock_store = fock_mo
    orbital_store = orbitals
    eig_store = eig
    strength_store = strengths

    if (electron_affinity) then
      write(iw,'(/,2x,"MRSF-EKT electron affinities (Dyson roots)")')
    else
      write(iw,'(/,2x,"MRSF-EKT ionization potentials (Dyson roots)")')
    end if
    write(iw,'(2x,"Rows index EKT Dyson roots/orbitals, not TDDFT excited states.")')
    ! eig holds the EKT eigenvalues epsilon of  W C = D C epsilon.  The
    ! electron binding energy (IP) of a detachment is -epsilon; print both
    ! the eigenvalue and the eBE in hartree and eV with the Dyson pole
    ! strength (norm of the Dyson vector, <= 1 for a physical root).
    write(iw,'(2x,"dyson",6x,"eig(ha)",8x,"eBE(ha)",8x,"eBE(eV)",7x,"metric",7x,"strength")')
    do i = 1, min(nroot, nactive)
      write(iw,'(2x,I5,2x,F14.6,2x,F14.6,2x,F12.4,2x,F12.6,2x,F12.6)') &
            i, eig(i), -eig(i), -eig(i)*hartree_to_ev, metric_norms(i), strengths(i)
    end do
    write(iw,'(/,2x,"--- EKT Dyson-root character diagnostics ---")')
    write(iw,'(2x,"dyson",2x,"dom_mo",2x,"mo_coeff",2x,"dom_no",2x,"no_coeff",2x,"no_occ",2x,"character",2x,"spurious")')
    do i = 1, min(nroot, nactive)
      write(iw,'(2x,I5,2x,I6,2x,F10.6,2x,I6,2x,F10.6,2x,F10.6,2x,A12,2x,L1)') &
            i, dom_mo_idx(i), dom_mo_coeff(i), dom_no_idx(i), dom_no_coeff(i), &
            dom_no_occ(i), ekt_root_character(dom_no_occ(i), electron_affinity), &
            ekt_is_spurious(-eig(i)*hartree_to_ev, metric_norms(i), strengths(i))
    end do
    write(iw,'(2x,"--------------------------------------")')

    if (electron_affinity) then
      write(iw,'(/,2x,"MRSF-EKT Dyson orbitals (AO-basis coefficients, GAMESS-style layout; columns = EA roots)")')
    else
      write(iw,'(/,2x,"MRSF-EKT Dyson orbitals (AO-basis coefficients, GAMESS-style layout; columns = IP roots)")')
    end if
    call print_dyson_orbitals(basis, mo_a, orbitals, eig, strengths, nbf, min(nroot, nactive))
    call flush(iw)

    call measure_time(print_total=1, log_unit=iw)
    close(iw)

  end subroutine tdhf_mrsf_ekt

  !> @brief Print Dyson orbitals in a GAMESS-like block.
  !> @detail The EKT solver returns Dyson-orbital coefficients in the OpenQP
  !>         MO basis.  For direct log-level comparison with GAMESS, transform
  !>         them back to the AO basis and print five roots per block with
  !>         ENERGY and STRENGTH header rows followed by AO basis-function
  !>         labels and coefficients.  This routine is output-only: it does not
  !>         change eigenvalues, strengths, root selection, or stored MO-basis
  !>         Dyson orbitals.
  subroutine print_dyson_orbitals(basis, mo_a, dyson_mo, eig, strengths, nbf, nroot)
    use precision, only: dp
    use io_constants, only: iw
    use basis_tools, only: basis_set

    implicit none

    type(basis_set), intent(in) :: basis
    integer, intent(in) :: nbf, nroot
    real(kind=dp), intent(in) :: mo_a(nbf,nbf), dyson_mo(nbf,*), eig(*), strengths(*)

    integer, parameter :: mxlen = 5
    integer :: i, j, k, i0, i1
    real(kind=dp), allocatable :: dyson_ao(:,:)

    allocate(dyson_ao(nbf,nroot), source=0.0_dp)
    do i = 1, nroot
      do k = 1, nbf
        do j = 1, nbf
          dyson_ao(j,i) = dyson_ao(j,i) + mo_a(j,k)*dyson_mo(k,i)
        end do
      end do
    end do

    write(iw,'(/,10x,46("-"))')
    write(iw,'(12x,"EKT DYSON ORBITALS, ENERGIES AND NORMS")')
    write(iw,'(12x,"OpenQP AO-basis coefficients; GAMESS-style layout")')
    write(iw,'(10x,46("-"))')

    do i0 = 1, nroot, mxlen
      i1 = min(nroot, i0+mxlen-1)
      write(iw,'(/,14x,5i11)') (i, i=i0, i1)
      write(iw,'(5x,"ENERGY",3x,5f11.6)') (eig(i), i=i0, i1)
      write(iw,'(5x,"STRENGTH",1x,5f11.6)') (strengths(i), i=i0, i1)
      write(iw,'(14x,5(10x,a1))') ("A", i=i0, i1)
      do j = 1, nbf
        write(iw,'(i5,2x,a8,5f11.6)') j, basis%bf_label(j), (dyson_ao(j,i), i=i0, i1)
      end do
    end do

    deallocate(dyson_ao)

  end subroutine print_dyson_orbitals

  !> @brief Fetch the packed AO overlap matrix S (OQP_SM) into smat(nbf2).
  subroutine get_overlap_matrix(infos, nbf, smat)
    use types, only: information
    use oqp_tagarray_driver
    use precision, only: dp
    use messages, only: show_message, with_abort
    implicit none
    type(information), target, intent(inout) :: infos
    integer, intent(in) :: nbf
    real(kind=dp), intent(out) :: smat(:)
    real(kind=dp), contiguous, pointer :: s_p(:)
    character(len=*), parameter :: sn = "get_overlap_matrix"
    call data_has_tags(infos%dat, (/ character(len=80) :: OQP_SM /), &
         module_name, sn, WITH_ABORT)
    call tagarray_get_data(infos%dat, OQP_SM, s_p)
    smat(1:nbf*(nbf+1)/2) = s_p(1:nbf*(nbf+1)/2)
  end subroutine get_overlap_matrix

  !> @brief Metric-covariant AO->MO transform of a one-particle density:
  !>          D_MO = C^T (S P S) C
  !> @details A density is contravariant, so its MO (natural-orbital)
  !>          representation needs the overlap metric on both sides.  With
  !>          C^T S C = I this yields Tr(D_MO) = Tr(P S) = N_elec and MO
  !>          occupations bounded in [0,1] (per spin).  Plain C^T P C (used
  !>          for operators) is WRONG for a density.
  subroutine density_ao_to_mo(nbf, nbf2, p_ao_pack, s_pack, mo, d_mo)
    use precision, only: dp
    use mathlib, only: unpack_matrix, pack_matrix, orthogonal_transform_sym
    use messages, only: show_message, with_abort
    implicit none
    integer, intent(in) :: nbf, nbf2
    real(kind=dp), intent(in) :: p_ao_pack(:)   ! packed AO density (upper)
    real(kind=dp), intent(in) :: s_pack(:)      ! packed AO overlap (upper)
    real(kind=dp), intent(in) :: mo(nbf,nbf)    ! MO coefficients C (AO x MO)
    real(kind=dp), intent(out) :: d_mo(nbf,nbf) ! MO-basis density
    real(kind=dp), allocatable :: p_ao(:,:), s_ao(:,:), sps(:,:), tmp(:,:), sps_pack(:)
    integer :: ok
    allocate(p_ao(nbf,nbf), s_ao(nbf,nbf), sps(nbf,nbf), tmp(nbf,nbf), &
             sps_pack(nbf2), source=0.0_dp, stat=ok)
    if (ok /= 0) call show_message('Cannot allocate memory', WITH_ABORT)
    call unpack_matrix(p_ao_pack, p_ao, nbf, 'U')
    call unpack_matrix(s_pack, s_ao, nbf, 'U')
    tmp = matmul(s_ao, p_ao)        ! S P
    sps = matmul(tmp, s_ao)         ! S P S  (symmetric)
    sps = 0.5_dp*(sps + transpose(sps))
    call pack_matrix(sps, sps_pack, 'U')
    call orthogonal_transform_sym(nbf, nbf, sps_pack, mo, nbf, sps_pack)
    call unpack_matrix(sps_pack, d_mo, nbf, 'U')
    deallocate(p_ao, s_ao, sps, tmp, sps_pack)
  end subroutine density_ao_to_mo

  !> @brief Acceptance/trace checks for the EKT relaxed density (printed).
  !> @details Verifies Tr(P_sigma S) = nelec_sigma in the AO basis and that
  !>          the MO-density trace matches, with bounded MO occupations.
  subroutine ekt_density_trace_checks(infos, nbf, nbf2, pa_ao_pack, pb_ao_pack, &
       s_pack, da_mo, db_mo, iw)
    use types, only: information
    use precision, only: dp
    use mathlib, only: unpack_matrix
    use messages, only: show_message, with_abort
    implicit none
    type(information), target, intent(inout) :: infos
    integer, intent(in) :: nbf, nbf2, iw
    real(kind=dp), intent(in) :: pa_ao_pack(:), pb_ao_pack(:), s_pack(:)
    real(kind=dp), intent(in) :: da_mo(nbf,nbf), db_mo(nbf,nbf)
    real(kind=dp), allocatable :: pa(:,:), pb(:,:), s_ao(:,:)
    real(kind=dp) :: trPaS, trPbS, trDa, trDb, occ_min, occ_max
    integer :: i, j, nea, neb, ok
    allocate(pa(nbf,nbf), pb(nbf,nbf), s_ao(nbf,nbf), source=0.0_dp, stat=ok)
    if (ok /= 0) call show_message('Cannot allocate memory', WITH_ABORT)
    call unpack_matrix(pa_ao_pack, pa, nbf, 'U')
    call unpack_matrix(pb_ao_pack, pb, nbf, 'U')
    call unpack_matrix(s_pack, s_ao, nbf, 'U')
    nea = int(infos%mol_prop%nelec_A)
    neb = int(infos%mol_prop%nelec_B)
    trPaS = 0.0_dp; trPbS = 0.0_dp
    do i = 1, nbf
      do j = 1, nbf
        trPaS = trPaS + pa(i,j)*s_ao(j,i)
        trPbS = trPbS + pb(i,j)*s_ao(j,i)
      end do
    end do
    trDa = 0.0_dp; trDb = 0.0_dp; occ_min = 1.0d30; occ_max = -1.0d30
    do i = 1, nbf
      trDa = trDa + da_mo(i,i)
      trDb = trDb + db_mo(i,i)
      occ_min = min(occ_min, da_mo(i,i), db_mo(i,i))
      occ_max = max(occ_max, da_mo(i,i), db_mo(i,i))
    end do
    write(iw,'(/,2x,"--- EKT relaxed-density acceptance checks ---")')
    write(iw,'(2x,"Tr(Pa S) = ",f12.6,"   Tr(Da_MO) = ",f12.6)') trPaS, trDa
    write(iw,'(2x,"Tr(Pb S) = ",f12.6,"   Tr(Db_MO) = ",f12.6)') trPbS, trDb
    write(iw,'(2x,"(AO and MO traces must agree; target-state singlet S0)")')
    write(iw,'(2x,"MO occupations range: [",f10.6,",",f10.6,"]")') occ_min, occ_max
    write(iw,'(2x,"--------------------------------------------",/)')
    call flush(iw)
    deallocate(pa, pb, s_ao)
  end subroutine ekt_density_trace_checks

  !> @brief EKT solve with natural-orbital occupation deflation.
  !> @details Solves  W X = D X eps  where the metric D (= ekt_metric) is the
  !>          relaxed one-particle density (IP) or its particle complement (EA),
  !>          and W (= ekt_operator) is the EKT Lagrangian.  Because D is
  !>          non-diagonal, deflation must be done on its NATURAL ORBITALS, not
  !>          its diagonal:
  !>            1. symmetrize D
  !>            2. diagonalize  D U = U n   (natural orbitals U, occupations n)
  !>            3. keep NOs with n_i > occ_tol
  !>            4. project  D_NO = U_keep^T D U_keep,  W_NO = U_keep^T W U_keep
  !>            5. solve the generalized problem in the retained NO space
  !>            6. metric-normalize eigenvectors:  X^T D_NO X = 1
  !>            7. back-transform Dyson orbitals  C = U_keep X
  !>            8. pole strength = X^T D_NO X  (<= 1 physical)
  subroutine solve_ekt_no_deflation(nbf, nroot, wmat, dmat, occ_tol, &
       eig, dyson_mo, strengths, metric_norms, dom_mo_idx, dom_mo_coeff, &
       dom_no_idx, dom_no_coeff, dom_no_occ, nkeep, iw)
    use precision, only: dp
    use eigen, only: diag_symm_full
    use messages, only: show_message, with_abort
    implicit none
    integer, intent(in) :: nbf, nroot, iw
    real(kind=dp), intent(in) :: wmat(nbf,nbf), dmat(nbf,nbf), occ_tol
    real(kind=dp), intent(out) :: eig(nroot)
    real(kind=dp), intent(out) :: dyson_mo(nbf,nroot)
    real(kind=dp), intent(out) :: strengths(nroot)
    real(kind=dp), intent(out) :: metric_norms(nroot)
    integer, intent(out) :: dom_mo_idx(nroot), dom_no_idx(nroot)
    real(kind=dp), intent(out) :: dom_mo_coeff(nroot), dom_no_coeff(nroot), dom_no_occ(nroot)
    integer, intent(out) :: nkeep
    real(kind=dp), allocatable :: dsym(:,:), nocc(:), uno(:,:), ukeep(:,:)
    integer, allocatable :: kept_idx(:)
    real(kind=dp), allocatable :: d_no(:,:), w_no(:,:), tmp(:,:)
    real(kind=dp), allocatable :: xvec(:,:), eps(:), cdys(:,:), d_times_x(:)
    real(kind=dp) :: tr_disc, occ_min, occ_max, xnorm, str, coeff_abs, best_abs
    integer :: i, j, m, ierr, ok, nout, best_idx, out_idx, n_skipped
    ! EKT eigenvalue print/retention threshold.  GAMESS skips near-zero EKT
    ! eigenvalues (|eps| <= thresh) when transferring solved EKT roots to its
    ! output (sfgrad.src:3411, thresh = 1.0d-2 Ha); mirror that here so the
    ! printed/stored Dyson ladder aligns with GAMESS without a one-root shift.
    ! Kept separate from the NO occupation tolerance (ekt_occ_tol) which only
    ! governs natural-orbital deflation.
    real(kind=dp), parameter :: ekt_eig_tol = 1.0e-2_dp

    ! (1) symmetrize the metric
    allocate(dsym(nbf,nbf), nocc(nbf), uno(nbf,nbf), source=0.0_dp, stat=ok)
    if (ok /= 0) call show_message('Cannot allocate memory', WITH_ABORT)
    do j = 1, nbf
      do i = 1, nbf
        dsym(i,j) = 0.5_dp*(dmat(i,j) + dmat(j,i))
      end do
    end do

    ! (2) diagonalize the metric -> natural orbitals U, occupations n
    uno = dsym
    call diag_symm_full(1, nbf, uno, nbf, nocc, ierr)
    if (ierr /= 0) call show_message('EKT: NO diagonalization failed', WITH_ABORT)

    ! (3) keep NOs with physical occupation n_i > occ_tol
    m = 0
    tr_disc = 0.0_dp
    do i = 1, nbf
      if (nocc(i) > occ_tol) then
        m = m + 1
      else
        tr_disc = tr_disc + nocc(i)
      end if
    end do
    if (m == 0) call show_message('EKT: no NOs above occupation tolerance', WITH_ABORT)

    allocate(ukeep(nbf,m), source=0.0_dp, stat=ok)
    if (ok /= 0) call show_message('Cannot allocate memory', WITH_ABORT)
    allocate(kept_idx(m), source=0, stat=ok)
    if (ok /= 0) call show_message('Cannot allocate memory', WITH_ABORT)
    j = 0
    occ_min = 1.0d30; occ_max = -1.0d30
    do i = 1, nbf
      if (nocc(i) > occ_tol) then
        j = j + 1
        ukeep(:,j) = uno(:,i)
        kept_idx(j) = i
        occ_min = min(occ_min, nocc(i))
        occ_max = max(occ_max, nocc(i))
      end if
    end do

    ! (4) deflation diagnostics
    write(iw,'(/,2x,"--- EKT natural-orbital deflation ---")')
    write(iw,'(2x,"occ_tol = ",es10.2,"   retained NOs = ",i0," / ",i0)') occ_tol, m, nbf
    write(iw,'(2x,"discarded trace = ",f12.6,"   retained occ range = [",f10.6,",",f10.6,"]")') &
          tr_disc, occ_min, occ_max
    write(iw,'(2x,"eig_tol = ",es10.2," Ha   (skip EKT roots with |eps| <= eig_tol; GAMESS sfgrad.src)")') &
          ekt_eig_tol
    write(iw,'(2x,"-------------------------------------",/)')

    ! (5) project D and W into the retained NO space
    allocate(d_no(m,m), w_no(m,m), tmp(nbf,m), source=0.0_dp, stat=ok)
    if (ok /= 0) call show_message('Cannot allocate memory', WITH_ABORT)
    tmp  = matmul(dsym, ukeep)
    d_no = matmul(transpose(ukeep), tmp)
    tmp  = matmul(wmat, ukeep)
    w_no = matmul(transpose(ukeep), tmp)
    d_no = 0.5_dp*(d_no + transpose(d_no))
    w_no = 0.5_dp*(w_no + transpose(w_no))

    ! (6) solve the generalized problem in the retained NO space
    allocate(xvec(m,m), eps(m), source=0.0_dp, stat=ok)
    if (ok /= 0) call show_message('Cannot allocate memory', WITH_ABORT)
    call solve_symmetric_generalized(w_no, d_no, eps, xvec, m, ierr)
    if (ierr /= 0) call show_message('(A,I0)', &
         'EKT: NO-space generalized solve failed; info=', ierr, WITH_ABORT)

    ! (7) metric-normalize eigenvectors:  X^T D_NO X = 1
    do j = 1, m
      xnorm = 0.0_dp
      do i = 1, m
        xnorm = xnorm + xvec(i,j)*dot_product(d_no(i,:), xvec(:,j))
      end do
      if (xnorm > 1.0e-14_dp) xvec(:,j) = xvec(:,j)/sqrt(xnorm)
    end do

    ! (8) back-transform Dyson orbitals to the MO basis  C = U_keep X.
    !     Metric norms are X^T D_NO X after normalization (=1 for retained
    !     physical roots).  EKT pole strengths follow the GAMESS/Pomogaev
    !     definition X^T D_NO^2 X; weak Dyson roots therefore retain small
    !     strengths instead of being forced to 1 by metric normalization.
    allocate(cdys(nbf,m), d_times_x(m), source=0.0_dp, stat=ok)
    if (ok /= 0) call show_message('Cannot allocate memory', WITH_ABORT)
    cdys = matmul(ukeep, xvec)

    eig = 0.0_dp
    dyson_mo = 0.0_dp
    strengths = 0.0_dp
    metric_norms = 0.0_dp
    dom_mo_idx = 0
    dom_no_idx = 0
    dom_mo_coeff = 0.0_dp
    dom_no_coeff = 0.0_dp
    dom_no_occ = 0.0_dp
    ! Transfer solved EKT roots to the output arrays, skipping near-zero EKT
    ! eigenvalues the GAMESS way (|eps| <= ekt_eig_tol).  The source eigenpair
    ! index is j (xvec/cdys/eps); out_idx is a separate compressed counter into
    ! the output arrays so the printed/stored ladder matches GAMESS.
    out_idx = 0
    n_skipped = 0
    do j = 1, m
      if (abs(eps(j)) <= ekt_eig_tol) then
        n_skipped = n_skipped + 1
        cycle
      end if
      out_idx = out_idx + 1
      if (out_idx > nroot) exit
      eig(out_idx) = eps(j)
      dyson_mo(:,out_idx) = cdys(:,j)
      str = 0.0_dp
      do i = 1, m
        str = str + xvec(i,j)*dot_product(d_no(i,:), xvec(:,j))
      end do
      metric_norms(out_idx) = str
      d_times_x = matmul(d_no, xvec(:,j))
      strengths(out_idx) = dot_product(d_times_x, d_times_x)

      best_abs = -1.0_dp
      best_idx = 0
      do i = 1, nbf
        coeff_abs = abs(cdys(i,j))
        if (coeff_abs > best_abs) then
          best_abs = coeff_abs
          best_idx = i
        end if
      end do
      dom_mo_idx(out_idx) = best_idx
      if (best_idx > 0) dom_mo_coeff(out_idx) = cdys(best_idx,j)

      best_abs = -1.0_dp
      best_idx = 0
      do i = 1, m
        coeff_abs = abs(xvec(i,j))
        if (coeff_abs > best_abs) then
          best_abs = coeff_abs
          best_idx = i
        end if
      end do
      if (best_idx > 0) then
        dom_no_idx(out_idx) = kept_idx(best_idx)
        dom_no_coeff(out_idx) = xvec(best_idx,j)
        dom_no_occ(out_idx) = nocc(kept_idx(best_idx))
      end if
    end do
    nout = min(out_idx, nroot)
    ! The caller uses nkeep as nactive when printing rows/Dyson orbitals, so it
    ! must reflect the number of stored output roots (nout), not the retained
    ! NO-space dimension (m); otherwise a zero row is printed for skipped roots.
    nkeep = nout
    write(iw,'(2x,"EKT eigenvalue filter: skipped ",i0," near-zero root(s) ", &
         &"(|eps| <= ",es10.2," Ha); stored ",i0," / ",i0," roots")') &
          n_skipped, ekt_eig_tol, nout, nroot

    deallocate(dsym, nocc, uno, ukeep, kept_idx, d_no, w_no, tmp, xvec, eps, cdys, d_times_x)
  end subroutine solve_ekt_no_deflation

  pure logical function ekt_is_spurious(ebe_ev, metric_norm, pole_strength) result(flag)
    use precision, only: dp
    implicit none
    real(kind=dp), intent(in) :: ebe_ev, metric_norm, pole_strength
    flag = (abs(ebe_ev) > 1.0e4_dp) .or. (metric_norm < 0.0_dp) .or. &
           (metric_norm > 1.05_dp) .or. (pole_strength < 0.0_dp) .or. &
           (pole_strength > 1.05_dp)
  end function ekt_is_spurious

  pure function ekt_root_character(no_occ, electron_affinity) result(label)
    use precision, only: dp
    implicit none
    real(kind=dp), intent(in) :: no_occ
    logical, intent(in) :: electron_affinity
    character(len=12) :: label
    if (electron_affinity) then
      if (no_occ < 0.20_dp) then
        label = 'virtual'
      else if (no_occ > 0.80_dp) then
        label = 'occupied'
      else
        label = 'mixed'
      end if
    else
      if (no_occ > 0.80_dp) then
        label = 'occupied'
      else if (no_occ < 0.20_dp) then
        label = 'virtual'
      else
        label = 'mixed'
      end if
    end if
  end function ekt_root_character

  !> @brief Solve the EKT generalized eigenproblem  op * C = metric * C * lambda
  !> @details The EKT working equation (eq 1 of Park et al., JCTC 2024) is a
  !>          symmetric generalized eigenproblem in which the metric is the
  !>          relaxed one-particle density (IP) or its particle complement (EA).
  !>          That metric is symmetric positive semidefinite but in general NOT
  !>          diagonal: the MRSF Z-vector relaxation introduces off-diagonal
  !>          occupation couplings in the ground-state MO basis.  We therefore
  !>          reduce the problem by Loewdin (symmetric) orthogonalization built
  !>          from the full eigendecomposition of the metric, discarding the
  !>          null-space directions whose occupation eigenvalues fall below a
  !>          tolerance.  This is the rigorous replacement for a diagonal-only
  !>          metric scaling, which is exact only when the metric is diagonal.
  subroutine solve_symmetric_generalized(operator, metric, eigenvalues, eigenvectors, n, ierr)
    use precision, only: dp
    use eigen, only: diag_symm_full
    use messages, only: show_message, with_abort
    implicit none
    integer, intent(in) :: n
    real(kind=dp), intent(inout) :: operator(n,n)
    real(kind=dp), intent(in) :: metric(n,n)
    real(kind=dp), intent(out) :: eigenvalues(n), eigenvectors(n,n)
    integer, intent(out) :: ierr
    integer :: i, j, m, ok
    real(kind=dp), parameter :: metric_eval_tol = 1.0e-10_dp
    real(kind=dp), allocatable :: smat(:,:), seval(:), xorth(:,:)
    real(kind=dp), allocatable :: opsym(:,:), tmp(:,:), reduced(:,:)
    real(kind=dp), allocatable :: vec(:,:), eval_red(:)

    ierr = 0
    eigenvalues = 0.0_dp
    eigenvectors = 0.0_dp
    if (n <= 0) return

    ! Eigendecomposition of the (symmetrized) metric:  S = U diag(seval) U^T
    allocate(smat(n,n), seval(n), source=0.0_dp, stat=ok)
    if (ok /= 0) call show_message('Cannot allocate memory', WITH_ABORT)
    do j = 1, n
      do i = 1, n
        smat(i,j) = 0.5_dp*(metric(i,j) + metric(j,i))
      end do
    end do
    call diag_symm_full(1, n, smat, n, seval, ierr)
    if (ierr /= 0) return

    ! Keep only directions with positive occupation (rank of the metric).
    m = 0
    do i = 1, n
      if (seval(i) > metric_eval_tol) m = m + 1
    end do
    if (m == 0) then
      ierr = -1
      return
    end if

    ! Symmetric orthogonalizer  X = U_kept * diag(seval_kept)^(-1/2)   (n x m)
    allocate(xorth(n,m), source=0.0_dp, stat=ok)
    if (ok /= 0) call show_message('Cannot allocate memory', WITH_ABORT)
    j = 0
    do i = 1, n
      if (seval(i) > metric_eval_tol) then
        j = j + 1
        xorth(:,j) = smat(:,i) / sqrt(seval(i))
      end if
    end do

    ! Reduced operator in the orthonormal metric basis:  reduced = X^T op X
    allocate(opsym(n,n), tmp(n,m), reduced(m,m), source=0.0_dp, stat=ok)
    if (ok /= 0) call show_message('Cannot allocate memory', WITH_ABORT)
    do j = 1, n
      do i = 1, n
        opsym(i,j) = 0.5_dp*(operator(i,j) + operator(j,i))
      end do
    end do
    tmp = matmul(opsym, xorth)
    reduced = matmul(transpose(xorth), tmp)
    reduced = 0.5_dp*(reduced + transpose(reduced))

    ! Standard symmetric eigenproblem on the reduced operator.
    allocate(vec(m,m), eval_red(m), source=0.0_dp, stat=ok)
    if (ok /= 0) call show_message('Cannot allocate memory', WITH_ABORT)
    vec = reduced
    call diag_symm_full(1, m, vec, m, eval_red, ierr)
    if (ierr /= 0) return

    ! Map results back:  eigenvalues are lambda, eigenvectors C = X * vec.
    eigenvalues(1:m) = eval_red(1:m)
    eigenvectors(:,1:m) = matmul(xorth, vec)
  end subroutine solve_symmetric_generalized

end module tdhf_mrsf_ekt_mod
