import math


def harmonic_oscillator_heat_capacity(T, freq):
    """
    Return the heat capacity in J/mol*K at the given set of temperatures `Tlist`
    in K for the harmonic oscillator with a frequency `freq` in cm^-1.
    """
    x = freq / (0.695039 * T)  # kB = 0.695039 cm^-1/K
    exp_x = math.exp(x)
    one_minus_exp_x = 1.0 - exp_x
    return x * x * exp_x / one_minus_exp_x / one_minus_exp_x