compute_partition_functions Subroutine

   pure subroutine compute_partition_functions(total_mass, moments, frequencies, n_freqs, &
                                               temperature, pressure, sigma, is_linear, &
                                               q_trans, q_rot, q_vib)
      !! Compute partition functions for translation, rotation, and vibration.
      real(dp), intent(in) :: total_mass        !! Total mass in amu
      real(dp), intent(in) :: moments(3)        !! Principal moments in amu*Angstrom^2
      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(in) :: pressure          !! Pressure in atm
      integer, intent(in) :: sigma              !! Symmetry number
      logical, intent(in) :: is_linear          !! True if linear molecule
      real(dp), intent(out) :: q_trans          !! Translational partition function
      real(dp), intent(out) :: q_rot            !! Rotational partition function
      real(dp), intent(out) :: q_vib            !! Vibrational partition function

      real(dp) :: mass_kg, T, P_pa
      real(dp) :: lambda, V_molar
      real(dp) :: theta_rot(3), u
      integer :: i

      T = temperature
      mass_kg = total_mass*AMU_TO_KG
      P_pa = pressure*ATM_TO_PA

      ! Translational partition function: q_trans = (2*pi*m*k*T/h^2)^(3/2) * V
      ! where V = kT/P for ideal gas (per molecule)
      lambda = H_SI/sqrt(2.0_dp*PI*mass_kg*KB_SI*T)  ! thermal de Broglie wavelength
      V_molar = KB_SI*T/P_pa  ! volume per molecule
      q_trans = V_molar/(lambda**3)

      ! Rotational partition function
      ! theta_rot = h^2 / (8*pi^2*I*k_B) = ROTTEMP_AMUA2_TO_K / I (for I in amu*Angstrom^2)
      do i = 1, 3
         if (moments(i) > 1.0e-6_dp) then
            theta_rot(i) = ROTTEMP_AMUA2_TO_K/moments(i)
         else
            theta_rot(i) = 0.0_dp
         end if
      end do

      if (is_linear) then
         ! Linear: q_rot = T / (sigma * theta_rot)
         if (theta_rot(3) > 0.0_dp) then
            q_rot = T/(real(sigma, dp)*theta_rot(3))
         else
            q_rot = 1.0_dp
         end if
      else
         ! Nonlinear: q_rot = sqrt(pi) * T^(3/2) / (sigma * sqrt(theta_A * theta_B * theta_C))
         if (theta_rot(1) > 0.0_dp .and. theta_rot(2) > 0.0_dp .and. theta_rot(3) > 0.0_dp) then
            q_rot = sqrt(PI)*(T**1.5_dp)/ &
                    (real(sigma, dp)*sqrt(theta_rot(1)*theta_rot(2)*theta_rot(3)))
         else
            q_rot = 1.0_dp
         end if
      end if

      ! Vibrational partition function: q_vib = Product_i [1 / (1 - exp(-u_i))]
      ! where u_i = h*nu/(k*T) = theta_vib/T = 1.4388 * nu(cm^-1) / T
      q_vib = 1.0_dp
      do i = 1, n_freqs
         if (frequencies(i) > 10.0_dp) then  ! Skip near-zero frequencies
            u = CM1_TO_KELVIN*frequencies(i)/T
            if (u < 100.0_dp) then  ! Avoid overflow
               q_vib = q_vib/(1.0_dp - exp(-u))
            end if
         end if
      end do

   end subroutine compute_partition_functions