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