xtb_calc_hessian Subroutine

   subroutine xtb_calc_hessian(this, fragment, result)
      !! Calculate Hessian using finite differences of gradients
      !!
      !! Since tblite does not natively support analytic Hessians, this routine
      !! computes the Hessian numerically via central finite differences:
      !!   H[i,j] = (grad_j(x_i + h) - grad_j(x_i - h)) / (2h)
      !!
      !! This requires 6N gradient calculations (forward and backward for each coordinate)
      use mqc_finite_differences, only: generate_perturbed_geometries, displaced_geometry_t, &
                                        finite_diff_hessian_from_gradients, finite_diff_dipole_derivatives, &
                                        DEFAULT_DISPLACEMENT
      use pic_logger, only: logger => global_logger
      use pic_io, only: to_char

      class(xtb_method_t), intent(in) :: this
      type(physical_fragment_t), intent(in) :: fragment
      type(calculation_result_t), intent(out) :: result

      type(displaced_geometry_t), allocatable :: forward_geoms(:), backward_geoms(:)
      real(dp), allocatable :: forward_gradients(:, :, :)  ! (n_displacements, 3, n_atoms)
      real(dp), allocatable :: backward_gradients(:, :, :)  ! (n_displacements, 3, n_atoms)
      real(dp), allocatable :: forward_dipoles(:, :)   ! (n_displacements, 3) for IR intensities
      real(dp), allocatable :: backward_dipoles(:, :)  ! (n_displacements, 3) for IR intensities
      type(calculation_result_t) :: disp_result
      real(dp) :: displacement
      integer :: n_atoms, n_displacements, i
      logical :: compute_dipole_derivs

      n_atoms = fragment%n_atoms
      n_displacements = 3*n_atoms
      displacement = DEFAULT_DISPLACEMENT

      if (this%verbose) then
         call logger%info("XTB: Computing Hessian via finite differences")
         call logger%info("  Method: Central differences of gradients")
         call logger%info("  Atoms: "//to_char(n_atoms))
         call logger%info("  Gradient calculations needed: "//to_char(2*n_displacements))
         call logger%info("  Finite difference step size: "//to_char(displacement)//" Bohr")
      end if

      ! Generate all perturbed geometries
      call generate_perturbed_geometries(fragment, displacement, forward_geoms, backward_geoms)

      ! Allocate storage for gradients at displaced geometries
      allocate (forward_gradients(n_displacements, 3, n_atoms))
      allocate (backward_gradients(n_displacements, 3, n_atoms))

      ! Allocate storage for dipoles at displaced geometries (for IR intensities)
      allocate (forward_dipoles(n_displacements, 3))
      allocate (backward_dipoles(n_displacements, 3))
      forward_dipoles = 0.0_dp
      backward_dipoles = 0.0_dp
      compute_dipole_derivs = .true.  ! Will be set to false if any dipole is missing

      ! Compute gradients at all forward-displaced geometries
      if (this%verbose) then
         call logger%info("  Computing forward-displaced gradients...")
      end if
      do i = 1, n_displacements

         ! Forward
         call this%calc_gradient(forward_geoms(i)%geometry, disp_result)
         if (disp_result%has_error .or. .not. disp_result%has_gradient) then
            call result%error%set(ERROR_GENERIC, "Failed to compute gradient for forward displacement "//to_char(i))
            result%has_error = .true.
            call disp_result%destroy()
            return
         end if
         forward_gradients(i, :, :) = disp_result%gradient
         ! Capture dipole for IR intensity calculation
         if (disp_result%has_dipole) then
            forward_dipoles(i, :) = disp_result%dipole
         else
            compute_dipole_derivs = .false.
         end if
         call disp_result%destroy()

         ! Backward
         call this%calc_gradient(backward_geoms(i)%geometry, disp_result)
         if (disp_result%has_error .or. .not. disp_result%has_gradient) then
            call result%error%set(ERROR_GENERIC, "Failed to compute gradient for backward displacement "//to_char(i))
            result%has_error = .true.
            call disp_result%destroy()
            return
         end if
         backward_gradients(i, :, :) = disp_result%gradient
         ! Capture dipole for IR intensity calculation
         if (disp_result%has_dipole) then
            backward_dipoles(i, :) = disp_result%dipole
         else
            compute_dipole_derivs = .false.
         end if
         call disp_result%destroy()

      end do
      if (this%verbose) then
         call logger%info("  Forward and backward gradient calculations complete ")
      end if
      ! Compute Hessian from finite differences
      if (this%verbose) then
         call logger%info("  Assembling Hessian matrix...")
      end if
      call finite_diff_hessian_from_gradients(fragment, forward_gradients, backward_gradients, &
                                              displacement, result%hessian)

      ! Compute dipole derivatives for IR intensity calculation
      if (compute_dipole_derivs) then
         call finite_diff_dipole_derivatives(n_atoms, forward_dipoles, backward_dipoles, &
                                             displacement, result%dipole_derivatives)
         result%has_dipole_derivatives = .true.
         if (this%verbose) then
            call logger%info("  Dipole derivatives computed for IR intensities")
         end if
      end if

      ! Also compute energy and gradient at reference geometry for completeness
      call this%calc_gradient(fragment, disp_result)
      if (disp_result%has_error) then
         result%error = disp_result%error
         result%has_error = .true.
         call disp_result%destroy()
         return
      end if
      result%energy = disp_result%energy
      result%has_energy = disp_result%has_energy
      if (disp_result%has_gradient) then
         allocate (result%gradient(3, n_atoms))
         result%gradient = disp_result%gradient
         result%has_gradient = .true.
      end if
      call disp_result%destroy()

      result%has_hessian = .true.

      if (this%verbose) then
         call logger%info("  Hessian calculation complete")
      end if

      ! Cleanup
      deallocate (forward_gradients, backward_gradients)
      deallocate (forward_dipoles, backward_dipoles)
      do i = 1, n_displacements
         call forward_geoms(i)%destroy()
         call backward_geoms(i)%destroy()
      end do
      deallocate (forward_geoms, backward_geoms)

   end subroutine xtb_calc_hessian