"This works for organic molecules, but what happens when the excited states are closer in energy to the ground state, for example in open-shell molecules or in atoms?"
If there's excited states close to the ground state, the approximation you said Gaussian uses, where excited state contributions are neglected, seems not to be such a great idea anymore. You would now want to include the close-lying excited states. There's plenty of methods that will give you excited state energies, for example EOM-CC and STEOM-CC.
"This is where everything becomes confusing to me. If I understand correctly, when the first or higher excited states are closer in energy to the ground state, the Hartree-Fock single determinant solution is not good enough."
When it comes to calculating energies, a single HF determinant is not good enough, even when the first excited state is far away from the ground state. This is because of the correlation energy missed by using HF alone.
"We use multi-configurational methods (such as CASSCF) where we generate excited determinants and consider the molecule's electronic state as a superposition of the different determinants (weighted by their energies?)."
It's true that CASSCF is an example of a multi-configurational method, and that the molecular state after a CASSCF calculation will be a superposition of various determinants, but they are not weighted by their energies, they are weighted by coefficients that are optimized by the variational principal.
There's more than one component of a CASSCF calculation, there's the CAS-CI (CAS-configuration-interaction) calculation where you calculate FCI (full configuration interaction) within an active space (for example all possible excitations from the HF determinant, involving $n$ electrons and $m$ orbitals); then you do an SCF iteration where you rotate the orbitals to get an even lower energy, then you do another CAS-CI calculation with the rotated orbitals, then you rotate the orbitals again and so on, until you get convergence to your desired precision.
"So this means the solution itself contains contributions from the low-lying excited states."
To get the ground state energy as accurately as possible from a variational calculation, you'll want your wavefunction ansatz to be as flexible as possible, which means including as many determinants as possible if you're modeling it as a linear combination of Slater determinants. Therefore, accurate ground state energies calculated this way, will use a "model" for the wavefunction which certainly contain more than one determinant.
"I might be wrong here, but I think the energy levels obtained from CASSCF calculations are obtained as an average of all the determinants,"
The word "average" is not the best word here. I'd say that the wavefunction is a linear combination of Slater determinants though. This is also not only true for CASSCF but also for other methods such as CISD which is considered a "single reference method".
"and have partial occupation numbers reflecting the excited state contributions."
The CASSCF (and CISD) density matrices will have non-zero occupation numbers for Slater determinants which include "excitations" from the Hartree-Fock (or reference) configuration.
"So, what happens when we do the thermochemistry calculation on a multi-configurational system? Does that mean the electronic partition function term is already included in the SCF energy?"
The phrase "multiconfigurational system" is maybe not the best to use here. All molecules have some multiconfigurational character when modeling the wavefunction as a sum of Slater determinants, unless perhaps you had some magical basis.
In your question you gave a formula for the partition function. All you need in order to evaluate that formula, are the energies and number of states at each energy level (the latter being what you call you the "degeneracies"). You can calculate all of the energies using something like EOM-CCSD or STEOM-CCSD (and this method will also allow you to determine the "degeneracies"). When you get the energies from a EOM-CCSD or STEOM-CCSD calculation, each of these energies will include contributions from several Slater determinants.
"Does this type of thermochemistry calculation actually make sense?"
I think you were over-complicating things for yourself, which lead to some confusion. The formula you presented is simple, and to evaluate it, you just need the energies from a method such as EOM-CCSD or STEOM-CCSD, which calculates several excited state energies. Each of these excited state energies will be calculated based on a wavefunction ansatz which involves more than one Slater determinant, but that's okay because even ground state calculations will need more than one Slater determinant if we want a decent accuracy (especially when "correlation energy" is a large contribution to the overall energy).
