Consider a Slater determinant of a two-electron system with spin orbitals $\phi$ and $\overline{\phi}$, i.e.
$$ \Psi = \frac{1}{\sqrt{2}}\left| \begin{matrix} \phi_{1s}(x_1) & \overline{\phi_{1s}}(x_1) \\ \phi_{1s}(x_2) & \overline{\phi_{1s}}(x_2) \\ \end{matrix} \right| $$
with an overbar indicating $\beta$ spin, and no overbar $\alpha$, and both electrons occupy a 1s state.
I am interested in simply calculating a value for $\Psi$ at some specific set of electron coordinates $x_1$ and $x_2$. However I must somehow deal with the spin functions contained in these orbitals, which to my understanding aren't really functions after all but rather symbols to indicate orthogonality. (Disclaimer: This question is coming from a background of quantum chemistry where concepts and formalisms may not be taught as stringently as considered necessary from a physicist's perspective.)
Apparently, in Quantum Monte Carlo the wavefunction is split into spin-up and spin-down Slater determinants, $\Psi = \psi^{\alpha}\psi^{\beta}$. (See also this unanswered question.) How exactly this is done? Does $\psi^{\alpha}$ contain only the spin orbitals with $\alpha$ spin, $\psi^{\beta}$ only those with $\beta$ spin? Then our two-electron wavefunction would be $\Psi = \phi_1(x_1) \overline{\phi_1}(x_2)$, which does not obey the Pauli principle.
So, how would one calculate numeric values of wavefunctions containing spin orbitals at specific sets of electron coordinates? Is this possible at all?