compute_vibrational_frequencies Subroutine

   subroutine compute_vibrational_frequencies(hessian, element_numbers, frequencies, &
                                              eigenvalues_out, eigenvectors, &
                                              coordinates, project_trans_rot, &
                                              projected_hessian_out)
      !! Compute vibrational frequencies from the Hessian matrix.
      !!
      !! Algorithm:
      !! 1. Mass-weight the Hessian: H_mw = M^{-1/2} * H * M^{-1/2}
      !! 2. Optionally project out translation/rotation modes
      !! 3. Diagonalize H_mw to get eigenvalues
      !! 4. Convert eigenvalues to frequencies in cm⁻¹
      !!
      !! Negative eigenvalues produce negative frequencies (imaginary modes,
      !! indicating transition states or saddle points).
      real(dp), intent(in) :: hessian(:, :)
         !! Hessian matrix in Hartree/Bohr² (3*N x 3*N)
      integer, intent(in) :: element_numbers(:)
         !! Atomic numbers for each atom (N atoms)
      real(dp), allocatable, intent(out) :: frequencies(:)
         !! Vibrational frequencies in cm⁻¹ (3*N modes, or 3*N-6 if projected)
      real(dp), allocatable, intent(out), optional :: eigenvalues_out(:)
         !! Raw eigenvalues from diagonalization (Hartree/Bohr²/amu)
      real(dp), allocatable, intent(out), optional :: eigenvectors(:, :)
         !! Normal mode eigenvectors (3*N x 3*N), columns are modes
      real(dp), intent(in), optional :: coordinates(:, :)
         !! Atomic coordinates in Bohr (3, N) - required for projection
      logical, intent(in), optional :: project_trans_rot
         !! If true, project out translation/rotation modes (requires coordinates)
      real(dp), allocatable, intent(out), optional :: projected_hessian_out(:, :)
         !! Mass-weighted Hessian after trans/rot projection (before diagonalization)

      real(dp), allocatable :: mw_hessian(:, :)
      real(dp), allocatable :: eigenvalues(:)
      integer :: n_coords, info, i
      logical :: do_projection

      n_coords = size(hessian, 1)

      ! Check if projection is requested
      do_projection = .false.
      if (present(project_trans_rot)) then
         if (project_trans_rot) then
            if (.not. present(coordinates)) then
               ! Cannot project without coordinates - fall back to no projection
               call logger%warning("Missing coordinates, not projecting out tran/rot motions")
               do_projection = .false.
            else
               do_projection = .true.
            end if
         end if
      end if

      ! Mass-weight the Hessian
      call mass_weight_hessian(hessian, element_numbers, mw_hessian)

      ! Optionally project out translation/rotation modes
      if (do_projection) then
         call project_translation_rotation(mw_hessian, coordinates, element_numbers)
      end if

      ! Return projected Hessian if requested (before diagonalization destroys it)
      if (present(projected_hessian_out)) then
         allocate (projected_hessian_out(n_coords, n_coords))
         projected_hessian_out = mw_hessian
      end if

      ! Allocate eigenvalue storage
      allocate (eigenvalues(n_coords))

      ! Diagonalize the mass-weighted Hessian
      ! pic_syev overwrites mw_hessian with eigenvectors (if jobz='V', default)
      call pic_syev(mw_hessian, eigenvalues, info=info)

      if (info /= 0) then
         ! Eigenvalue decomposition failed
         call logger%error("Eigenvalue decomposition in vibrational frequencies failed")
         allocate (frequencies(n_coords))
         frequencies = 0.0_dp
         return
      end if

      ! Convert eigenvalues to frequencies in cm⁻¹
      allocate (frequencies(n_coords))
      do i = 1, n_coords
         if (eigenvalues(i) >= 0.0_dp) then
            ! Real frequency
            frequencies(i) = sqrt(eigenvalues(i)*AU_TO_CM1)
         else
            ! Imaginary frequency (negative eigenvalue) - report as negative
            frequencies(i) = -sqrt(abs(eigenvalues(i))*AU_TO_CM1)
         end if
      end do

      ! Return eigenvalues if requested
      if (present(eigenvalues_out)) then
         allocate (eigenvalues_out(n_coords))
         eigenvalues_out = eigenvalues
      end if

      ! Return eigenvectors if requested
      if (present(eigenvectors)) then
         allocate (eigenvectors(n_coords, n_coords))
         eigenvectors = mw_hessian
      end if

      deallocate (eigenvalues, mw_hessian)

   end subroutine compute_vibrational_frequencies