How can I implement the Cartesian harmonics?

Question

I'm studying the quantum chemistry calculation and now trying to implement the basic methods.

For example, given a water molecule, $M = {(\ce{O}, R_\ce{O}), (\ce{H}, R_\ce{H}), (\ce{H}, R_\ce{H})}$, where $R = [X, Y, Z]$ is the atomic position, the Cartesian Gaussian-type orbital $\psi$ on $r$, where $r = [x, y, z]$ is an arbitrary position, can be written as follows: $$\psi(r) = N (x-X)^l (y-Y)^m (z-Z)^n \exp(-\alpha ||r-R||^2),$$ where $N$ is the normalization term, $\alpha$ is the orbital exponent, and $||r-R||$ is the distance between $r$ and $R$. The term of $(x-X)^l (y-Y)^m (z-Z)^n$ is the Cartesian harmonics.

However, I believe that, because a water molecule can be described by any Cartesian coordinates if all the atomic distances (i.e., $\ce{O-H}$, $\ce{O-H}$, and $\ce{H-H}$) are correct, $\psi(r)$ can not be determined to a single value.

For example, the following two Cartesian coordinates represent the same water molecule.

O 0.00000 0.00000 0.11779
H 0.00000 0.75545 -0.47116
H 0.00000 -0.75545 -0.47116

O -0.034 0.978 0.0076
H 0.065 0.021 0.0015
H 0.87 1.30 0.00069

However, $\psi(r)$ is different because the above equation has the term of $(x-X)^l (y-Y)^m (z-Z)^n$. If the equation has only the distance term (i.e., $\exp(-\alpha ||r-R||^2)$), $\psi(r)$ has always the same value and I can implement this very easily only using the distance matrix.

How can I deal with this problem and implement the Cartesian harmonics?

Postscript: Of course, I know the spherical harmonics but I believe that even if the angles, θ and φ, are used, the problem is the same; how can I determine the θ and φ or the orientation of the 2px (2py or 2pz) orbital of each atom in a molecule? For example, in the following figure, specifically how can I determine the 2px of the H orbital in a water molecule? Of course, I know the 2px (2py or 2pz) is along with the H-O bond in this case; however, how can I determine the θ, φ, and orientation of 2px, 2py, and 2pz orbitals of each atom in arbitrary other more complex molecule and calculate the LCAO?

Answer 1

The question is not very well-posed. However, the question appears to be about whether rotational invariance is honored by LCAO calculations.

The old-time approach is to use Cartesian basis functions. Translating our origin to the nucleus, we can write the Cartesian basis functions in the form $\psi({\bf r}) = N x^k y^l z^m \exp(-\alpha r^2)$, where $\lambda:=k+l+m$ is the angular momentum of the function.

Basis functions are added one shell at a time: there is one cartesian S function with $(k,l,m)=(0,0,0)$, three P functions with $(1,0,0)$, $(0,1,0)$, and $(0,0,1)$, and six D functions: $(2,0,0)$, $(0,2,0)$, $(0,0,2)$, $(1,1,0)$, $(1,0,1)$, $(0,1,1)$.

While Cartesian functions were used in the early days of quantum chemistry, and form part of e.g. the infamous 6-31G* basis set, they have mostly fallen out of use, since they contain angular degeneracies, the number of which increases rapidly in the angular momentum.

What are these degeneracies? We know that spherical harmonics $Y_l^m$, or equivalently their real counterparts $Y_{lm}$ form a complete set of angular distributions. For S states we have $Y_{00}$, for P states we have $Y_{1,-1}$, $Y_{1,0}$ and $Y_{1,1}$, and for D states we have $Y_{2,-2}$, $Y_{2,-1}$, $Y_{2,0}$, $Y_{2,1}$ and $Y_{2,2}$. That is, for S or P functions it doesn't matter if you use cartesians or pure spherical harmonics; but for D it starts to matter: we have 6 cartesians but only 5 pure angular functions.

If you write out $x$, $y$, and $z$ in spherical harmonics (you can find the expressions in the literature), you can use the closure properties to derive what the angular form for $x^k y^l z^m$ is. However, it's pretty simple to see that if you combine three of the six D functions $(2,0,0) + (0,2,0) + (0,0,2)$ you get $x^2 + y^2 + z^2 = r^2$ which is an S function, $Y_{00}$. This means that the cartesian D shell is formed of the 5 pure D functions, and a S-type contaminant.

In quantum chemical calculations, the pure basis turns out to be of the form $\psi({\bf r}) = N r^l Y_{lm} \exp(-\alpha r^2)$. This is very convenient, since one can express the solid harmonics $Y_{lm}$ in terms of the cartesian functions. You can find a table on Wikipedia, but it is also quite simple to program up the bare equations to get the translation for arbitrary angular momentum. The transformation has also been discussed by Schlegel and Frisch in Int. J. Quantum Chem. 54, 83 (1995).

So, finally.. why are LCAO electronic structure calculations rotationally invariant?$^*$ No matter how you rotate the molecule, the spherical harmonics span the same angular degrees of freedom: e.g. when you rotate $Y_{10}$ into a new coordinate system, it is still expressible in terms of $Y_{10}$, $Y_{1,-1}$ and $Y_{11}$; the expansion coefficients are given by the so-called Wigner's $D$ matrix. This means you will always get the same energy regardless of how you rotate your system; the individual components of the wave function will just change.

$^*$ Rotational invariance doesn't hold fully for e.g. density functional calculations, since the quadrature grids employed to integrate the density functional are not rotationally invariant. The error can, however, be made insignificantly small by using a large enough quadrature grid. This error is also removed by default by transforming the molecule into a standard orientation; see Chem. Phys. Lett. 209, 506 (1993).

Addendum based on postscript: sounds like your problem is just not understanding how the molecular orbitals are defined. You don't need to try to align the orbitals; they arise automatically by solving the self-consistent field equations. The molecular orbitals are linear combinations of the atomic orbitals, so they'll automatically figure out what is the best combination. There are no atomic orbitals in a molecule, they're just basis functions used in the expansion, the coefficients of which are found by minimizing the total energy of the resulting wave function.