# How can I implement the Cartesian harmonics?

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? • The wave function has translational and rotational invariance and you can therefore define it uniquely.
– Bertram
Mar 5 at 20:37
• My question is not the definition of the wave function; my question is how to calculate the same wave function from the above two different input data of Cartesian coordinates. Quantum chemical simulation softwares (e.g., Gaussian) can output the same result from two different input data if the data represent the same molecule. How can the simulation softwares deal with such input data how can I implement it? Mar 5 at 23:22
• It may be useful to remember that quantum chemistry suites use the basis functions mostly in integrals. Since the integration is done over all space (once or twice, depending on the integral), the orientation of the molecule does not matter (within the limits of numerical accuracy). Mar 6 at 14:15
• I see, but for example when we consider the LCAO, the molecular orbital ψ is calculated from the atomic basis functions φ as ψ(r) = c1φ1(r-R) + c2φ2(r-R) ... and then the electron density will be calculated as ρ(r) = sum |ψi(r)|^2. In these calculations, the integral does not appear. Additionally, if the atomic basis function φ has the terms of (x-X)^l(y-Y)^m(z-Z)^n that depend on the coordinates themselves (of course the distance does not depend on the coordinates), the density ρ(r) will be different. I may have some mistakes but I do not understand yet. Mar 7 at 1:04
• I found the following description in Stack Exchange Chemistry. chemistry.stackexchange.com/questions/73177/… When specifying a Gaussian basis set in quantum chemistry, one talks only about the radial part, independent of the angular part. The angular part is really considered an implementation detail, because its form usually doesn't have any significant effect on the results.'' Is it really? Mar 7 at 7:52

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.

• I edited the first question and added the figure. 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 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 is along with the H-O bond in this case; however, how can I determine the θ, φ, and orientation of 2px orbital in arbitrary other more complex molecules and calculate the LCAO? Mar 16 at 12:19
• @LinLin I've extended the answer. Shortly: you don't. There is no 2px H orbital in a water molecule. You just build your molecular basis set; a triple-zeta basis typically includes S, P, and D functions on the hydrogens and S, P, D, and F functions on the oxygen, and then you find the molecular orbitals that minimize the Hartree-Fock or density functional energy. Mar 16 at 16:21
• Thanks for answering this Susi! @LinLin if you want a longer answer, I would recommend that you ask your follow-up questions in separate posts. I think that Susi's 10 paragraph answer in this case is quite substantial and sufficiently answers the original question! Mar 16 at 18:34
• @Susi Lehtola, Thank you very much for the answer. It is very helpful and I keep studying this topic. Mar 17 at 1:37