  The variational inclusion of spin-orbit coupling in self-consistent field
(SCF) calculations requires a generalised two-component framework, which
permits the single-determinant wave function to completely break spin symmetry.
The individual components of the molecular g-tensor are commonly obtained from
separate SCF solutions that align the magnetic moment along one of the three
principal tensor axes. However, this strategy raises the question if energy
differences between solutions are relevant, or how convergence is achieved if
the principal axis system is not determined by molecular symmetry. The present
work resolves these issues by a simple two-step procedure akin to the generator
coordinate method (GCM). First, a few generalised Hartree Fock (GHF) solutions
are converged, applying, where needed, a constraint to the orientation of the
magnetic-moment or spin vector. Then, superpositions of GHF determinants are
formed through non-orthogonal configuration interaction. This procedure yields
a Kramers doublet for the calculation of the complete g-tensor. Alternatively,
for systems with weak spin-orbit effects, diagonalisation in a basis spanned by
spin rotations of a single GHF determinant affords qualitatively correct
g-tensors by eliminating errors related to spin contamination. For small
first-row molecules, these approaches are evaluated against experimental data
and full configuration interaction results. It is further demonstrated for two
systems (a fictitious tetrahedral CH4+ species, and a CuF4(2-) complex) that a
GCM strategy, in contrast to alternative mean-field methods, can correctly
describe the spin-orbit splitting of orbitally-degenerate ground states, which
causes large g-shifts and may lead to negative g-values.
