Do LCAO wavefunctions form a basis for states of diatomic molecules? I've been trying to gain understanding for diatomic molecules by visualizing the electron wavefunctions of their ground and excited states. As, for this purpose, quantitatively correct energies, dipole moments, etc., are not so important, I decided to start with a very simple case - molecular orbitals calculated with LCAO for a one-electron molecule, $H_2^+$ or $HD^+$ (where $D$ stands for deuterium).
To be specific, let's say we have hydrogen nuclei $a, b$, which can be in the atomic states $\psi_{1s}, \psi_{2s}, \psi_{2p}$, etc. Then, in the example of the ground state of $H_2^+$, we have $\psi = c_1 \psi_{1s}^a + c_2 \psi_{1s}^b$. In order to conform to the symmetries of the problem, $c_1 = \pm c_2$ and we get the bonding and antibonding $\sigma_{1s}, \sigma_{1s}^*$ molecular orbitals. This makes sense.
However, going to higher excited states, I'm starting to have trouble. For example, I think that the combination $\psi = c_1 \psi_{1s}^a + c_2\psi_{2p}^b$ does not exist as it would violate inversion symmetry, but is this always the case when combining two different atomic states?
Going to $HD^+$, I get even more confused. I think that a combination of two different orbitals like $\psi = c_1 \psi_{1s}^a + c_2\psi_{2s}^b$ now seems reasonable, because inversion symmetry is broken, but what about $\psi = c_1 \psi_{1s}^a + c_2\psi_{2p}^b$? Naively, it seems that this combination would have to violate rotational symmetry around the internuclear axis, but couldn't there be some superposition of such states that fixes this problem?
So: How do I know which linear combinations are allowed and which are not, and how do I know that I have taken enough combinations into account to actually cover all the existing excited states? How can I make sure that my molecular orbitals indeed form a basis of (approximate) eigenstates of the potential of the two nuclei, both in the case of homonuclear and heteronuclear diatomic molecules?
 A: LCAO is essentially a variational method, where all we are doing is guessing the answer and calculating the energies that result.  These energies may or may not be close to the actual experimental answers, which we should always refer to.

Let me briefly derive the answer for $\text{H}_2^+$ using the full method... As you say, we start by saying we would like a solution of the form $\psi=c_1\psi_\text{1s}^a+c_2\psi_\text{1s}^b$.  In the variational method, we say that the quantity given by
$$\widetilde{E}=\frac{\int\psi^*\hat{H}\psi\;\text{d}\tau}{\int\psi^*\psi\;\text{d}\tau}$$
is the energy of our wavefunction. $\hat{H}$ is the Hamiltonian for the system, and $\text{d}\tau$ indicates we should integrate over all space.  Sticking in $\psi$ and after a bit of rearranging, we get:
$$\widetilde{E}=\frac{c_1^2\int\psi_\text{1s}^a\hat{H}\psi_\text{1s}^a\;\text{d}\tau+c_2^2\int\psi_\text{1s}^b\hat{H}\psi_\text{1s}^b\;\text{d}\tau+2c_1c_2\int\psi_\text{1s}^a\hat{H}\psi_\text{1s}^b\;\text{d}\tau}{c_1^2\int\psi_\text{1s}^a\psi_\text{1s}^a\;\text{d}\tau+c_2^2\int\psi_\text{1s}^b\psi_\text{1s}^b\;\text{d}\tau+2c_1c_2\int\psi_\text{1s}^a\psi_\text{1s}^b\;\text{d}\tau}$$
Note I've assumed that everything is real, and that for example $\int\psi_\text{1s}^b\hat{H}\psi_\text{1s}^a\;\text{d}\tau=\int\psi_\text{1s}^a\hat{H}\psi_\text{1s}^b\;\text{d}\tau$.  Let's define some constants to make the integrals go away:
\begin{align}
\alpha&=\int\psi_\text{1s}^a\hat{H}\psi_\text{1s}^a\;\text{d}\tau=\int\psi_\text{1s}^b\hat{H}\psi_\text{1s}^b\;\text{d}\tau\\
\beta&=\int\psi_\text{1s}^b\hat{H}\psi_\text{1s}^a\;\text{d}\tau=\int\psi_\text{1s}^a\hat{H}\psi_\text{1s}^b\;\text{d}\tau\\
S&=\int\psi_\text{1s}^a\psi_\text{1s}^b\;\text{d}\tau
\end{align}
These can be interpreted as the self-energy ($\alpha$), the interaction energy ($\beta$) and the orbital overlap ($S$).  This leaves us with
$$\widetilde{E}=\frac{(c_1^2+c_2^2)\alpha+2c_1c_2\beta}{c_1^2+c_2^2+2c_1c_2S}$$
Let's assume for simplicity that $S=0$:
$$\widetilde{E}=\frac{(c_1^2+c_2^2)\alpha+2c_1c_2\beta}{c_1^2+c_2^2}$$
We now have to minimise $\widetilde{E}$ to get an upper bound on the ground state energy of the system described by $\hat{H}$ - this is what the variational principle says.  To do this, we must differentiate with respect to $c_1$ and $c_2$, and after doing a bit more algebra and setting $\frac{\partial\widetilde{E}}{\partial c_1}=\frac{\partial\widetilde{E}}{\partial c_2}=0$, we get what are called the secular equations:
$$\begin{pmatrix}\alpha-\widetilde{E}&\beta\\\beta&\alpha-\widetilde{E}\end{pmatrix}\begin{pmatrix}c_1\\c_2\end{pmatrix}=\begin{pmatrix}0\\0\end{pmatrix}$$
Hopefully it's not too difficult to see that this has non-zero solutions for $c_1$ and $c_2$ when $c_1=\pm c_2$, so long as $\widetilde{E}=\alpha\pm\beta$ matching your expectation.

