One post tagged with "linear algebra"

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, $\mathbf{O}$, 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 $K$ is the degree of contraction (i.e. number of primitive Gaussian functions), $N$ is the normalization constant we're going to ignore for now, $\alpha_k$ is the exponent for each primitive basis function and $c_k$ 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 $\alpha_k$, and the right column being $c_k$

H     0 
S   3   1.00
      3.42525091             0.15432897       
      0.62391373             0.53532814       
      0.16885540             0.44463454   

With the exception of the $f(\mathbf{r} - \mathbf{r_k})$ term, $\chi$ is purely radial. However, there are two common choices of $f$: Cartesian and spherical.

Cartesian polynomials​

In a Cartesian Gaussian basis set, $f$ takes the following form based on angular momentum $l$:

where $i + j + k = l$ and all of $i,j,k$ are whole numbers.

Because multiplication is commutative i.e. $x^2 y = y x^2$, for each angular momenta, $l$, there are a number of 'unique' functions e.g.

and so on...

Spherical/Pure polynomials​

In a pure (or spherical) Gaussian basis set, $f$ usually corresponds to the real regular solid harmonics, which are directly related to the usual spherical harmonics $Y_l^m$ These are a bit more of a rabbit hole..

Spherical harmonics ($Y_l^m$)

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 $N$ is yet another different normalization factor, $\theta$ and $\phi$ are the polar angle/colatitude $[0, \pi]$ and the azimuthal angle/longitude $[0, 2 \pi]$ respectively, and $P_l^m$ is an associated Legendre Polynomial

Regular solid harmonics ($R_l^m$)

Given the previous definition of the spherical harmonics, the regular solid harmonics take the form:

or, simplified:

Real regular solid harmonics ($R_{l,m}$)

Finally, the functions $f$ that are actually used are of the form:

where

This means that there are $2 l + 1$ components for a pure Gaussian basis with angular momentum $l$ e.g.:

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 $l$ rank tensor e.g. a $3 \times 3$ matrix for a $d$ function, a $3 \times 3 \times 3$ cube for an $f$ 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 $\mathbf{c}$, a transformation matrix which will convert from Cartesian coefficients to pure coefficients. But what of the inverse transform? $\mathbf{c}$ 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 $\mathbf{S}$ is the matrix of overlap between normalized Cartesian gaussians of the same total angular momentum, where (if $i,j,k$ are powers of $x,y,z$ respectively) $\mathbf{S}$ can be derived from the following horrible expression:

where all of $i_1 + i_2$, $j_1 + j_2$ and $k_1 + k_2$ are even, $0$ otherwise.

The neat part after all of this, is that we can simply transform $\mathbf{R}$, the Cartesian Gaussian rotation matrix into $\mathbf{R}'$, 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.