Atom indexes of density matrix

1

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 = psi4.energy("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())
Da.shape

(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.

image
image1021×370 91 KB

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
2

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
3

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"))
and
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.