What happens to the complex part of the wave function when doing LCAO?
I understand that in general the wave function can take complex numbers [...]


It is not that it can, rather it should: by the very postulates of quantum mechanics wave function is a complex-valued function. So, both the spin orbitals (that are one-electron wave functions) and the Slater determinant built out of them are in principle complex-valued functions.

However, in solving the HF equations by the SCF procedure it is quite typical to impose some constraints on the spin orbitals, for instance, restrict them to be real-valued functions rather then complex-valued ones. As with any other constraint this restriction might lead to what is called SCF instabilities, i.e. situations when relaxing the constraints lead to a different variational solution of lower energy.

To elaborate a bit more let us mention few other well-known (as well as lesser-known) constraints for the SCF procedure:


The infamous restricted Hartree-Fock (RHF) model with the requirement that spin orbitals come in pairs: two spin orbitals corresponding to two different pure spin states are constructed from the same spatial orbital:
$$
\begin{aligned}
  \psi_{2i-1}(1) &amp;= \phi_{i}(1) \alpha(1) \, , \\ 
  \psi_{2i}(1)   &amp;= \phi_{i}(1) \beta(1)  \, .      
\end{aligned}
$$
In the unrestricted Hartree-Fock (UHF) model the requirement above is relaxed and we exclusively use spatial orbitals from one set to construct $\alpha$ spin orbitals and spatial orbitals from another set to construct $\beta$ spin orbitals:
$$
\begin{aligned}
  \psi_{2i-1}(1) &amp;= \phi_{i}^{\alpha}(1) \alpha(1) \, , \\ 
  \psi_{2i}(1)   &amp;= \phi_{i}^{\beta}(1) \beta(1)   \, .      
\end{aligned}
$$
A bit less known is the fact that the unrestricted Hartree-Fock method is not so unrestricted. In fact, there is a constraint here: each and every spin orbital describes an electron in a pure spin state, either $\alpha$ or $\beta$, while in general (in accordance with the postulates of quantum mechanics) electron can be in superposition of this states:
$$
  \psi_{i}(1) = \phi_{i}^{\alpha}(1) \alpha(1) + \phi_{i}^{\beta}(1) \beta(1)
$$
This most general setting is known as the general Hartree-Fock (GHF) method.


And all these three methods (RHF, UHF, GHF) exist in two variants: the real one, in which spin orbitals additionally required to be real-valued functions, and the more general complex one. All this give rise to six variants of the HF method with many possible instabilities between them discussed it some details in the seminal paper by Schlegel &amp; McDouall1:



There is also a very similar earlier paper by Seeger &amp; Pople.2 which suddenly is not freely available.

Its important to realize that any constraint of the mentioned above can only raise the electronic energy: with a great deal of certainty, we may expect that if any constraint is relaxed, the variational procedure will result in a lower energy due to greater variational freedom. Thus, for instance, for the most usual different real variants of the Hartree-Fock method we have
$$
    E_\mathrm{e}(\mathrm{RGHF}) \leq E_\mathrm{e}(\mathrm{RUHF}) \leq E_\mathrm{e}(\mathrm{RRHF}) \, .
$$
The same, of course, is true for the corresponding complex variants of these methods,
$$
    E_\mathrm{e}(\mathrm{CGHF}) \leq E_\mathrm{e}(\mathrm{CUHF}) \leq E_\mathrm{e}(\mathrm{CRHF}) \, .
$$
And for any of the three formalism 
$$
    E_\mathrm{e}(\mathrm{C}x\mathrm{HF}) \leq E_\mathrm{e}(\mathrm{R}x\mathrm{HF}) \, , \quad \text{where} \quad x = \mathrm{G}, \mathrm{U}, \mathrm{R} \, .
$$



RHF/UHF/GHF and correct spin symmetry

The possible rise of electronic energy that accompanies introduction of more and more constraints in the GHF -> UHF -> RHF sequence seems to be quite contrary to the goal of the variation method. However, UHF and RHF constraints are nothing but symmetry constraints: they arise by requiring an approximate electronic wave function to have the same spin symmetry as the exact non-relativistic one, i.e. to be an eigenfunction of spin operators, the total spin-squared operator $\hat{S}^2$ and the $z$-component of the total spin operator $\hat{S}_{z}$. 


With respect to $\hat{S}_{z}$ operator it is well known that any Slater determinant built out of spin orbitals corresponding to pure spin states is an eigenfunction of $\hat{S}_{z}$. Thus, RHF and UHF wave functions are eigenfunctions of $\hat{S}_{z}$, but GHF wave function is not.
The situation with $\hat{S}^{2}$ is a bit more complicated, but all in all only in RHF (but not in UHF or in GHF) formalism it is possible to construct an approximate electronic wave function which is an eigenfunction of $\hat{S}^{2}$.3


To conclude, as it was pointed out by Löwdin, we always face with a dilemma here: should we seek a solution that is a true variational minimum or should we seek a solution with the correct spin symmetry? 



1) H. B. Schlegel and J. J. W. McDouall, Do You Have SCF Stability and Convergence Problems? in Computational Advances in Organic Chemistry: Molecular Structure and Reactivity, Springer Netherlands, 1991, pp. 167-185. DOI: 10.1007/978-94-011-3262-6_2. Free PDF from wayne.edu.

2) Seeger, R., &amp; Pople, J. A., Self‐consistent molecular orbital methods. XVIII. Constraints and stability in Hartree–Fock theory. The Journal of Chemical Physics, 66(7), 1977, 3045-3050. DOI: 10.1063/1.434318.

3) This can be done even for an open-shell electron configuration, although it would require a linear combination of Slater determinants with constant coefficients, none of which alone is an eigenfunction of $\hat{S}^{2}$.