I’ve been thinking lately about the transformations between the cartesian and spherical basis. For instance, in the SAD guess one is supposed to symmetrize the densities of each atom. For a spherical basis this is very simple: you just set the [(lm),(l’m’)] components to zero if (lm)!=(l’m’), and for the diagonal you use an average over the m values. For a cartesian basis it is not nearly as simple… However, if you assume that the contaminants can be forgotten (e.g. the 1 S spherical function on the cartesian D manifold), then you should be able to do the spherical average for cartesian basis sets just by doing cartesian -> spherical -> manipulate -> back to cartesian.

However, the back transformation doesn’t appear to be trivial. If I take a matrix in the AO basis, M_AO, and transform it to the SO basis with

M_SO = aotoso^T * M_AO * aotoso

and then back-transform it as

M_AO, new = sotoao^T * M_SO * sotoao

I should be getting back something close to the original matrix. Say, if the original M_AO is a unit matrix, then I should be seeing something close to a unit matrix, with just the spherical contaminants removed.

However, this appears not to be the case - the resulting matrix contains some significantly larger values as well as off-diagonal elements - and the same thing also happens if I start with the SO matrix.

What am I doing wrong? I’m running in C1 symmetry and I’ve used the aotoso and sotoao functions in psi4.core.Wavefunction.