Atom indexes of density matrix

Hi psi4 developers & community,

I was wondering if there is any way to identify the corresponding atoms (and their orbitals) from a molecule that map to the rows and columns of the electron density matrix extracted from the wfn object after running scf_e, scf_wfn ="scf", return_wfn=True) [“basis”: “sto-3g”].

For an example with a glycine molecule (10 atoms, 30 atomic orbitals).
Da = np.asarray(scf_wfn.Da())

(30, 30)

which matches the total number of atomic orbitals contributed by Cs, Hs, Ns, & Os for a glycine molecule.

But I was looking for a way to map the rows and columns of Da to the atomic orbitals which contributed to them.

Perhaps something like this as an example.

Does anyone know if this is possible within psi4’s capabilities?

I appreciate any help I can get on this:)

1 Like

First, you need to realize that Da uses symmetry-adapted AOs rather than AOs. Use Da_subset("AO") rather than Da to get atomic orbitals.

AOs are ordered as follows: loop over atoms, loop over quantum number l (angular), loop over quantum number n (principal quantum), loop over those by m_l in pattern 0, +1, -1, +2, -2, …

A quick but non-programmatic way to match an AO index to a physical AO is to use set print_mos true and look for the “Molecular Orbitals” section of the output file. The row labels will contain the information you’re looking for. For example, minimal basis water:

 1    O1 s0         0.9967095    0.2233723    0.0000000    0.0932516    0.0000000
 2    O1 s0         0.0153898   -0.8532247    0.0000000   -0.5183826    0.0000000
 3    O1 p0         0.0029549   -0.1154667    0.0000000    0.7687640    0.0000000
 4    O1 p+1        0.0000000    0.0000000    0.0000000    0.0000000    1.0000000
 5    O1 p-1        0.0000000    0.0000000    0.6088096    0.0000000    0.0000000
 6    H2 s0        -0.0034491   -0.1510180   -0.4451103    0.2948036    0.0000000
 7    H3 s0        -0.0034491   -0.1510180    0.4451103    0.2948036    0.0000000
1 Like

Thank you for pointing that out.
I ran both wfn.Da() and wfn.Da_subset("AO") as:
Da = np.asarray(ccsdt_wfn.Da())
Da_AO = np.asarray(ccsdt_wfn.Da_subset("AO"))
Da == Da_AO
Return True for all elements.
I think the amino acids I am working with do not have symmetry and these two are hence equal.

Thank you for the suggestion to

loop over atoms, loop over quantum number l (angular), loop over quantum number n (principal quantum)

Do you know if there is any way to pipe the output from the scf run into Python itself as a dict or some other psi4 internal object instead of writing it to disk as a file? It would be much simpler to access something in Python directly rather than having to parse a file.