compute_rotational_thermo Subroutine

   pure subroutine compute_rotational_thermo(moments, temperature, symmetry_number, is_linear, E, S, Cv)
      !! Compute rotational contributions to thermodynamic properties.
      !!
      !! Uses rigid rotor approximation (classical limit, high T).
      !! For nonlinear: E = 3/2 RT, Cv = 3/2 R
      !! For linear: E = RT, Cv = R
      real(dp), intent(in) :: moments(3)      !! Principal moments in amu*Angstrom^2
      real(dp), intent(in) :: temperature     !! Temperature in K
      integer, intent(in) :: symmetry_number  !! Rotational symmetry number
      logical, intent(in) :: is_linear        !! True if linear molecule
      real(dp), intent(out) :: E              !! Thermal energy in Hartree
      real(dp), intent(out) :: S              !! Entropy in cal/(mol*K)
      real(dp), intent(out) :: Cv             !! Heat capacity in cal/(mol*K)

      real(dp) :: theta_rot(3)
      real(dp) :: T, sigma
      real(dp) :: qr
      integer :: i

      T = temperature
      sigma = real(symmetry_number, dp)

      ! Calculate rotational temperatures: theta_rot = h^2 / (8*pi^2*I*k_B)
      ! For I in amu*Angstrom^2:
      !   I_SI = I * 1.66054e-47 kg*m^2
      !   theta_rot = h^2 / (8*pi^2 * I_SI * k_B)
      !             = (6.62607e-34)^2 / (8 * pi^2 * 1.66054e-47 * 1.38065e-23 * I)
      !             = ROTTEMP_AMUA2_TO_K / I  [K]
      do i = 1, 3
         if (moments(i) > LINEAR_THRESHOLD) 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 molecule: 2 rotational degrees of freedom
         E = R_HARTREE*T
         Cv = R_CALMOLK

         ! S = R * [1 + ln(T / (sigma * theta_rot))]
         ! Use the average of the two non-zero theta values
         if (theta_rot(3) > 0.0_dp) then
            qr = T/(sigma*theta_rot(3))
            S = R_CALMOLK*(1.0_dp + log(qr))
         else
            S = 0.0_dp
         end if
      else
         ! Nonlinear molecule: 3 rotational degrees of freedom
         E = 1.5_dp*R_HARTREE*T
         Cv = 1.5_dp*R_CALMOLK

         ! S = R * [3/2 + ln(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
            qr = sqrt(PI)*(T**1.5_dp)/(sigma*sqrt(theta_rot(1)*theta_rot(2)*theta_rot(3)))
            S = R_CALMOLK*(1.5_dp + log(qr))
         else
            S = 0.0_dp
         end if
      end if

   end subroutine compute_rotational_thermo