module nmr_giao_shielding_mod

  use precision, only: dp
  implicit none

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

  private
  public nmr_giao_shielding_debug

contains

  subroutine nmr_giao_shielding_debug_C(c_handle) bind(C, name="nmr_giao_shielding_debug")
    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 nmr_giao_shielding_debug(inf)
  end subroutine nmr_giao_shielding_debug_C

  !> Production entry point for GIAO NMR shielding (RHF/UHF/ROHF).
  subroutine nmr_giao_shielding_C(c_handle) bind(C, name="nmr_giao_shielding")
    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 nmr_giao_shielding_debug(inf)
  end subroutine nmr_giao_shielding_C

!> @brief Test-only/debug emitter for native GIAO (London-orbital) NMR shielding.
!> @details Computes the GIAO paramagnetic nuclear magnetic shielding tensor for
!>  RHF/UHF/ROHF (and the corresponding DFT) references.  Open-shell support:
!>  the diamagnetic term uses the total density D_alpha+D_beta; the paramagnetic
!>  term is solved per spin channel (independent same-spin exchange response) and
!>  summed.  ROHF orbitals are semicanonicalized (from the ground-state spin Fock
!>  matrices) so the UHF-like CPHF is well defined; the closed-shell ROHF limit
!>  reproduces RHF, and open-shell UHF/ROHF match an independent GIAO reference to ~1e-4
!>  ppm (OH, CH3 radicals).  Built from the validated native GIAO building blocks:
!>   - first-order GIAO core Hamiltonian h10 (one-electron, giao_h10_core),
!>   - first-order GIAO two-electron Fock derivative (giao_h10_twoe_matrix),
!>   - first-order GIAO overlap derivative S10 (giao_overlap_derivative),
!>   - the PSO operator at each nucleus (pso_integrals).
!>  It assembles the magnetic first-order Hamiltonian/overlap in the MO basis,
!>  solves the GIAO CPHF/CPKS first-order equation (uncoupled and coupled; the
!>  coupled response is exchange-only for the imaginary/antisymmetric first-order
!>  density, scaled by the exact-exchange fraction c_x), and contracts the
!>  resulting first-order density with the PSO operator to form the paramagnetic
!>  shielding.
!>
!>  Hartree-Fock (RHF/UHF/ROHF) is validated against an independent GIAO reference.
!>
!>  DFT (RKS/UKS/ROKS) is supported: the closed-shell h1 two-electron exchange is
!>  scaled by c_x, and the London (GIAO) derivative of the exchange-correlation
!>  potential (the "vxc_giao" term -- an essential contribution, ~tens of ppm) is
!>  added via a GIAO-weighted XC grid integration (mod_dft_gridint_giao::giao_vxc;
!>  the London derivative of the AO values and gradients).  Validated vs an independent
!>  reference (replicating get_vxc_giao): H2O/PBE and /PBE0 reproduce the
!>  reference to ~1e-3 ppm.  vxc_giao is grid-sensitive, so DFT GIAO benefits from
!>  a converged integration grid.
!>
!>  Results are written as machine-parseable records to the log so
!>  the native GIAO path can be validated against an independent GIAO reference WITHOUT
!>  ungating production nmr_gauge=giao.
!>
!>  Diamagnetic GIAO term (both pieces reference-validated):
!>   - a11part (London diamagnetic): validated CGO diamagnetic at gauge origin 0
!>     plus the Hellmann-Feynman field correction weighted by the ket center.
!>   - a01gp (GIAO gauge correction = London derivative of the PSO): cvec x M,
!>     M^{(col)}_b = <mu| r_b PSO_col |nu> with r referenced to the molecular
!>     origin (R0I = (r-R_bra) bra-raise + R_bra*base, the libcint convention).
!>     Full-tensor agreement with libcint int1e_a01gp to ~3e-8 (e.g. CH4).
!>  Total GIAO shielding matches an independent GIAO reference to ~1e-4 ppm for HF
!>  (grid-free) and to cross-code DFT-SCF/grid level (~0.03 ppm) for DFT, for all
!>  geometries and angular momenta tested (He, H2, HF, CO2, CH4, H2O).
!>  SG/SC/SA are sign/scale conventions fixed against the oracle (SA folds the
!>  factor-2 normalization of the native a01gp vs the libcint int1e_a01gp).
!>
!>  References (GIAO/London-orbital NMR shielding methodology):
!>   - F. London, J. Phys. Radium 8, 397 (1937).
!>   - R. Ditchfield, Mol. Phys. 27, 789 (1974).
!>   - K. Wolinski, J. F. Hinton, P. Pulay, J. Am. Chem. Soc. 112, 8251 (1990).
!>   - T. Helgaker, M. Jaszunski, K. Ruud, Chem. Rev. 99, 293 (1999).
  subroutine nmr_giao_shielding_debug(infos)
    use io_constants, only: iw
    use oqp_tagarray_driver
    use basis_tools, only: basis_set
    use messages, only: show_message, with_abort
    use types, only: information
    use constants, only: tol_int
    use int1, only: giao_h10_core, giao_overlap_derivative, pso_integrals, &
                    nmr_dia_shielding, giao_a11part_corr, giao_a01gp_contract
    use nmr_giao_debug_mod, only: giao_h10_twoe_matrix
    use dft, only: dft_initialize
    use mod_dft_molgrid, only: dft_grid_t
    use mod_dft_gridint_giao, only: giao_vxc

    implicit none

    character(len=*), parameter :: subroutine_name = "nmr_giao_shielding_debug"
    real(kind=dp), parameter :: ALPHA = 1.0d0/137.035999084d0
    real(kind=dp), parameter :: a2ppm = ALPHA*ALPHA*1.0d6
    real(kind=dp), parameter :: ha2ppm = 0.5d0*ALPHA*ALPHA*1.0d6
    ! Calibration signs for the GIAO diamagnetic pieces (fixed vs the reference).
    real(kind=dp), parameter :: SG = -1.0d0, SC = 1.0d0, SA = -0.5d0
    real(kind=dp), parameter :: SX = -1.0d0   ! London-XC (vxc_giao): h1 -= vxc_giao

    type(information), target, intent(inout) :: infos
    type(basis_set), pointer :: basis

    integer :: nbf, nbf2, nat, nocc, nmo, nvir, nocc_b
    integer :: i, j, m, c, t, s, ok, iat
    integer(4) :: status
    logical :: is_dft, open_shell, iw_open
    real(kind=dp) :: tol, scale_exch

    real(kind=dp), allocatable :: h10p(:,:), s10p(:,:)            ! packed (nbf2,3)
    real(kind=dp), allocatable :: twoe(:,:,:), twoe2(:,:,:), vj(:,:,:), vk(:,:,:), vkb(:,:,:)
    real(kind=dp), allocatable :: dm(:,:), dmp(:), dm_b(:,:)
    real(kind=dp), allocatable :: h1ao(:,:,:), h1ao_b(:,:,:), s1ao(:,:,:) ! (nbf,nbf,3)
    real(kind=dp), allocatable :: coords(:,:), zq(:)
    real(kind=dp), allocatable :: sig_u(:,:,:), sig_c(:,:,:)      ! (3,3,nat)
    real(kind=dp), allocatable :: vxa(:,:,:), vxb(:,:,:)          ! London-XC (vxc_giao)
    type(dft_grid_t) :: molGrid
    integer :: mxAngMom
    real(kind=dp), allocatable :: gdia0(:,:,:), corrpre(:,:,:)    ! GIAO dia pieces
    real(kind=dp), allocatable :: a01(:,:,:)                      ! a01gp contracted
    real(kind=dp), allocatable :: sig_dia(:,:,:), sig_tot(:,:,:)  ! (3,3,nat)
    real(kind=dp) :: trg0, trc, o0(3)

    real(kind=dp), contiguous, pointer :: dmat_a(:), dmat_b(:)
    real(kind=dp), contiguous, pointer :: mo_a(:,:), mo_b(:,:)
    real(kind=dp), contiguous, pointer :: mo_e(:), mo_e_b(:)
    real(kind=dp), contiguous, pointer :: fock_a(:), fock_b(:)
    real(kind=dp), contiguous, pointer :: nmrout(:)
    real(kind=dp), allocatable :: ca_sc(:,:), cb_sc(:,:), ea_sc(:), eb_sc(:)

    basis => infos%basis
    basis%atoms => infos%atoms
    nbf = basis%nbf
    nbf2 = nbf*(nbf+1)/2
    nat = ubound(basis%atoms%zn,1)
    tol = log(10.0d0)*tol_int

    ! Connect the log unit early so guard aborts and CPHF warnings land in the
    ! log instead of an orphan fort.* file.
    inquire(unit=iw, opened=iw_open)
    if (.not. iw_open) open(unit=iw, file=infos%log_filename, position="append")

    call tagarray_get_data(infos%dat, OQP_DM_A, dmat_a, status)
    call check_status(status, module_name, subroutine_name, OQP_DM_A)
    call tagarray_get_data(infos%dat, OQP_VEC_MO_A, mo_a, status)
    call check_status(status, module_name, subroutine_name, OQP_VEC_MO_A)
    call tagarray_get_data(infos%dat, OQP_E_MO_A, mo_e, status)
    call check_status(status, module_name, subroutine_name, OQP_E_MO_A)

    ! scftype: 1=RHF (closed shell), 2=UHF, 3=ROHF
    open_shell = infos%control%scftype == 2 .or. infos%control%scftype == 3
    nmo  = size(mo_e)
    if (open_shell) then
      nocc   = int(infos%mol_prop%nelec_A)
      nocc_b = int(infos%mol_prop%nelec_B)
      call tagarray_get_data(infos%dat, OQP_DM_B, dmat_b, status)
      call check_status(status, module_name, subroutine_name, OQP_DM_B)
      call tagarray_get_data(infos%dat, OQP_VEC_MO_B, mo_b, status)
      call check_status(status, module_name, subroutine_name, OQP_VEC_MO_B)
      call tagarray_get_data(infos%dat, OQP_E_MO_B, mo_e_b, status)
      call check_status(status, module_name, subroutine_name, OQP_E_MO_B)
    else
      nocc   = int(infos%mol_prop%nocc)
      nocc_b = nocc
    end if
    nvir = nmo - nocc

    is_dft = infos%control%hamilton == 20
    scale_exch = 1.0d0
    if (is_dft) scale_exch = infos%dft%HFscale

    ! Not-implemented classes must abort instead of silently producing wrong
    ! shieldings:
    !  - CAM/range-separated hybrids: the coupled magnetic response and the
    !    GIAO two-electron derivative use the global exchange fraction only;
    !    the range-separation attenuation (alpha/beta/mu) is not wired in.
    !  - meta-GGAs: the tau channel of the London-XC (vxc_giao) term is not
    !    implemented (mod_dft_gridint_giao evaluates LDA/GGA ingredients only).
    !  - ECP: no effective-core magnetic-derivative term is implemented.
    if (is_dft) then
      if (infos%dft%cam_flag) then
        call show_message('GIAO NMR shielding with range-separated (CAM) &
          &functionals is not implemented', with_abort)
      end if
      if (infos%functional%needtau) then
        call show_message('GIAO NMR shielding with meta-GGA (tau-dependent) &
          &functionals is not implemented', with_abort)
      end if
    end if
    if (allocated(infos%basis%ecp_zn_num)) then
      if (any(infos%basis%ecp_zn_num /= 0)) then
        call show_message('GIAO NMR shielding with ECP basis sets is not &
          &implemented', with_abort)
      end if
    end if

    allocate(coords(3,nat), zq(nat))
    do iat = 1, nat
      coords(:,iat) = basis%atoms%xyz(:,iat)
    end do
    zq = infos%atoms%zn - infos%basis%ecp_zn_num

    ! --- Total density (full + packed).  For RHF OQP_DM_A is already the total
    !     (closed-shell) density; for UHF/ROHF total = D_alpha + D_beta. ---
    allocate(dm(nbf,nbf), dmp(nbf2), source=0.0d0)
    if (open_shell) then
      allocate(dm_b(nbf,nbf), source=0.0d0)
      dmp = dmat_a + dmat_b
      call unpack_sym(dmat_a, dm, nbf)        ! D_alpha
      call unpack_sym(dmat_b, dm_b, nbf)      ! D_beta
    else
      dmp = dmat_a
      call unpack_sym(dmp, dm, nbf)           ! D_total (closed shell)
    end if

    ! --- One-electron GIAO Hamiltonian h10(1e) and overlap derivative S10 ---
    allocate(h10p(nbf2,3), s10p(nbf2,3), source=0.0d0)
    call giao_h10_core(basis, coords, zq, h10p, debug=.false., logtol=tol)
    call giao_overlap_derivative(basis, s10p, debug=.false., logtol=tol)

    allocate(h1ao(nbf,nbf,3), s1ao(nbf,nbf,3), source=0.0d0)
    do c = 1, 3
      call expand_antisym(h10p(:,c), h1ao(:,:,c), nbf)
      call expand_antisym(s10p(:,c), s1ao(:,:,c), nbf)
    end do

    ! --- London (GIAO) derivative of the XC potential (vxc_giao); DFT only ---
    allocate(vxa(3,nbf,nbf), vxb(3,nbf,nbf), source=0.0d0)
    if (is_dft) then
      mxAngMom = maxval(basis%am) + 2
      call dft_initialize(infos, basis, molGrid)
      block
        real(kind=dp), allocatable :: ca(:,:), cb(:,:)
        allocate(ca(nbf,nbf), cb(nbf,nbf))
        ca = mo_a(:,1:nbf)
        if (open_shell) then
          cb = mo_b(:,1:nbf)
          call giao_vxc(basis, molGrid, infos, ca, cb, .true., &
                        vxa, vxb, mxAngMom, nbf, infos%dft%grid_density_cutoff)
        else
          cb = ca
          call giao_vxc(basis, molGrid, infos, ca, cb, .false., &
                        vxa, vxb, mxAngMom, nbf, infos%dft%grid_density_cutoff)
        end if
        deallocate(ca, cb)
      end block
    end if

    ! --- Two-electron GIAO Fock derivative ---
    allocate(twoe(3,nbf,nbf), twoe2(3,nbf,nbf), vj(3,nbf,nbf), vk(3,nbf,nbf), source=0.0d0)
    allocate(sig_u(3,3,nat), sig_c(3,3,nat), source=0.0d0)
    if (open_shell) then
      ! Spin-resolved: h1_sigma = h10(1e) + J[D_tot] - cx*K[D_sigma].  giao_h10_
      ! twoe_matrix returns (vj=J, vk=K, h10) for its input density; call it once
      ! per density, using a throwaway 'twoe' for the unused J/h10 outputs.
      allocate(vkb(3,nbf,nbf), h1ao_b(nbf,nbf,3), source=0.0d0)
      call giao_h10_twoe_matrix(basis, infos, dm+dm_b, vj,  twoe, twoe2) ! vj  = J[D_tot]
      call giao_h10_twoe_matrix(basis, infos, dm,      twoe, vk,  twoe2) ! vk  = K[D_a]
      call giao_h10_twoe_matrix(basis, infos, dm_b,    twoe, vkb, twoe2) ! vkb = K[D_b]
      do c = 1, 3
        do i = 1, nbf
          do j = 1, nbf
            h1ao_b(i,j,c) = h1ao(i,j,c) + vj(c,i,j) - scale_exch*vkb(c,i,j) + SX*vxb(c,i,j)
            h1ao(i,j,c)   = h1ao(i,j,c) + vj(c,i,j) - scale_exch*vk(c,i,j)  + SX*vxa(c,i,j)
          end do
        end do
      end do
      ! Two independent spin channels (same-spin exchange), each occ_factor = 1.
      if (infos%control%scftype == 3) then
        ! ROHF: a single orbital set with effective-Fock eigenvalues.  Build
        ! semicanonical spin orbitals/energies from the ground-state spin Fock
        ! matrices so the UHF-like CPHF/Delta-e is well defined.
        call tagarray_get_data(infos%dat, OQP_FOCK_A, fock_a, status)
        call check_status(status, module_name, subroutine_name, OQP_FOCK_A)
        call tagarray_get_data(infos%dat, OQP_FOCK_B, fock_b, status)
        call check_status(status, module_name, subroutine_name, OQP_FOCK_B)
        allocate(ca_sc(nbf,nmo), cb_sc(nbf,nmo), ea_sc(nmo), eb_sc(nmo), source=0.0d0)
        call semicanon_orbitals(fock_a, mo_a, nbf, nmo, nocc,   ca_sc, ea_sc)
        call semicanon_orbitals(fock_b, mo_b, nbf, nmo, nocc_b, cb_sc, eb_sc)
        call giao_para_channel(infos, basis, ca_sc, ea_sc, nocc,   nmo, nbf, nat, &
                               coords, h1ao,   s1ao, scale_exch, 1.0d0, sig_u, sig_c)
        call giao_para_channel(infos, basis, cb_sc, eb_sc, nocc_b, nmo, nbf, nat, &
                               coords, h1ao_b, s1ao, scale_exch, 1.0d0, sig_u, sig_c)
      else
        call giao_para_channel(infos, basis, mo_a, mo_e,   nocc,   nmo, nbf, nat, &
                               coords, h1ao,   s1ao, scale_exch, 1.0d0, sig_u, sig_c)
        call giao_para_channel(infos, basis, mo_b, mo_e_b, nocc_b, nmo, nbf, nat, &
                               coords, h1ao_b, s1ao, scale_exch, 1.0d0, sig_u, sig_c)
      end if
    else
      ! Closed shell: h1 = h10(1e) + (J - 0.5 K)[D_tot] (twoe), single channel.
      call giao_h10_twoe_matrix(basis, infos, dm, vj, vk, twoe)
      do c = 1, 3
        do i = 1, nbf
          do j = 1, nbf
            h1ao(i,j,c) = h1ao(i,j,c) + vj(c,i,j) - 0.5d0*scale_exch*vk(c,i,j) + SX*vxa(c,i,j)
          end do
        end do
      end do
      call giao_para_channel(infos, basis, mo_a, mo_e, nocc, nmo, nbf, nat, &
                             coords, h1ao, s1ao, scale_exch, 2.0d0, sig_u, sig_c)
    end if
    sig_u = sig_u * a2ppm
    sig_c = sig_c * a2ppm

    ! --- Diamagnetic shielding (GIAO) ---
    !   a11part = cg_a11part(O=0) + 0.5 field_a R_nu,b  (verified vs libcint).
    !   e11_pre_{t,s} = 0.5*gdia0_{s,t} + corrpre_{t,s};  e11 = e11_pre - I*tr;
    !   sigma_dia = e11 * alpha^2 * 1e6.  (a01gp gauge-correction: TODO.)
    o0 = 0.0d0
    allocate(gdia0(3,3,nat), corrpre(3,3,nat), a01(3,3,nat), &
             sig_dia(3,3,nat), sig_tot(3,3,nat), source=0.0d0)
    call nmr_dia_shielding(basis, dmp, o0, coords, nat, gdia0, logtol=tol)
    call giao_a11part_corr(basis, dmp, coords, nat, corrpre, logtol=tol)
    call giao_a01gp_contract(basis, dmp, coords, nat, a01, logtol=tol)
    do iat = 1, nat
      trg0 = gdia0(1,1,iat)+gdia0(2,2,iat)+gdia0(3,3,iat)
      trc  = corrpre(1,1,iat)+corrpre(2,2,iat)+corrpre(3,3,iat)
      do t = 1, 3
        do s = 1, 3
          ! a11part (trace-corrected) + a01gp (raw, per the standard diamagnetic decomposition)
          sig_dia(t,s,iat) = ( SG*0.5d0*gdia0(s,t,iat) + SC*corrpre(t,s,iat) &
                 - merge(SG*0.5d0*trg0 + SC*trc, 0.0d0, t==s) &
                 + SA*a01(t,s,iat) ) * a2ppm
          sig_tot(t,s,iat) = sig_dia(t,s,iat) + sig_c(t,s,iat)
        end do
      end do
    end do

    ! --- Store isotropic shielding (ppm) to a tagarray for JSON output ---
    !     rows: dia, para_uncoupled, para_coupled, total_uncoupled, total_coupled
    !     stored atom-major (flat): atom a occupies entries 5*(a-1)+1 .. 5*a.
    call infos%dat%alloc_or_die(OQP_nmr_shielding, (/ 5*nat /), nmrout, description=OQP_nmr_shielding_comment)
    do iat = 1, nat
      nmrout(5*(iat-1)+1) = (sig_dia(1,1,iat)+sig_dia(2,2,iat)+sig_dia(3,3,iat))/3.0d0
      nmrout(5*(iat-1)+2) = (sig_u(1,1,iat)+sig_u(2,2,iat)+sig_u(3,3,iat))/3.0d0
      nmrout(5*(iat-1)+3) = (sig_c(1,1,iat)+sig_c(2,2,iat)+sig_c(3,3,iat))/3.0d0
      nmrout(5*(iat-1)+4) = nmrout(5*(iat-1)+1) + nmrout(5*(iat-1)+2)
      nmrout(5*(iat-1)+5) = nmrout(5*(iat-1)+1) + nmrout(5*(iat-1)+3)
    end do

    ! --- Emit parseable records ---
    inquire(unit=iw, opened=iw_open)
    if (.not. iw_open) open(unit=iw, file=infos%log_filename, position="append")
    write(iw,'(/,A)') 'GIAO_SHIELDING_DEBUG_BEGIN native-giao shielding (ppm)'
    write(iw,'(A,1X,I0)') 'GIAO_SHIELDING_DEBUG_NATOM', nat
    write(iw,'(A,1X,F10.6)') 'GIAO_SHIELDING_DEBUG_CX', scale_exch
    do iat = 1, nat
      do t = 1, 3
        do s = 1, 3
          write(iw,'(A,1X,I0,1X,I0,1X,I0,1X,ES24.16)') 'GIAO_SHIELDING_DEBUG_PARA_UNC', &
            iat, t, s, sig_u(t,s,iat)
          write(iw,'(A,1X,I0,1X,I0,1X,I0,1X,ES24.16)') 'GIAO_SHIELDING_DEBUG_PARA_CPL', &
            iat, t, s, sig_c(t,s,iat)
          write(iw,'(A,1X,I0,1X,I0,1X,I0,1X,ES24.16)') 'GIAO_SHIELDING_DEBUG_DIA', &
            iat, t, s, sig_dia(t,s,iat)
          write(iw,'(A,1X,I0,1X,I0,1X,I0,1X,ES24.16)') 'GIAO_SHIELDING_DEBUG_TOTAL', &
            iat, t, s, sig_tot(t,s,iat)
        end do
      end do
      write(iw,'(A,1X,I0,4(1X,ES24.16))') 'GIAO_SHIELDING_DEBUG_ISO', iat, &
        (sig_u(1,1,iat)+sig_u(2,2,iat)+sig_u(3,3,iat))/3.0d0, &
        (sig_c(1,1,iat)+sig_c(2,2,iat)+sig_c(3,3,iat))/3.0d0, &
        (sig_dia(1,1,iat)+sig_dia(2,2,iat)+sig_dia(3,3,iat))/3.0d0, &
        (sig_tot(1,1,iat)+sig_tot(2,2,iat)+sig_tot(3,3,iat))/3.0d0
    end do
    write(iw,'(A)') 'GIAO_SHIELDING_DEBUG_END'

    ! --- Human-readable production table (GIAO isotropic shielding) ---
    write(iw,'(2/)')
    write(iw,'(4x,a)') '========================================'
    write(iw,'(4x,a)') 'NMR nuclear magnetic shielding (GIAO)'
    write(iw,'(4x,a)') '========================================'
    write(iw,'(4x,a)') 'Gauge formulation: GIAO (London) atomic orbitals; gauge-origin independent.'
    if (open_shell) then
      write(iw,'(4x,a)') 'Reference: open shell (UHF/ROHF); diamagnetic from total density, '// &
        'paramagnetic summed over spin channels.'
    end if
    write(iw,'(/4x,a,f8.4,a)') 'Isotropic shielding (GIAO, ppm)   [exact-exchange c_x =', &
           scale_exch, ']'
    write(iw,'(4x,a)')  '   Atom    Z   sigma_dia   para_uncoupled   para_coupled'// &
           '   total_uncoupled   total_coupled'
    do iat = 1, nat
      write(iw,'(4x,i6,f6.1,5f16.6)') iat, basis%atoms%zn(iat), &
        (sig_dia(1,1,iat)+sig_dia(2,2,iat)+sig_dia(3,3,iat))/3.0d0, &
        (sig_u(1,1,iat)+sig_u(2,2,iat)+sig_u(3,3,iat))/3.0d0, &
        (sig_c(1,1,iat)+sig_c(2,2,iat)+sig_c(3,3,iat))/3.0d0, &
        (sig_dia(1,1,iat)+sig_dia(2,2,iat)+sig_dia(3,3,iat))/3.0d0 + &
        (sig_u(1,1,iat)+sig_u(2,2,iat)+sig_u(3,3,iat))/3.0d0, &
        (sig_tot(1,1,iat)+sig_tot(2,2,iat)+sig_tot(3,3,iat))/3.0d0
    end do
    close(iw)

    deallocate(gdia0, corrpre, a01, sig_dia, sig_tot)
    deallocate(h10p, s10p, twoe, twoe2, vj, vk, dm, dmp, h1ao, s1ao)
    if (allocated(vxa)) deallocate(vxa)
    if (allocated(vxb)) deallocate(vxb)
    deallocate(coords, zq, sig_u, sig_c)
    if (allocated(vkb))    deallocate(vkb)
    if (allocated(h1ao_b)) deallocate(h1ao_b)
    if (allocated(dm_b))   deallocate(dm_b)

  end subroutine nmr_giao_shielding_debug

!> Expand packed lower-triangular antisymmetric matrix to full form.
  subroutine expand_antisym(packed, full, n)
    real(kind=dp), intent(in)  :: packed(:)
    real(kind=dp), intent(out) :: full(:,:)
    integer, intent(in) :: n
    integer :: p, q, idx
    full = 0.0d0
    do p = 1, n
      do q = 1, p
        idx = q + p*(p-1)/2
        full(p,q) =  packed(idx)
        full(q,p) = -packed(idx)
      end do
    end do
  end subroutine expand_antisym

  subroutine unpack_sym(packed, full, n)
    real(kind=dp), intent(in) :: packed(:)
    real(kind=dp), intent(out) :: full(:,:)
    integer, intent(in) :: n
    integer :: i, j, ij
    full = 0.0d0
    do i = 1, n
      do j = 1, i
        ij = j + i*(i-1)/2
        full(i,j) = packed(ij)
        full(j,i) = packed(ij)
      end do
    end do
  end subroutine unpack_sym

!> MO transform keeping occupied ket columns: m(p,i) = sum_mn C(m,p) a(m,n) C(n,i),
!>  p = 1..nmo, i = 1..nocc.
  subroutine ao_to_mo_occ(a_ao, c, m_mo, nbf, nmo, nocc)
    real(kind=dp), intent(in)  :: a_ao(:,:), c(:,:)
    real(kind=dp), intent(out) :: m_mo(:,:)
    integer, intent(in) :: nbf, nmo, nocc
    real(kind=dp), allocatable :: tmp(:,:)
    allocate(tmp(nbf,nocc))
    tmp = matmul(a_ao, c(:,1:nocc))
    m_mo(1:nmo,1:nocc) = matmul(transpose(c(:,1:nmo)), tmp)
    deallocate(tmp)
  end subroutine ao_to_mo_occ

!> Uncoupled first-order equation (standard uncoupled first-order equation):
!>   hs = h1 - s1*e_i ;  mo1[vir,i] = -hs[vir,i]/(e_a-e_i) ;
!>   mo1[occ,i] = -0.5*s1[occ,i].
  subroutine solve_mo1_uncoupled(h1mo, s1mo, e, nocc, nmo, mo1)
    real(kind=dp), intent(in)  :: h1mo(:,:,:), s1mo(:,:,:), e(:)
    integer, intent(in) :: nocc, nmo
    real(kind=dp), intent(out) :: mo1(:,:,:)
    integer :: x, p, i
    real(kind=dp) :: hs
    mo1 = 0.0d0
    do x = 1, 3
      do i = 1, nocc
        do p = 1, nmo
          hs = h1mo(p,i,x) - s1mo(p,i,x)*e(i)
          if (p > nocc) then
            mo1(p,i,x) = -hs/(e(p)-e(i))
          else
            mo1(p,i,x) = -0.5d0*s1mo(p,i,x)
          end if
        end do
      end do
    end do
  end subroutine solve_mo1_uncoupled

!> Coupled GIAO CPHF/CPKS: fixed-point of
!>   mo1[vir,i] = -(hs[vir,i] + v1[vir,i])/(e_a-e_i),  mo1[occ,i] = -0.5 s1[occ,i],
!>  with v1 the exact-exchange response (scaled by c_x) of the imaginary
!>  antisymmetric first-order density built from the full mo1.  For c_x = 0 the
!>  loop is skipped and the result equals the uncoupled solution.
  subroutine solve_mo1_coupled(infos, basis, mo, h1mo, s1mo, e, nocc, nmo, &
                               scale_exch, mo1)
    use int2_compute, only: int2_compute_t
    use tdhf_lib, only: int2_td_data_t, mntoia
    use types, only: information
    use basis_tools, only: basis_set
    use messages, only: show_message
    real(kind=dp), intent(in) :: mo(:,:), h1mo(:,:,:), s1mo(:,:,:), e(:)
    integer, intent(in) :: nocc, nmo
    real(kind=dp), intent(in) :: scale_exch
    real(kind=dp), intent(out) :: mo1(:,:,:)
    type(information), target, intent(inout) :: infos
    type(basis_set), intent(in) :: basis

    integer, parameter :: maxit = 100
    real(kind=dp), parameter :: tol = 1.0d-9
    integer :: nbf, nvir, x, p, i, a, k, it
    real(kind=dp) :: diff, hs
    type(int2_compute_t) :: int2_driver
    type(int2_td_data_t), target :: kdat
    real(kind=dp), allocatable, target :: pa(:,:,:)
    real(kind=dp), allocatable :: gxv(:), mo1x(:,:), prev(:,:), gao(:,:)

    nbf = basis%nbf
    nvir = nmo - nocc

    ! Start from the uncoupled solution.
    call solve_mo1_uncoupled(h1mo, s1mo, e, nocc, nmo, mo1)
    if (abs(scale_exch) <= 1.0d-12) return

    allocate(pa(nbf,nbf,1), gxv(nocc*nvir), mo1x(nmo,nocc), &
             prev(nmo,nocc), gao(nbf,nbf), source=0.0d0)

    call int2_driver%init(basis, infos)
    ! NMR uses the native Rys ERI path only (the GIAO two-electron derivative
    ! and the reference data are Rys-based); never route through libint.
    int2_driver%rys_only = .true.
    call int2_driver%set_screening()
    kdat = int2_td_data_t(d2=pa, int_apb=.false., int_amb=.true., &
                          tamm_dancoff=.false., scale_exchange=scale_exch)

    do x = 1, 3
      mo1x = mo1(:,:,x)
      do it = 1, maxit
        prev = mo1x
        ! Imaginary antisymmetric AO first-order density from the full MO
        ! response (occ + vir rows); CGO-consistent normalization (no x2).
        call giao_pb_density(mo, mo1x, pa(:,:,1), nbf, nmo, nocc)
        call int2_driver%run(kdat)
        gao = 0.5d0*kdat%amb(:,:,1,1)
        ! Exchange response projected to the occ-vir block (i fast):
        ! gxv(i+(a-1)*nocc) = (C^T gao C)[i_occ, a_vir].
        call mntoia(gao, gxv, mo, mo, nocc, nocc)
        do i = 1, nocc
          do a = 1, nvir
            p = nocc + a
            k = i + (a-1)*nocc
            hs = h1mo(p,i,x) - s1mo(p,i,x)*e(i)
            ! mntoia returns the [occ,vir] block gxv(i,a); the response element
            ! needed here is the [vir,occ] entry v1(a,i) = -gxv(i,a) (the
            ! exchange image is antisymmetric).
            mo1x(p,i) = -(hs - gxv(k))/(e(p)-e(i))
          end do
          do p = 1, nocc
            mo1x(p,i) = -0.5d0*s1mo(p,i,x)
          end do
        end do
        diff = maxval(abs(mo1x - prev))
        if (diff < tol) exit
      end do
      if (diff >= tol) then
        call show_message('WARNING: GIAO coupled magnetic response (CPHF) did &
          &not converge within the iteration limit; shieldings may be inaccurate')
      end if
      mo1(:,:,x) = mo1x
    end do

    call int2_driver%clean()
    deallocate(pa, gxv, mo1x, prev, gao)
  end subroutine solve_mo1_coupled

!> Imaginary antisymmetric AO first-order density from a full MO response vector
!>  (CGO-consistent normalization, no double-occupancy factor):
!>   D = C mo1 orbo^T ;  pa = D - D^T.
  subroutine giao_pb_density(mo, mo1x, pa, nbf, nmo, nocc)
    real(kind=dp), intent(in) :: mo(:,:), mo1x(:,:)
    real(kind=dp), intent(out) :: pa(:,:)
    integer, intent(in) :: nbf, nmo, nocc
    real(kind=dp), allocatable :: dleft(:,:)
    allocate(dleft(nbf,nbf))
    dleft = matmul(mo(:,1:nmo), matmul(mo1x, transpose(mo(:,1:nocc))))
    pa = dleft - transpose(dleft)
    deallocate(dleft)
  end subroutine giao_pb_density

!> Paramagnetic shielding tensor for one nucleus (standard paramagnetic convention):
!>   dm10(a,b,x) = 2 sum_{p,i} C(a,p) mo1(p,i,x) C(b,i) ;
!>   sigma_para[x,y] = 2 sum_{a,b} dm10(a,b,x) * h01i(b,a,y).
  !> Paramagnetic shielding contribution from one spin channel: MO transform,
  !> CPHF (uncoupled + coupled), and PSO contraction, ACCUMULATED into sig_u/sig_c.
  !> occ_factor = 2 for RHF (closed shell), 1 for each UHF spin channel.
  subroutine giao_para_channel(infos, basis, mo, e, nocc, nmo, nbf, nat, coords, &
                               h1ao, s1ao, scale_exch, occ_factor, sig_u, sig_c)
    use types, only: information
    use basis_tools, only: basis_set
    use int1, only: pso_integrals
    type(information), target, intent(inout) :: infos
    type(basis_set), intent(in) :: basis
    real(kind=dp), intent(in) :: mo(:,:), e(:), coords(:,:)
    real(kind=dp), intent(in) :: h1ao(:,:,:), s1ao(:,:,:), scale_exch, occ_factor
    integer, intent(in) :: nocc, nmo, nbf, nat
    real(kind=dp), intent(inout) :: sig_u(:,:,:), sig_c(:,:,:)
    real(kind=dp), allocatable :: h1mo(:,:,:), s1mo(:,:,:), mo1u(:,:,:), mo1c(:,:,:)
    real(kind=dp), allocatable :: pso(:,:,:), st(:,:)
    integer :: c, iat

    allocate(h1mo(nmo,nocc,3), s1mo(nmo,nocc,3), mo1u(nmo,nocc,3), mo1c(nmo,nocc,3), &
             pso(nbf,nbf,3), st(3,3), source=0.0d0)
    do c = 1, 3
      call ao_to_mo_occ(h1ao(:,:,c), mo, h1mo(:,:,c), nbf, nmo, nocc)
      call ao_to_mo_occ(s1ao(:,:,c), mo, s1mo(:,:,c), nbf, nmo, nocc)
    end do
    call solve_mo1_uncoupled(h1mo, s1mo, e, nocc, nmo, mo1u)
    call solve_mo1_coupled(infos, basis, mo, h1mo, s1mo, e, nocc, nmo, scale_exch, mo1c)
    do iat = 1, nat
      call pso_integrals(basis, coords(:,iat), pso)
      call para_tensor(mo1u, mo, pso, nbf, nmo, nocc, st, occ_factor)
      sig_u(:,:,iat) = sig_u(:,:,iat) + st
      call para_tensor(mo1c, mo, pso, nbf, nmo, nocc, st, occ_factor)
      sig_c(:,:,iat) = sig_c(:,:,iat) + st
    end do
    deallocate(h1mo, s1mo, mo1u, mo1c, pso, st)
  end subroutine giao_para_channel

  !> Semicanonicalize ROHF orbitals for spin sigma: diagonalize the occ-occ and
  !> vir-vir blocks of C^T F_sigma C, returning semicanonical orbitals (csc) and
  !> orbital energies (esc).  ROHF stores a single orbital set with effective-Fock
  !> eigenvalues; the UHF-like CPHF needs proper spin orbital energies (the
  !> ground-state spin Fock matrices F_a/F_b from OQP_FOCK_A/B).
  subroutine semicanon_orbitals(fock_p, c0, nbf, nmo, nocc, csc, esc)
    use eigen, only: diag_symm_full
    real(kind=dp), intent(in) :: fock_p(:), c0(:,:)
    integer, intent(in) :: nbf, nmo, nocc
    real(kind=dp), intent(out) :: csc(:,:), esc(:)
    real(kind=dp), allocatable :: fao(:,:), fmo(:,:), blk(:,:), eb(:)
    integer :: nvir, ierr
    nvir = nmo - nocc
    allocate(fao(nbf,nbf), fmo(nmo,nmo))
    call unpack_sym(fock_p, fao, nbf)
    fmo = matmul(transpose(c0(:,1:nmo)), matmul(fao, c0(:,1:nmo)))
    csc = c0(:,1:nmo); esc = 0.0d0
    if (nocc > 0) then
      allocate(blk(nocc,nocc), eb(nocc))
      blk = fmo(1:nocc,1:nocc)
      call diag_symm_full(1, nocc, blk, nocc, eb, ierr)
      csc(:,1:nocc) = matmul(c0(:,1:nocc), blk)
      esc(1:nocc) = eb(1:nocc)
      deallocate(blk, eb)
    end if
    if (nvir > 0) then
      allocate(blk(nvir,nvir), eb(nvir))
      blk = fmo(nocc+1:nmo, nocc+1:nmo)
      call diag_symm_full(1, nvir, blk, nvir, eb, ierr)
      csc(:,nocc+1:nmo) = matmul(c0(:,nocc+1:nmo), blk)
      esc(nocc+1:nmo) = eb(1:nvir)
      deallocate(blk, eb)
    end if
    deallocate(fao, fmo)
  end subroutine semicanon_orbitals

  subroutine para_tensor(mo1, mo, h01i, nbf, nmo, nocc, sig, occ_factor)
    real(kind=dp), intent(in) :: mo1(:,:,:), mo(:,:), h01i(:,:,:)
    integer, intent(in) :: nbf, nmo, nocc
    real(kind=dp), intent(out) :: sig(:,:)
    real(kind=dp), intent(in), optional :: occ_factor
    integer :: x, y, a, b
    real(kind=dp), allocatable :: dm10(:,:,:)
    real(kind=dp) :: acc, ofac
    ofac = 2.0d0                 ! RHF closed-shell occupation; UHF per spin = 1
    if (present(occ_factor)) ofac = occ_factor
    allocate(dm10(nbf,nbf,3))
    do x = 1, 3
      dm10(:,:,x) = ofac*matmul(mo(:,1:nmo), matmul(mo1(:,:,x), transpose(mo(:,1:nocc))))
    end do
    do x = 1, 3
      do y = 1, 3
        acc = 0.0d0
        do b = 1, nbf
          do a = 1, nbf
            acc = acc + dm10(a,b,x)*h01i(b,a,y)
          end do
        end do
        ! OpenQP pso_integrals stores the negative of libcint int1e_prinvxp
        ! (h01i); the leading 2 is the +c.c. factor of the the standard paramagnetic routine.
        sig(x,y) = -2.0d0*acc
      end do
    end do
    deallocate(dm10)
  end subroutine para_tensor

end module nmr_giao_shielding_mod
