print_vibrational_analysis Subroutine

   subroutine print_vibrational_analysis(frequencies, reduced_masses, force_constants, &
                                         cartesian_displacements, element_numbers, &
                                         force_constants_mdyne, print_displacements, &
                                         n_atoms, ir_intensities, &
                                         coordinates, electronic_energy, &
                                         temperature, pressure)
      !! Print vibrational analysis results in a formatted table.
      !!
      !! Output format is similar to Gaussian, with frequencies grouped in columns.
      !! Optionally prints Cartesian displacement vectors for each mode.
      !! If coordinates and electronic_energy are provided, also computes and prints
      !! thermochemistry using the RRHO approximation.
      real(dp), intent(in) :: frequencies(:)
         !! Vibrational frequencies in cm⁻¹
      real(dp), intent(in) :: reduced_masses(:)
         !! Reduced masses in amu
      real(dp), intent(in) :: force_constants(:)
         !! Force constants in Hartree/Bohr² (or mdyne/Å if force_constants_mdyne provided)
      real(dp), intent(in) :: cartesian_displacements(:, :)
         !! Cartesian displacement vectors (3*N x 3*N)
      integer, intent(in) :: element_numbers(:)
         !! Atomic numbers for each atom
      real(dp), intent(in), optional :: force_constants_mdyne(:)
         !! Force constants in mdyne/Å (if provided, these are printed instead)
      logical, intent(in), optional :: print_displacements
         !! If true, print Cartesian displacement vectors (default: true)
      integer, intent(in), optional :: n_atoms
         !! Number of atoms (if not provided, derived from size of element_numbers)
      real(dp), intent(in), optional :: ir_intensities(:)
         !! IR intensities in km/mol
      real(dp), intent(in), optional :: coordinates(:, :)
         !! Atomic coordinates (3, n_atoms) in Bohr - needed for thermochemistry
      real(dp), intent(in), optional :: electronic_energy
         !! Electronic energy in Hartree - needed for thermochemistry
      real(dp), intent(in), optional :: temperature
         !! Temperature for thermochemistry in K (default: 298.15)
      real(dp), intent(in), optional :: pressure
         !! Pressure for thermochemistry in atm (default: 1.0)

      integer :: n_modes, n_at, n_groups, igroup, imode, iatom, icoord, k
      integer :: mode_start, mode_end, modes_in_group
      logical :: do_print_disp
      character(len=512) :: line
      character(len=16) :: freq_str, mass_str, fc_str, ir_str
      character(len=2) :: elem_sym
      character(len=3) :: coord_label
      real(dp) :: fc_value

      n_modes = size(frequencies)
      if (present(n_atoms)) then
         n_at = n_atoms
      else
         n_at = size(element_numbers)
      end if

      do_print_disp = .true.
      if (present(print_displacements)) do_print_disp = print_displacements

      call logger%info(" ")
      call logger%info("============================================================")
      call logger%info("                  VIBRATIONAL ANALYSIS")
      call logger%info("============================================================")
      call logger%info(" ")

      ! Print in groups of 3 modes (like Gaussian)
      n_groups = (n_modes + 2)/3

      do igroup = 1, n_groups
         mode_start = (igroup - 1)*3 + 1
         mode_end = min(igroup*3, n_modes)
         modes_in_group = mode_end - mode_start + 1

         ! Mode numbers header
         line = "                    "
         do k = mode_start, mode_end
            write (freq_str, '(i12)') k
            line = trim(line)//freq_str
         end do
         call logger%info(trim(line))

         ! Frequencies (show "i" only for significant imaginary frequencies)
         line = " Frequencies --  "
         do k = mode_start, mode_end
            if (frequencies(k) < 0.0_dp .and. abs(frequencies(k)) > 10.0_dp) then
               ! Significant imaginary frequency - show with "i"
               write (freq_str, '(f12.4,a)') abs(frequencies(k)), "i"
            else
               ! Real or near-zero frequency
               write (freq_str, '(f12.4)') abs(frequencies(k))
            end if
            line = trim(line)//freq_str
         end do
         call logger%info(trim(line))

         ! Reduced masses
         line = " Red. masses --  "
         do k = mode_start, mode_end
            write (mass_str, '(f12.4)') reduced_masses(k)
            line = trim(line)//mass_str
         end do
         call logger%info(trim(line))

         ! Force constants
         if (present(force_constants_mdyne)) then
            line = " Frc consts  --  "
            do k = mode_start, mode_end
               write (fc_str, '(f12.4)') force_constants_mdyne(k)
               line = trim(line)//fc_str
            end do
         else
            line = " Frc consts  --  "
            do k = mode_start, mode_end
               write (fc_str, '(f12.6)') force_constants(k)
               line = trim(line)//fc_str
            end do
         end if
         call logger%info(trim(line))

         ! IR intensities (if provided)
         if (present(ir_intensities)) then
            line = " IR Intens  --  "
            do k = mode_start, mode_end
               write (ir_str, '(f12.4)') ir_intensities(k)
               line = trim(line)//ir_str
            end do
            call logger%info(trim(line))
         end if

         ! Cartesian displacements
         if (do_print_disp) then
          call logger%info(" Atom          X         Y         Z       X         Y         Z       X         Y         Z")

            do iatom = 1, n_at
               elem_sym = element_number_to_symbol(element_numbers(iatom))

               ! Build line with atom info and displacements for each mode
               write (line, '(i4,1x,a2)') iatom, elem_sym

               do k = mode_start, mode_end
                  do icoord = 1, 3
                     write (freq_str, '(f10.5)') cartesian_displacements(3*(iatom - 1) + icoord, k)
                     line = trim(line)//freq_str
                  end do
               end do
               call logger%info(trim(line))
            end do
         end if

         call logger%info(" ")
      end do

      ! Summary statistics
      call logger%info("------------------------------------------------------------")
      call logger%info(" Summary:")

      ! Count real vs imaginary frequencies
      block
         integer :: n_real, n_imag, n_zero
         real(dp) :: zero_thresh
         zero_thresh = 10.0_dp  ! frequencies below 10 cm⁻¹ considered "zero"

         n_real = 0
         n_imag = 0
         n_zero = 0
         do k = 1, n_modes
            if (abs(frequencies(k)) < zero_thresh) then
               n_zero = n_zero + 1
            else if (frequencies(k) < 0.0_dp) then
               n_imag = n_imag + 1
            else
               n_real = n_real + 1
            end if
         end do

         write (line, '(a,i5)') "   Total modes:              ", n_modes
         call logger%info(trim(line))
         write (line, '(a,i5)') "   Real frequencies:         ", n_real
         call logger%info(trim(line))
         write (line, '(a,i5)') "   Imaginary frequencies:    ", n_imag
         call logger%info(trim(line))
         write (line, '(a,i5)') "   Near-zero (trans/rot):    ", n_zero
         call logger%info(trim(line))

         if (n_imag > 0) then
            call logger%warning("System has imaginary frequencies - may be a transition state")
         end if
      end block

      call logger%info("============================================================")
      call logger%info(" ")

      ! Compute and print thermochemistry if coordinates and energy are provided
      if (present(coordinates) .and. present(electronic_energy)) then
         block
            type(thermochemistry_result_t) :: thermo_result
            call compute_thermochemistry(coordinates, element_numbers, frequencies, &
                                         n_at, n_modes, thermo_result, &
                                         temperature=temperature, pressure=pressure)
            call print_thermochemistry(thermo_result, electronic_energy)
         end block
      end if

   end subroutine print_vibrational_analysis