compute_vibrational_thermo Subroutine

   pure subroutine compute_vibrational_thermo(frequencies, n_freqs, temperature, E, S, Cv)
      !! Compute vibrational contributions to thermodynamic properties.
      !!
      !! Uses quantum harmonic oscillator partition function.
      !! E_vib = R * sum(theta_v * [1/(exp(u)-1)]) where u = theta_v/T
      !! S_vib = R * sum([u/(exp(u)-1) - ln(1-exp(-u))])
      !! Cv_vib = R * sum(u^2 * exp(u) / (exp(u)-1)^2)
      !!
      !! Note: ZPE is NOT included here (computed separately).
      real(dp), intent(in) :: frequencies(:)  !! Frequencies in cm^-1
      integer, intent(in) :: n_freqs          !! Number of frequencies
      real(dp), intent(in) :: temperature     !! Temperature in K
      real(dp), intent(out) :: E              !! Thermal energy in Hartree (excluding ZPE)
      real(dp), intent(out) :: S              !! Entropy in cal/(mol*K)
      real(dp), intent(out) :: Cv             !! Heat capacity in cal/(mol*K)

      integer :: i
      real(dp) :: T, freq, theta_v, u
      real(dp) :: exp_u, exp_neg_u
      real(dp) :: E_sum, S_sum, Cv_sum

      T = temperature
      E_sum = 0.0_dp
      S_sum = 0.0_dp
      Cv_sum = 0.0_dp

      do i = 1, n_freqs
         freq = frequencies(i)

         ! Skip imaginary and near-zero frequencies
         if (freq <= IMAG_FREQ_THRESHOLD) cycle
         if (freq < 10.0_dp) cycle  ! Skip very low frequencies (likely trans/rot residuals)

         ! Vibrational temperature: theta_v = h*c*nu / k = 1.4388 * nu (cm^-1)
         theta_v = CM1_TO_KELVIN*freq

         ! Reduced temperature ratio
         u = theta_v/T

         ! Avoid numerical issues for very large u (very low T or high freq)
         if (u > VIB_CLASSICAL_LIMIT) then
            ! Classical limit: modes are frozen out
            cycle
         end if

         exp_u = exp(u)
         exp_neg_u = exp(-u)

         ! Energy contribution (excluding ZPE): theta_v / (exp(u) - 1)
         E_sum = E_sum + theta_v/(exp_u - 1.0_dp)

         ! Entropy contribution: u/(exp(u)-1) - ln(1-exp(-u))
         S_sum = S_sum + u/(exp_u - 1.0_dp) - log(1.0_dp - exp_neg_u)

         ! Heat capacity contribution: u^2 * exp(u) / (exp(u)-1)^2
         Cv_sum = Cv_sum + (u*u*exp_u)/((exp_u - 1.0_dp)**2)
      end do

      ! Convert to proper units
      E = KB_HARTREE*E_sum           ! Hartree
      S = R_CALMOLK*S_sum            ! cal/(mol*K)
      Cv = R_CALMOLK*Cv_sum          ! cal/(mol*K)

   end subroutine compute_vibrational_thermo