Should one consider the spatial degeneracies of the electronic state while calculating the Gibbs free energy?

Question

The Gibbs free energy of a molecule is generally expressed as a sum of translational, rotational, vibrational, electronic and nuclear contributions. The electronic contribution $G_{elec}$ is formally a sum of Boltzmann factors of all possible electronic states of the system, but since the first electronic excitation energy is usually much larger than $kT$, many programs make the approximation $G_{elec} = -RT \ln g_0$, where $g_0$ is the degree of degeneracy of the ground state. For non-relativistic calculations of spatially non-degenerate states, $g_0$ is simply the spin multiplicity $2S+1$. In this case the program can easily choose the correct $g_0$ automatically, without intervention by the user.

However, during our recent development of the BDF software (http://182.92.69.169:7226/Introduction), I realized that the determination of $g_0$ in BDF is much more complicated than simply equating it to $2S+1$. The reason is that BDF is known for its accurate relativistic methods and full consideration of spatial symmetry, including double-group symmetry. Consequently, compared to other programs, a disproportionately high number of BDF users routinely perform calculations that take into account spin-orbit coupling, where the degree of degeneracy is $2J+1$ instead of $2S+1$, and a considerable number of users also calculate spatially degenerate electronic states.

My question is about the latter case. If the electronic state is spatially degenerate, say we are calculating the Gibbs free energy of the triply degenerate $1^3P$ state of the carbon atom, should we also take the spatial degeneracy of the state into account (i.e. in this case we have $g_0=9$ instead of $g_0=3$)?

Theoretically it seems that the answer is yes, but if the electronic structure method used in calculating the electronic energy of the system breaks the spatial symmetry, then the consideration of spatial symmetry may lead to inconsistencies. For example, if the aforementioned carbon atom is computed at the DFT level, then the spherical rotational symmetry of the atom is broken, so that if we compute the excited states of the carbon atom at the same level using TDDFT, we will find that the other two components of the triply degenerate $1^3P$ state lie quite a bit above the ground state, and a sum of Boltzmann factors will yield a nearly zero spatial contribution to the electronic Gibbs free energy, instead of the expected $-RT \ln 3$. Therefore, choosing $g_0=9$ in this case amounts to taking into account a degeneracy that exists in reality but does not exist in the approximate computational method that we are using, which does not sound entirely satisfactory. Moreover, the consideration of spatial symmetry in this case makes the automatic determination of $g_0$ very difficult, as the program has to deduce that the true spectral term of the carbon atom is $1^3P$, even though the DFT wavefunction has only $D_{\infty h}$ symmetry.

On the practical side, the latter problem means that in all quantum chemistry programs that I know of, a $1^3P$ carbon atom computed at the DFT level is treated as if it is spatially non-degenerate (at least the degeneracy does not enter the partition function). But I don't know if the user is supposed to manually correct for this. If no, then what about the calculations where the atom is treated by a method that does respect the spatial symmetry, like state-averaged CASSCF? To the best of my knowledge, many (if not all) quantum chemistry programs do not take into account the spatial degree of degeneracy even in this case, but again I don't know if the user has the responsibility to manually add the spatial contribution to the electronic Gibbs free energy.

Thank you very much in advance.

Answer 1

First, let's look from the most fundamental point of view. For this part of the answer I'll be referencing McQuarrie's Statistical Mechanics [1]. As the Gibbs free energy, $G$, can be calculated from the molecular partition function, $q_{\text {molecule }}$, we will start with $q_{\text {molecule }}$:

$$q_{\text {molecule }}=q_{\text {translational }} q_{\text {rotational }} q_{\text {vibrational }} q_{\text {electronic }} \cdots$$

Here, we assume that molecular Hamiltonian can be approximated by a sum of Hamiltonians for the various degrees of freedom of the molecule (which is not always true).

Now, focusing on the electronic partition function, we have $$q_{\mathrm{elect}}=\sum \omega_{e i} e^{-\beta \varepsilon_{i}} \quad\quad(1)$$ where $\omega_{e i}$ is the degeneracy, and $\varepsilon_i$ the energy of the $i$th electronic level. $\beta = 1/k_\mathrm{B}T$. Setting the arbitrary zero of energy such that $\varepsilon_1 = 0$, the electronic contribution to $q_{\text {molecule }}$ can then be written as $$q_{\mathrm{elect}}=\omega_{e 1}+\omega_{e 2} e^{-\beta \Delta \varepsilon_{12}}+\cdots$$ where $\Delta \varepsilon_{1 j}$ is the energy of the $j$th electronic level relative to the ground state. These are typically of the order of electron volts, so $\beta\Delta\varepsilon$ is typically large at ordinary temperature and only the first term in the summation for $q_{\mathrm{elect}}$ is significantly different from zero.

As your question mentions the carbon atom and as electronic energies of atoms are well tabulated [2], let's take a closer look at that.

atom	electronic configuration	term symbol	degeneracy $g = 2J + 1$	energy (cm$^{-1})$
C	2s$^2$2p$^2$	$^3P_0$	1	0.0000000
		$^3P_1$	3	16.4167130
		$^3P_2$	5	43.4134567

Using eq. 1, we can calculate $q_{\mathrm{elect}}$ at different temperatures:

$T(\mathrm{K})$	$q_{\mathrm{elect}}(T)$
0	1
298.15	7.8264551
2000	8.811036
20000	8.9808679

We can see that at absolute zero the partition function is the degeneracy of the ground state and that it increases as the temperature increases, going to the limit of 9 at infinite temperature as all available states would be accessible and the spin-orbit coupling would be negligible (note that we would have to include more states in eq. 1 at high temperatures). If one neglected spin-orbit coupling from the beginning, the result would be $$ q_{\mathrm{elect}}(T) = 9 = (2S+1)(2L+1) = \sum_J (2J+1)$$ for any temperature.

Now, one has to think how to make a general implementation of something like this in an electronic structure code. And here is my opinion: you don't.

If one does a non relativistic calculation for the ground state of the carbon atom, assuming that one has to work with an Abelian point group, one can use the $D_{2h}$ point group where a $P$ state would have components in the $B_{1u}$, $B_{2u}$, and $B_{3u}$ irreducible representations. Individually, each single one of them is non-degenerate, so I don't think that there is a way for the code to know that it must account for other components that it is not computing. Therefore, it seems reasonable to consider only the spin part as other codes do, for example Gaussian [3].

If the user wants to include excited states or consider the full degeneracy of the electronic ground state in a thermochemical calculation, it must be dealt case by case, manually. If one has near degenerate electronic states, including a few terms in the electronic partition function is the least of their problems.

One thing that you can to in your code to facilite for the user is to allow them to override the standard $g_0 = (2S+1)$ in an input option. For example, if I want to use $g_0 = 7.8264551$ in my thermochemical calculation for the carbon atom, I could provide that in my input.

And of course, everything could be well documented in the manual and, even, in the output file.

Here are some examples of rate constant calculations that make use of experimental data to include the spin-orbit coupling in the calculation of the electronic partition function of some reactants.

In the study of the O($^3P$) + HBr reaction [4], the potential energy surface was calculated at the MRCI/CBS+SO level of theory and the electronic partition function of the oxygen atom, in its ground state was calculated as $$ q_\mathrm{el} = 5+3 \mathrm{e}^{-227.8 / \mathrm{T}}+\mathrm{e}^{-326.6 / T}$$.

In the study of the OH($^2\Pi$) + HBr [5], the potential energy surface was determined at the UCCSD(T)/cc-pVDZ-F12a level of theory and the electronic partition function of the hydroxyl radical, in its ground state was calculated as $$g_{\mathrm{el}}(T)=2\left(1+\mathrm{e}^{-|A| h c / k_{\mathrm{B}} T}\right)$$ where $A$ is the experimental spin-orbit coupling constant.

In the context of chemical kinetics, this procedure is related to the fact that multiple potential energy surfaces arise from the reaction of species in degenerate electronic states and only a fraction of the collisions will follow the ground state surface, which is the most relevant for low energy reactive scattering [6].

1. McQuarrie, Donald A. Statistical mechanics. New York: Harper & Row, 1975.
2. Kramida, A., Ralchenko, Yu., Reader, J., and NIST ASD Team (2021). NIST Atomic Spectra Database (ver. 5.9), [Online]. Available: https://physics.nist.gov/asd [2021, December 5]. National Institute of Standards and Technology, Gaithersburg, MD. DOI: https://doi.org/10.18434/T4W30F
3. https://gaussian.com/thermo/
4. A. G. S. de Oliveira-Filho, F. R. Ornellas, K. A. Peterson, S. L. Mielke, Thermal Rate Constants for the O(3P) + HBr and O(3P) + DBr Reactions: Transition-State Theory and Quantum Mechanical Calculations. J. Phys. Chem. A. 117, 12703–12710 (2013).
5. A. G. S. de Oliveira-Filho, F. R. Ornellas, J. M. Bowman, Energy Disposal and Thermal Rate Constants for the OH + HBr and OH + DBr Reactions: Quasiclassical Trajectory Calculations on an Accurate Potential Energy Surface. J. Phys. Chem. A. 118, 12080–12088 (2014).
6. D. G. Truhlar, Multiple Potential Energy Surfaces for Reactions of Species in Degenerate Electronic States, J. Chem. Phys. 56, 3189–3190 (1972).