If you're a masochist like me, and you're writing or have written a quantum chemistry program, you'll want to rotate molecular orbital (MO) coefficients to save yourself recalculating wavefunctions at different positions/orientations. Now that I've sufficiently reduced the interested audience for this to the point we could fit in the seats of a social-distanced sedan, let me take you on a painful but hopefully helpful journey.
Throughout this I'm going to assume that all rotations are about some common origin, , but of course you're welcome to double the pain and mix in rotations about arbitrary points in space if you wish to be more general. I'm also going to put aside the various different ordering conventions, which will probably make you pull your hair out if you try to interoperate between many different programs.
But first we need a quick refresher on the difference between Cartesian and spherical/pure Gaussian basis sets like those used in quantum chemistry. Note that throughout this I'm going to ignore elaboration of the different normalizations within contracted Gaussians, because frankly that just makes the equations more verbose...
Contracted Gaussian basis setsโ
Contracted Gaussian functions are of the form
where is the degree of contraction (i.e. number of primitive Gaussian functions), is the normalization constant we're going to ignore for now, is the exponent for each primitive basis function and the contraction coefficient.
Typically these sorts of basis sets will be in some library read in by the QM program of choice,
e.g. this is an excerpt from the STO-3G basis set in .gbs
format, with the left column being
the , and the right column being
H 0
S 3 1.00
3.42525091 0.15432897
0.62391373 0.53532814
0.16885540 0.44463454
With the exception of the term, is purely radial. However, there are two common choices of : Cartesian and spherical.
Cartesian polynomialsโ
In a Cartesian Gaussian basis set, takes the following form based on angular momentum :
where and all of are whole numbers.
Because multiplication is commutative i.e. , for each angular momenta, , there are a number of 'unique' functions e.g.
Label | Angular momentum | # | Components |
---|---|---|---|
3 | , . | ||
6 | , , , , | ||
10 | , , , , , , , , , |
and so on...
Spherical/Pure polynomialsโ
In a pure (or spherical) Gaussian basis set, usually corresponds to the real regular solid harmonics, which are directly related to the usual spherical harmonics These are a bit more of a rabbit hole..
Spherical harmonics ()
A much better overview of the spherical harmonics and their derivation can be found on wikipedia or in any number of other places.
Nevertheless, they take the form:
where is yet another different normalization factor, and are the polar angle/colatitude and the azimuthal angle/longitude respectively, and is an associated Legendre Polynomial
Regular solid harmonics ()
Given the previous definition of the spherical harmonics, the regular solid harmonics take the form:
or, simplified:
Real regular solid harmonics ()
Finally, the functions that are actually used are of the form:
where
This means that there are components for a pure Gaussian basis with angular momentum e.g.:
Label | Angular momentum | # | Components |
---|---|---|---|
3 | , , | ||
5 | , , , , | ||
7 | , , , , , , |
This should give a good idea why they find use in QM programs (fewer basis functions = fewer integrals to compute).
Rotating Cartesian Gaussiansโ
One way to imagine the molecular orbital coefficients associated with a Cartesian Gaussian is as an rank tensor e.g. a matrix for a function, a cube for an etc.
This obviously is very wasteful, as only the upper triangle (or pyramid) is permutationally unique. As such, they're typically stored as vectors (or columns of a matrix) of length
This usually renders rotation of these coefficients in real space to be done by hand coded routines rather than a simple matrix multiply.
However, we can convert the rotation matrix into an appropriate matrix to perform the conversion via the following algorithm implemented in C++:
/*
* Result should be R: an MxM rotation matrix for
* P: a MxN set of coordinates
* giving results P' = R P
*
* Assume that the poorly named power_index_arrays gives
* arrays of length l that look like the following:
* xy = {0, 1}
* xyz = {0, 1, 2}
* xxy = {0, 0, 1}
* xyx = {0, 1, 0}
*/
template <int l>
Mat cartesian_gaussian_rotation_matrix(const Mat3x3 rotation) {
constexpr int num_moments = (l + 1) * (l + 2) / 2;
Mat result = Mat::Zero(num_moments, num_moments);
auto cg_powers = power_index_arrays<l>();
int p1_idx = 0;
for (const auto &p1: cg_powers) {
int p2_idx = 0;
// copy as we're permuting p2
for (auto p2: cg_powers) {
do {
double tmp{1.0};
for (int k = 0; k < l; k++) {
tmp *= rotation(p2[k], p1[k]);
}
result(p2_idx, p1_idx) += tmp;
} while (std::next_permutation(p2.begin(), p2.end()));
p2_idx++;
}
p1_idx++;
}
return result;
}
In essence this is just accounting for all the redundant terms we'd get that arise from the permutations of tensor indices. Going through the code line-by-line
So our rotated column vector (or block of a matrix) of molecular orbital coefficients will simply be
Rotating pure Gaussiansโ
If you've made it this far, well now we're in for the more painful bit.
Obviously we can't directly use the above rotation matrix for a set of pure Gaussian molecular orbital coefficients, but we can convert the coefficients to their equivalent Cartesian representation, perform the rotation, then convert back.
Frankly, I find this alternative preferable to dealing with Wigner D matrices, Clebsch-Gordan coefficients or any of the other painful aspects of rotating spherical harmonics...
Spherical to Cartesian transformations
Thankfully, I'm far from the first person to have wanted to perform this conversion, so there are tabulated transforms (along with implementations to generate them) available e.g., in HORTON. Much of this is based on the Schlegel & Frisch paper from the mid 90s.
This gives us , a transformation matrix which will convert from Cartesian coefficients to pure coefficients. But what of the inverse transform? is not an invertible matrix (obviously, as it's not square). Thankfully this is mentioned in the Schlegel & Frisch paper, and we can find the inverse transformation via:
Here is the matrix of overlap between normalized Cartesian gaussians of the same total angular momentum, where (if are powers of respectively) can be derived from the following horrible expression:
where all of , and are even, otherwise.
The neat part after all of this, is that we can simply transform , the Cartesian Gaussian rotation matrix into , the pure/spherical Gaussian rotation matrix as follows:
and perform the rotation directly on the spherical molecular orbital coefficients:
Whew... Hopefully I don't have any errors in the equations, but if I do feel free to name and shame me for being so sloppy.