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.
Hi,
More than a year later, but do you have any answer to your own question now ?
Actually I’m hitting the same issue: I want to get MO coefficients in the AO basis but starting from a SCF calculation in the SO one. So I also used the sotoao matrix from psi4.core.PetiteList.
I do M_AO = sotoao * M_SO * sotoao.T
and I compare this M_AO with the one obtained by doing the SCF calculation with puream false.
Both are clearly not close to each other.
Hi, I’m one of the developers.
I strongly suspect this is a bug in Psi. sotoao should be an orthogonal matrix, but it is not in some of the cases I’ve tried. I need this fixed for something I’m working on, so I’ll investigate and report back to you.
The matrices shouldn’t match unless your basis is limited to s and p functions.
Can you elaborate? Are you saying Psi’s M_AO shouldn’t be M_AO = sotoao * M_SO * sotoao.T as bsen was talking about that? That sotoao shouldn’t be orthogonal, as I was talking about?
If you have d, f, etc shells in the basis, then there are fewer spherical functions than cartesian ones. The calculation with puream false will have more variational freedom than the one with puream true, and thus will give different AO coefficients.
Yes I agree if you have d or f functions. So if we have d or f functions, can we still use the sotoao transformation or would it just be plain wrong ?
Another issue that I would like to talk about is the spherical averaging. If you take just the scf calculation of an oxygen or carbon atom for instance, you will not have the same MO coefficients wrt the p_x, p_y and p_z functions. Of course I could do this averaging by hand but I’d like something general, maybe there is an option for this in psi4 ?
Coming back to the original problem, I enclose a little script using Psi4.Numpy on the oxygen atom to check the sotoao transformation. Do you agree that in this case (sto-3g), the coefficient matrices should be the same ?
import psi4
import sys
import numpy as np
psi4.set_options({'basis': 'sto-3g',
'reference': 'rohf',
'puream':'true'})
oxygen = psi4.geometry("""
O
symmetry c1
""")
scf_e, scf_wfn = psi4.energy('scf', return_wfn=True)
coeff_a_spherical = scf_wfn.Ca().np
sotoao = scf_wfn.sobasisset().petite_list().sotoao().np
psi4.set_options({'basis': 'sto-3g',
'reference': 'rohf',
'puream':'false'})
scf_e, scf_wfn = psi4.energy('scf', return_wfn=True)
coeff_a = scf_wfn.Ca().np
coeff_a_cartesian = np.dot(np.dot(sotoao,coeff_a_spherical),sotoao.T)
print("Are arrays equal ?", np.allclose(coeff_a, coeff_a_cartesian))
I want to go one thing at a time, since there are lots of small details to get hung up over regarding your example:
Oxygen is a poor example because even after you’ve chosen your AOs, the MO coefficients are not uniquely determined. You can rotate the occupied p orbitals however you wish and get the same result. Let’s change our test case to water, which doesn’t have this problem. As a nice bonus, we can also use RHF rather than ROHF. To eliminate numerical noise, let’s change our SCF to use exact integrals and increase the convergence threshholds by a few orders of magnitude.
With these modifications, I get the following input file.
import psi4
import sys
import numpy as np
psi4.set_options({'basis': 'sto-3g',
'puream':'true',
'e_convergence': 12,
'scf_type': 'pk'})
oxygen = psi4.geometry("""
O
H 1 1.0
H 1 1.0 2 104.5
symmetry c1
""")
scf_e, scf_wfn = psi4.energy('scf', return_wfn=True)
coeff_a_spherical = scf_wfn.Ca().np
sotoao = scf_wfn.sobasisset().petite_list().sotoao().np
psi4.set_options({'basis': 'sto-3g',
'puream':'false'})
scf_e, scf_wfn = psi4.energy('scf', return_wfn=True)
coeff_a = scf_wfn.Ca().np
coeff_a_cartesian = np.dot(np.dot(sotoao,coeff_a_spherical),sotoao.T)
print("Are arrays equal ?", np.allclose(coeff_a, coeff_a_cartesian))
This has indeed demonstrated a problem in the AO -> SO conversion. However, the problem is not in Psi4 but in your conversion formula.
The coefficient matrix is (ao x occupied orbitals). We have no need to transform the second index, only the first. We want to go from ao to so, so the actual transformation is np.dot(sotoao.T, coeff_a_spherical). When you update coeff_a_cartesian accordingly, the test passes.
Aw, yes indeed thanks for the correction. It works also on my side now.
But now as @susilethola said, we cannot do this transformation if d of f functions are there (so just doing the same example as before with the cc-pVDZ basis will lead to an error).
Is it then impossible to transform MO coefficients express in terms of spherical AO to the cartesian one, and vice-versa ?
I think so. The transformation is not one-to-one; you lose the contaminants.
It would be different if the transformation was unitary, i.e. it would keep the S function on the D shell and the P functions on the F shell etc…
1 Like
It’s not obvious to me how you’d do such a transform. You could maybe define some pseudoinverse? I’ll keep thinking about it and let you know if I have anything.
At the very least, there needs to be a warning about this problem…
2 Likes
Status update: I’ve talked with the other developers. We were already updating our integrals package from libint1 to libint2. As part of the transition, we’ll be changing the PetiteList so it will always transform either spherical to spherical or cartesian to cartesian. Then ao2so and so2ao will work without any problems at all.
This is a high priority for us, but a change this big will take time, so I don’t have an estimate on when it will be fixed.
2 Likes