Figure out ordering conventions of basis functions in other codes

[tovrstra]

We should know of other codes (and their file formats) how they order segmented SP basis functions (in case they are not working with generalized contractions). See #39 for an example. PySCF reorders segmented SP basis functions by angular momentum after: All S functions are put first, followed by all P functions. How is this working in other codes: PSI4, ORCA, Gaussian, MOLPRO, NWChem, ...

In the case of Gaussian, generalized contractions are not segmented, i.e. functions of an SP shell are not reordered because they form one generalized shell, which saves some computing time because the integrals over S and P primitives are computed in one pass.

This may also affect wavefunction files written by these codes, e.g. Molden and Molekel files, which is our primary means to test this: compute a HF wavefunction of water with a 3-21G basis and write out a molden or mkl file. What is the ordering of the coefficients? Second way to test this is by writing out integrals. MOLPRO, PySCF, Gaussian (to some extent) and PSI4 can do this.

[linminhtoo]

hello @tovrstra , it's been many years since this issue was made, but I was wondering if you could share any of your insights from this. I'm building a machine learning model to predict things in the atomic orbital basis and I need to ascertain the ordering of these basis functions in programmes (specifically psi4). any pointers would be greatly appreciated :)

[tovrstra]

Sorry, I'm not too actively following up the issues here. The newer IOData has an improved API to represent the ordering of the orbitals when loading a wavefunction from disk. Many people have contributed loaders for different wavefunction formats, and both Molden anf FCHK files generated with recent Psi4 versions should work well.

A small comment on your research: keep in mind that the expansion coefficients (just like e.g. Mulliken charges) may become ill-conditioned in the complete basis set limit. Predicting these with machine learning would amount to predicting practically random numbers (with only a few well-defined linear combinations). This means you have to make careful choices in terms of basis set.

[linminhtoo]

Sorry, I'm not too actively following up the issues here. The newer IOData has an improved API to represent the ordering of the orbitals when loading a wavefunction from disk. Many people have contributed loaders for different wavefunction formats, and both Molden anf FCHK files generated with recent Psi4 versions should work well.

A small comment on your research: keep in mind that the expansion coefficients (just like e.g. Mulliken charges) may become ill-conditioned in the complete basis set limit. Predicting these with machine learning would amount to predicting practically random numbers (with only a few well-defined linear combinations). This means you have to make careful choices in terms of basis set.

Hey, thanks so much for your reply! I'll have a look at those.

Regarding your 2nd point, could you elaborate a bit more, or perhaps share any references further explaining this phenomenon? I am trying to predict the converged molecular coefficient matrix at the cc-pVDZ basis set, would that fall into this problem?

[tovrstra]

I would guess cc-pVDZ is on the safe side. Things will become more problematic with diffuse functions, e.g. Aug-cc-pVDZ.

There are quite a few papers showing the deficiencies of Mulliken charges, which suggest the same problem arises with the underlying orbital coefficients. Just to name an old paper of mine: https://biblio.ugent.be/publication/1997717 In figure 7, Mulliken charges of a dipeptide vary rather wildly as a function of a torsional scan. Given that other methods show much more modest trends, these big swings are a bit suspect.

In the end, for your work, the condition number of the overlap matrix is the key parameter to check. When that condition number becomes large, the orbital coefficient will become very sensitive to small changes in anything (geometry, level of theory, integration grids, ...), simply because the basis set is degenerate. If the coefficients are not robust, predicting them will become hard for just numerical reasons. One way to circumvent this issue, is to train your ML model that for wavefunctions using robust observables as training data. That way, it is ok not to predict flaky expansion coefficients exactly, as long as the observables come out right.

[PaulWAyers]

@tovrstra has a good point: the overlap matrix will tell you whether the MO coefficients are extremely sensitive to small changes.

Note, @linminhtoo , that the problem you're describing is already ill-defined, because the MO coefficients change drastically (but immaterially) when the molecule is rotated in space. For example, the key contribution to $\sigma$ bonds in a diatomic will shift from the $p_z$ basis functions to the $p_x$ basis functions when the molecule is alligned along the $x$ axis instead of the $z$ axis. So I'm not even sure that orbital ordering is relevant, honestly, because changes in rotation effectively convert different orders into each other.

One could argue that the issue is even worse than this in Hartree-Fock and DFT, because unitary transformations of the occupied molecular orbitals do not change any observable property.

There are tricks to try to fix this (trying to make MO coefficients as sparse as possible) but I would endorse @tovrstra 's advice: use observables. I can't actually think of any case where MOs are directly useful where there isn't an observable, derived from te MOs, that achieves the same objectives.