Don't worry if you didn't follow all of that.  What's important is my introduction of this quantity $\beta$, which is the interaction strength between my two orbitals.  For two 1s orbitals, we could expect $\beta$ to be a good non-zero value.  For the case of $\text{H}$ with $\text{D}$, we might expect $\beta$ to be a little different as the atoms are slightly different.  Note that I think you should still be combining 1s and 1s orbitals together - they're the closest in energy and so will have the biggest overlap, $\beta$.  However, there is nothing stopping you from proposing a wavefunction that looks like $\psi=c_1\psi_\text{1s}^a+c_2\psi_\text{1s}^b+c_3\psi_\text{2s}^a+c_4\psi_\text{2s}^b$ and so on.  You'll just have a much bigger matrix to worry about.  The number of molecular orbitals you get out equals the number of atomic orbitals you put in.  In the case of using the above 4 orbitals, the system of equations you'll have to solve will look like
$$\begin{pmatrix}\alpha_\text{1s}-\widetilde{E}&\beta_\text{1s-1s}&\beta_\text{1s-2s}&\beta_\text{1s-2s}\\\beta_\text{1s-1s}&\alpha_\text{1s}-\widetilde{E}&\beta_\text{1s-2s}&\beta_\text{1s-2s}\\\beta_\text{1s-2s}&\beta_\text{1s-2s}&\alpha_\text{2s}-\widetilde{E}&\beta_\text{2s-2s}\\\beta_\text{1s-2s}&\beta_\text{1s-2s}&\beta_\text{2s-2s}&\alpha_\text{2s}-\widetilde{E}\end{pmatrix}\begin{pmatrix}c_1\\c_2\\c_3\\c_4\end{pmatrix}=\begin{pmatrix}0\\0\\0\\0\end{pmatrix}$$
Note here the $\beta$s are different, because we have different overlaps.  I suspect you'll find that $\beta_\text{1s-1s}$ and $\beta_\text{2s-2s}$ are strong, but $\beta_\text{1s-2s}$ is weak (since it's not a good energy match).  Also note that $\alpha_\text{2s}$ will be a lot bigger than $\alpha_\text{1s}$.  All of this will lead to you finding that the coefficients for the lower molecular orbitals will be very heavily weighted to the 1s orbitals only.

You can indeed also include 2p, 3d, 4f whatever contributions, but for the same reasons as above, you will find that orbitals wildly away in energy from the others will have large $\alpha$ and small $\beta$, and consequently will not show up very much in molecular orbitals you care about (the lower ones).  Further, if you wanted to try combining for example $\text{2s}$ and $\text{2p}_x$ along the $z$-axis, you would find that their overlap $\beta$ would be $0$, by symmetry, and so you would not get a $\text{2p-2s}$ contribution to your molecular orbitals (it'll just pop out as a non bonding orbital if nothing else is there to interact with).

So to answer your questions...

How do I know which linear combinations are allowed and which are not...

You calculate the overlap integral $$\beta=\int\psi_1\hat{H}\psi_2\;\text{d}\tau$$ which tells you mathematically which orbitals can interact and which cannot.

...how do I know that I have taken enough combinations into account to actually cover all the existing excited states?

You calculate the energy of the wavefunction you've already tried, and you compare it to experiment.  If you're close, you're on the right lines.  If you're not, you have to look at your model and think about what you've missed / approximated / simplified.  Yes, often adding more orbitals to your basis set gets you a more accurate answer, but that comes at the expense of your poor computer dealing with large matrices.
