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())
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:)
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
Thank you for pointing that out.
I ran both
Da = np.asarray(ccsdt_wfn.Da())
Da_AO = np.asarray(ccsdt_wfn.Da_subset("AO"))
Da == Da_AO
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.