!> @brief Thin interface between OpenQP and OpenTrustRegion.
!> @detail Provides callback glue so OpenTrustRegion’s generic trust-region
!>         optimizer can drive OpenQP’s TRAH SCF updates without requiring
!>         an external solver class. Exposes:
!>           - init_trah_solver: bind OpenQP state to module pointers
!>           - run_trah_solver : configure and invoke the OTR solver
!>           - update_orbs     : objective/gradient/Hdiag callback
!>           - hess_x_cb       : Hessian–vector product callback
!>           - obj_func        : energy-only evaluation (trial move)
!>           - logger          : forwards OTR log lines to OpenQP I/O
!> @author Mohsen Mazaherifar
!> @date August 2025
module otr_interface

  use, intrinsic :: iso_c_binding, only: c_bool
  use opentrustregion, only: solver, update_orbs_type,&
          obj_func_type, hess_x_type, logger_type, rp, ip, &
          solver_settings_type, stability_settings_type, &
          default_solver_settings
  use mathlib, only: unpack_matrix
  use scf_converger, only: trah_converger, scf_conv_trah_result, scf_conv_result
  use scf_addons,only: compute_energy,calc_fock
  use precision, only: dp
  use types, only:information
  use mod_dft_molgrid, only: dft_grid_t
  use basis_tools,      only: basis_set
  use guess, only:get_ab_initio_density
  use scf_addons, only: scf_energy_t

  implicit none

  ! Module-level state for callbacks
  class(information), pointer :: infos       ! OpenQP information object
  type(dft_grid_t), pointer :: molgrid
  type(trah_converger), pointer :: conv
  type(scf_energy_t), pointer :: energy
  integer :: iter_otr
  real(dp) :: grad_norm

  real(dp), allocatable :: work1(:,:), work2(:,:)


