I am not so certain about the terminology because I do not use electronic structure software, but since nobody has answered so far, I will try to see if my answer clarifies something.
In a closed-shell configuration you know that a single <PERSON> determinant built from atomic orbitals (or for that matter, any function with angular dependence in terms of spherical harmonics) will be an eigenstate of $\hat{L}^2$ and $\hat{S}^2$ with quantum numbers L = S = 0. Additionally in this case there is no spin-orbit coupling, so that J = 0 as well. Therefore the angular momentum properties of this wave function will be well defined.
In the simplest ab initio molecular orbital calculation you obtain your atomic orbitals solving the restricted Hartree-Fock equations (RHF). Why restricted? Because if you have n electrons and the dimension of your basis (whatever STO, GTO or more advance basis, assuming you expand your atomic orbitals in a basis) is L > n/2, you populate the first n/2 orbitals and assume each orbital is doubly occupied by two electrons with spins $\uparrow$ and $\downarrow$. You could use the unrestricted Hartree-Fock (UHF) method, populating the lowest n spin-orbitals, but then the density of electrons with spin $\uparrow$ will be in general different than the density of electrons with spin $\downarrow$. Thus, although the energy of your system will be lower (some restrictions in the variational procedure are lifted) $S$ will not be well defined.
In general, for an open shell configuration the angular momentum will not be well defined by a single Slater determinant wave function. How do you create the proper wave function? You can use several strategies.
The more general strategy is to use a Multi-Configurational Self Consistent Field (MCSCF) approach, in which you generate several Slater determinants with the L orbitals, in which both the orbitals and the coefficients multiplying the Slater determinants are variational parameters that you have to determine consistently within the MCSCF approach.
A simpler method including correlation (among the many that exists) is the Configurational Interaction (CI) in which you optimize only once the orbitals and then generate the <PERSON> determinants and then at the end of the day obtain the coefficients of these <PERSON> determinants by diagonalizing the Hamiltonian.
The previous methods are called method with dynamic correlation or beyond the HF approach. More simply, you can use the RHF scheme and after obtaining your n/2 populated orbitals force the symmetry by suitable projector operators that generate the other <PERSON> determinants needed so that your overall wave function is an eigenstate of $\hat{L}^2$ and $\hat{S}^2$ (and $\hat{L}_z$ and $\hat{S}_z$). This is good only in the so-called L-S coupling scheme when spin-orbit coupling is assumed small (small relativistic corrections). 976
But of course you can later diagonalize the Hamiltonian with the spin-orbit coupling terms to resolve the fine-structure of the atom. I believe this is what GAMESS calls the ROHF method but I am not sure. To be sure you have to look at the manual.
More commonly, you can use the UHF method and after obtaining you n populated spin-orbitals you force the right symmetry by applying projector operators to filter the "wrong" components that give an incorrect spin symmetry. This method is often called the PHF or Projected Hartree-Fock. Finally, you can also use the RHF equations but start from a number of Slater determinants with the coefficients fixed by the projector operator (that is, they are no longer variational parameters) and then variationally optimize the orbitals for this set of determinants. The equations become more difficult and I believe this method, called EHF or Extended Hartree-Fock is not often applied. If you want to work with several <PERSON> determinants you typically use CI or MCSCF.
Now going back to your examples.
In Lithium you have two microstates for the ground state that you can write assigning the 4 quantum numbers $n_j$, $l_j$, $m_j$ and $m_{s,j}$ for each electron $j$, assuming the usual <PERSON> principle. However these two states are degenerate. Therefore you can write two <PERSON> determinants which will give you the same ground state energy. In the absence of an external magnetic field, any of these <PERSON> determinants (or whatever linear combination of the two) is a good description of your ground state. 976