finite_diff_hessian_from_gradients Subroutine

   subroutine finite_diff_hessian_from_gradients(reference_geom, forward_gradients, backward_gradients, &
                                                 displacement, hessian)
      !! Compute Hessian matrix from finite differences of gradients
      !!
      !! Uses central finite differences: H_ij = (grad_i(+h) - grad_i(-h)) / (2h)
      !!
      !! Args:
      !!   reference_geom: Reference geometry (for dimensioning)
      !!   forward_gradients: Gradients at forward-displaced geometries (3*n_atoms, 3, n_atoms)
      !!   backward_gradients: Gradients at backward-displaced geometries (3*n_atoms, 3, n_atoms)
      !!   displacement: Step size used in Bohr
      !!   hessian: Output Hessian matrix (3*n_atoms, 3*n_atoms)
      type(physical_fragment_t), intent(in) :: reference_geom
      real(dp), intent(in) :: forward_gradients(:, :, :)   !! (n_displacements, 3, n_atoms)
      real(dp), intent(in) :: backward_gradients(:, :, :)  !! (n_displacements, 3, n_atoms)
      real(dp), intent(in) :: displacement
      real(dp), intent(out), allocatable :: hessian(:, :)  !! (3*n_atoms, 3*n_atoms)

      integer :: n_atoms, n_coords
      integer :: iatom, jatom, icoord, jcoord
      integer :: i_global, j_global
      integer :: disp_idx

      n_atoms = reference_geom%n_atoms
      n_coords = 3*n_atoms

      allocate (hessian(n_coords, n_coords))
      hessian = 0.0_dp

      ! Build Hessian using central differences
      ! H[i,j] = d²E/(dx_i dx_j) = (dE/dx_j at x_i+h - dE/dx_j at x_i-h) / (2h)

      disp_idx = 0
      do iatom = 1, n_atoms
         do icoord = 1, 3
            disp_idx = disp_idx + 1
            i_global = 3*(iatom - 1) + icoord

            ! For each displacement, compute derivatives of all gradient components
            do jatom = 1, n_atoms
               do jcoord = 1, 3
                  j_global = 3*(jatom - 1) + jcoord

                  ! Central difference: (grad_j(+h) - grad_j(-h)) / (2h)
                  hessian(i_global, j_global) = &
                     (forward_gradients(disp_idx, jcoord, jatom) - &
                      backward_gradients(disp_idx, jcoord, jatom))/(2.0_dp*displacement)
               end do
            end do
         end do
      end do

      ! Symmetrize the Hessian: H = (H + H^T) / 2
      ! This reduces numerical noise from finite differences
      do i_global = 1, n_coords
         do j_global = i_global + 1, n_coords
            hessian(i_global, j_global) = 0.5_dp*(hessian(i_global, j_global) + hessian(j_global, i_global))
            hessian(j_global, i_global) = hessian(i_global, j_global)
         end do
      end do

   end subroutine finite_diff_hessian_from_gradients