contains
  !> @brief Initialize the OTR–OpenQP bridge and working buffers.
  !> @detail Stores references to OpenQP objects (infos, molgrid, TRAH converger,
  !>         energy accumulator), allocates temporary work arrays, and zeros the
  !>         incremental Fock/Density buffers used for ΔD updates.
  !> @param[inout] infos_in   OpenQP information/control object (target).
  !> @param[in]    molgrid_in DFT molecular grid (target).
  !> @param[inout] conv_in    TRAH converger (provides MO/D/Fock buffers).
  !> @param[inout] energy_in  SCF energy structure to be updated.
  !> @author Mohsen Mazaherifar
  !> @date August 2025
  subroutine init_trah_solver(infos_in, molgrid_in, conv_in, energy_in)
    class(information), intent(inout), target :: infos_in
    type(dft_grid_t), intent(in), target :: molgrid_in
    class(trah_converger), intent(inout), target :: conv_in
    class(scf_energy_t), intent(inout), target :: energy_in
    type(basis_set), pointer :: basis
    ! Initialize module state
    infos   => infos_in
    molgrid => molgrid_in
    conv => conv_in
    energy => energy_in
    iter_otr = 0

    basis => infos%basis
    allocate(work1(conv%nbf,conv%nbf), work2(conv%nbf,conv%nbf))

    conv%f_old = 0.0_dp
    conv%d_old = 0.0_dp

  end subroutine init_trah_solver

  !> @brief Configure and run the OpenTrustRegion driver.
  !> @detail Wires the required callbacks (`update_orbs`, `obj_func`, `logger`),
  !>         maps OpenQP control flags to OTR options (stability, line-search,
  !>         Davidson/Jacobi–Davidson, trust-radius settings), executes the solve,
  !>         and returns iteration/error status in `res`.
  !> @param[inout] res  Output SCF converger result (TRAH-specific fields filled).
  !> @note Updates the active OpenQP buffers (MO/Fock) upon return.
  !> @author Mohsen Mazaherifar
  !> @date August 2025
  subroutine run_trah_solver(res)
    procedure(update_orbs_type), pointer :: p_update
    procedure(obj_func_type),   pointer :: p_obj
    procedure(logger_type),     pointer :: p_log
    class(scf_conv_result), intent(inout) :: res
    type(solver_settings_type) :: settings
    logical(kind=4) :: stability, line_search, davidson,&
                         jacobi_davidson, prefer_jacobi_davidson
    integer(ip) :: error, n_random_trial_vectors, n_micro,&
                        n_param, max_iter, verbose
    real(dp) :: start_trust_radius, global_red_factor,&
                          local_red_factor, conv_tol

    n_param = conv%n_param
    max_iter = int(infos%control%maxit, kind=ip)
    conv_tol = real(infos%control%conv, kind=rp)
    verbose  = int(3, kind=ip)
    settings = default_solver_settings
    settings%conv_tol = conv_tol
    settings%verbose = verbose
    settings%stability = (infos%control%trh_stab .eqv. .true._c_bool)
    settings%line_search = (infos%control%trh_ls   .eqv. .true._c_bool)
    select case (infos%control%trh_sub_solver)
    case (0)
        settings%subsystem_solver = "davidson"
    case (1)
        settings%subsystem_solver = "jacobi_davidson"
    case (2)
        settings%subsystem_solver = "tcg"
    case default
        error stop "Invalid trh_sub_solver value"
    end select
    settings%n_random_trial_vectors = int(infos%control%trh_nrtv, kind=ip)
    settings%start_trust_radius = real(infos%control%trh_r0, kind=ip)
    settings%jacobi_davidson_start = int(infos%control%trh_jd_start, kind=ip)
    settings%global_red_factor = real(infos%control%trh_gred, kind=rp)
    settings%local_red_factor = real(infos%control%trh_lred, kind=rp)
    settings%n_macro = max_iter
    settings%n_micro = int(infos%control%trh_nmic, kind=ip)

    settings%logger => logger

    call print_trah_settings(settings)
    ! Bind callbacks
    p_update => update_orbs
    p_obj    => obj_func

    call solver(p_update, p_obj, n_param, error, settings)

    conv%dat%buffer(conv%dat%slot)%mo_a = conv%mo_a
    conv%dat%buffer(conv%dat%slot)%focks = conv%fock_ao
    if (infos%control%scftype>1) then
      conv%dat%buffer(conv%dat%slot)%mo_b = conv%mo_b
    end if

    select type (res)
    class is (scf_conv_trah_result)
      res%iter = iter_otr
    end select

    if (error /= 0) then
      write(*,*) 'OpenTrustRegion solver failed.'
      res%ierr = 4
      select type (res)
      class is (scf_conv_trah_result)
        res%iter = max_iter
      end select
    else
      if(grad_norm>conv_tol) then
        write(*,*) 'Trust radius too small. Convergence criterion&
        is not fulfilled but calculation should be converged up to floating&
        point precision.'
        res%error = min(conv_tol*0.99,grad_norm)
      else
        res%error = grad_norm
      end if
    endif


    if (allocated(work1)) deallocate(work1)
    if (allocated(work2)) deallocate(work2)

  end subroutine run_trah_solver

  !> @brief Objective/gradient/Hessian-diagonal callback used by OTR.
  !> @detail Applies orbital rotations `kappa` to (α[,β]) MOs, rebuilds densities,
  !>         constructs Fock via `calc_fock` (using incremental ΔD/ΔF when available),
  !>         then forms orbital-rotation gradient and Hessian diagonal with
  !>         `conv%calc_g_h`. Also binds the Hessian–vector product callback.
  !> @param[in]   kappa    Packed rotation vector(s).
  !> @param[out]  func     Objective value (total electronic energy).
  !> @param[out]  grad     Objective gradient in rotation coordinates.
  !> @param[out]  h_diag   Diagonal of approximate Hessian in rotation space.
  !> @param[out]  hess_x_funptr Pointer to Hessian–vector product routine.
  !> @author Mohsen Mazaherifar
  !> @date August 2025
  subroutine update_orbs(kappa, func, grad, h_diag, hess_x_funptr, error)
    real(dp), intent(in), target :: kappa(:)
    real(dp), intent(out) :: func
    real(dp), intent(out), target :: grad(:), h_diag(:)
    procedure(hess_x_type), intent(out), pointer :: hess_x_funptr
    integer(ip), intent(out) :: error
    type(basis_set), pointer :: basis
    integer :: nschwz

    basis => infos%basis
    iter_otr = iter_otr + 1
    ! Rotate orbitals
    select case (infos%control%scftype)
    case (1)
      call conv%rotate_orbs(kappa, conv%nbf, conv%nocc_a, conv%mo_a)
      call get_ab_initio_density(conv%dens(:,1), conv%mo_a, conv%dens(:,1), conv%mo_a,infos,basis)
      call calc_fock(basis, infos, molgrid, conv%fock_ao, energy, conv%mo_a, conv%dens, conv%mo_b, nschwz, conv%f_old, conv%d_old)
      call conv%calc_g_h(grad, h_diag)
    case (2)
      call conv%rotate_orbs(kappa(1:conv%nocc_a*conv%nvir_a), conv%nbf, conv%nocc_a, conv%mo_a)
      call conv%rotate_orbs(kappa(conv%nocc_a*conv%nvir_a+1:), conv%nbf, conv%nocc_b, conv%mo_b)
      call get_ab_initio_density(conv%dens(:,1), conv%mo_a, conv%dens(:,2), conv%mo_b,infos,basis)
      call calc_fock(basis, infos, molgrid, conv%fock_ao, energy, conv%mo_a, conv%dens, conv%mo_b, nschwz, conv%f_old, conv%d_old)
      call conv%calc_g_h(grad, h_diag)
    case (3)
      call conv%rotate_orbs(kappa, conv%nbf, conv%nocc_a, conv%mo_a)
      conv%mo_b = conv%mo_a
      call get_ab_initio_density(conv%dens(:,1), conv%mo_a, conv%dens(:,2), conv%mo_b,infos,basis)
      call calc_fock(basis, infos, molgrid, conv%fock_ao, energy, conv%mo_a, conv%dens, conv%mo_b, nschwz, conv%f_old, conv%d_old)
      call conv%calc_g_h(grad, h_diag)
    end select
    grad_norm = sqrt(dot_product(grad, grad)/conv%n_param)
    func = compute_energy(energy)
    conv%etot = func
    hess_x_funptr => hess_x_cb
    h_diag = 2.0_dp * h_diag
    grad = 2.0_dp * grad
  end subroutine update_orbs

  !> @brief hess_x_cb.
  !> @author Mohsen Mazaherifar
  !> @date August 2025
  subroutine hess_x_cb(x, hx, error)
    real(dp), intent(in), target :: x(:)
    real(dp), intent(out), target :: hx(:)
    integer(ip), intent(out) :: error
    call conv%calc_h_op(infos, x, hx)
    hx = 2.0_dp * hx
  end subroutine hess_x_cb

  !> @brief Energy-only objective for a trial move (no gradient).
  !> @detail Rotates temporary copies of the MOs according to `kappa`, rebuilds
  !>         densities, recomputes Fock and energies, and returns the total energy.
  !>         Used by line-search/auxiliary steps in OTR.
  !> @param[in]  kappa  Packed rotation vector(s).
  !> @return     val    Total electronic energy at the trial point.
  !> @author Mohsen Mazaherifar
  !> @date August 2025
  function obj_func(kappa, error) result(val)
    real(dp), intent(in), target :: kappa(:)
    integer(ip), intent(out) :: error
    real(dp)             :: val
    type(basis_set), pointer :: basis
    integer :: nschwz
    basis => infos%basis
    select case(infos%control%scftype)
    case (1)
      work1 = conv%mo_a
      call conv%rotate_orbs(kappa, conv%nbf, conv%nocc_a, work1)
      call get_ab_initio_density(conv%dens(:,1), work1, conv%dens(:,1), work1,infos,basis)
      call calc_fock(basis, infos, molgrid, conv%fock_ao, energy, work1, conv%dens, work1, nschwz, conv%f_old, conv%d_old)
    case (2)
      work1 = conv%mo_a
      work2 = conv%mo_b
      call conv%rotate_orbs(kappa(1:conv%nvir_a*conv%nocc_a), conv%nbf, conv%nocc_a, work1)
      call conv%rotate_orbs(kappa(conv%nvir_a*conv%nocc_a+1:), conv%nbf, conv%nocc_b, work2)
      call get_ab_initio_density(conv%dens(:,1), work1, conv%dens(:,2), work2,infos,basis)
      call calc_fock(basis, infos, molgrid, conv%fock_ao, energy, work1, conv%dens, work2, nschwz, conv%f_old, conv%d_old)
    case (3)
      work1 = conv%mo_a
      work2 = conv%mo_b
      call conv%rotate_orbs(kappa, conv%nbf, conv%nocc_a, work1)
      work2 = work1
      call get_ab_initio_density(conv%dens(:,1), work1, conv%dens(:,2), work2,infos,basis)
      call calc_fock(basis, infos, molgrid, conv%fock_ao, energy, work1, conv%dens, work2, nschwz, conv%f_old, conv%d_old)
    end select
    val = compute_energy(energy)
  end function obj_func

  subroutine logger(message)
    use io_constants, only: IW
    implicit none
    character(*), intent(in) :: message
    write(IW, "(A)") trim(message)
  end subroutine
  subroutine print_trah_settings(settings)
      use io_constants, only : IW
      implicit none

      type(solver_settings_type), intent(in) :: settings

      write(IW, '(5X, a)') "----------------------------------------"
      write(IW, '(6X, a)') "TRAH / Trust-Region Augmented Hessian Settings"
      write(IW, '(5X, a)') "----------------------------------------"
      write(IW, '(7X, a, es12.5)') "conv_tol                : ", settings%conv_tol
      write(IW, '(7X, a, i0)')     "verbose                 : ", settings%verbose
      write(IW, '(7X, a, l1)')     "stability               : ", settings%stability
      write(IW, '(7X, a, l1)')     "line_search             : ", settings%line_search
      write(IW, '(7X, a, a)')      "subsystem_solver        : ", &
                                    trim(settings%subsystem_solver)
      write(IW, '(7X, a, i0)')     "n_random_trial_vectors  : ", &
                                    settings%n_random_trial_vectors
      write(IW, '(7X, a, es12.5)') "start_trust_radius      : ", &
                                    settings%start_trust_radius
      write(IW, '(7X, a, i0)')     "jacobi_davidson_start   : ", &
                                    settings%jacobi_davidson_start
      write(IW, '(7X, a, es12.5)') "global_red_factor       : ", &
                                    settings%global_red_factor
      write(IW, '(7X, a, es12.5)') "local_red_factor        : ", &
                                    settings%local_red_factor
      write(IW, '(7X, a, i0)')     "n_macro                 : ", settings%n_macro
      write(IW, '(7X, a, i0)')     "n_micro                 : ", settings%n_micro
      write(IW, '(5X, a)') "----------------------------------------"
      flush(IW)
  end subroutine print_trah_settings
end module otr_interface
