Obtaining excited state 1RDMs and 1TDMs in AO Basis for symmetries other than C1

1

Hello,
I am new to Psi4, and I really like the Python API! I would like to develop an interface with Psi4 for a software package I am developing which is aimed at describing metastable states using the complex absorbing potential method. We use a subspace projection scheme in which the CAP is projected onto a set of eigenvectors obtained from an excited electronic structure method such as CI,EOM-CC etc. My package computes the CAP matrix in atomic orbital basis, and then uses 1RDMs and TDMs (passed to my program from Psi4 as numpy matrices) to compute CAP matrix elements (a 1-electron operator) in the state basis. So far, I’ve gotten it to work with CI in Psi4, but only for C1 symmetry:

import psi4
import numpy as np
mol = psi4.geometry("""H 0.0000000000 0.0000000000 0.3705000000
H 0.0000000000 0.0000000000 -0.3705000000
Symmetry C1""")
E, wfn = psi4.energy('scf/aug-cc-pvqz', return_wfn=True)
mints = psi4.core.MintsHelper(wfn.basisset())
S_mat = np.asarray(mints.ao_overlap())
# verify overlap matrix in my program to check AO ordering, this always works
# s.check_overlap_mat(S_mat,"psi4")
nstates= 10
psi4.set_options({"opdm":True,"num_roots":nstates,"tdm":True,"basis":"aug-cc-pvqz"})
ci_energy, ci_wfn = psi4.energy('FCI', return_wfn=True)
mo_coeff = ci_wfn.Ca()
for i in range(0,nstates):
    for j in range(i,nstates):
        opdm_mo = ci_wfn.get_opdm(i, j, "SUM", True)
        opdm_ao = psi4.core.triplet(mo_coeff, opdm_mo, mo_coeff, False, False, True)
        opdm_ao = opdm_ao.to_array(dense=True)
        # pass tdm/rdm for each state/transition to my program, only works properly for C1
        # pc.add_tdm(opdm_ao,i,j,"psi4")

My program currently understands the GTO ordering that Psi4 uses (verified by checking overlap matrix), so when there’s just a single irrep, things work fine. But I haven’t been able to figure out how to perform this MO–>AO transformation (in the original AO ordering) when there are multiple irreps. I know the dense=True argument to the to_array function at returns a single matrix rather than a matrix for each irrep, but it seems to be sorted by irrep, and I’m not sure how to transform this into the original AO ordering, which is required for my program. So to summarize, I have two main questions:

  1. Is there a convenient way to perform MO–>AO transformation for 1RDMs and 1TDMs which recovers the original atomic orbital ordering (i.e. which matches the orbital ordering used in the overlap matrix) when there are multiple irreps?
  2. Are there functions similar to get_opdm for CI which provide direct access to 1RDMs and 1TDMs between excited states for other methods such as EOM-CC?

Any help would be greatly appreciated!

Here are the topic threads I have consulted:

I’ve also consulted the doc string in the Molden function here: psi4/driver.py at a548d758c52fd76d5f10e9592e36d84b127bb432 · psi4/psi4 · GitHub

2

I did not find anything in the respective code base.

In the case of symmetry you have an SO instead of AO basis.
Likely you need to to do MO->SO->AO.
The molden function has an SO->MO example, you’d need to reverse that and then you could try:
https://psicode.org/psi4manual/master/api/psi4.core.Matrix.html#psi4.core.Matrix.remove_symmetry

Do you need all the 1e and 2e integrals from mints = psi4.core.MintsHelper(wfn.basisset())? For the overlap you can also use wfn.S()

1 Like
3

Wonderful, remove_symmetry did the trick! The only other piece needed was grabbing the so2ao matrix from mints. The snippet now looks like this:

so2ao = mints.petite_list().sotoao()
for i in range(0,nstates):
    for j in range(i,nstates):
        opdm_mo = ci_wfn.get_opdm(i, j, "SUM", True)
        opdm_so = psi4.core.triplet(ci_wfn.Ca(), opdm_mo, ci_wfn.Ca(), False, False, True)
        opdm_ao = psi4.core.Matrix(nbasis,nbasis)
        opdm_ao.remove_symmetry(opdm_so,so2ao)

In regards to get_opdm, in case anybody was curious, ADCC can do it with state2state.transition_dm, but I couldn’t find it for anything else. I’ll stay tuned for if/when something similar gets implemented for other methods :slight_smile:

Thanks for the help!

4

Great that you figured it out!

If you need it for a specific method it is best to open a github issue with a feature request so that the respective module developers are aware of what users need!

5

Thanks for such a thorough report :slight_smile: I hadn’t seen Holger’s reply so I took a very slightly different approach that I’ll post here in case it offers a useful alternative for anybody following this thread. The wavefunction can directly return the MO coefficients in the AO (unsymmetrized) basis via a helper function. Under the hood, it’s simply applying the AO2SO matrix for you. I just took this Ca matrix and cast it to Numpy, doing the same for the density matrices. The transformations can then be handled by Numpy itself, albeit inefficiently because the symmetry is removed at this point.

import psi4
import numpy as np
mol = psi4.geometry("""H 0.0000000000 0.0000000000 0.3705000000
H 0.0000000000 0.0000000000 -0.3705000000""")
E, wfn = psi4.energy('scf/aug-cc-pvqz', return_wfn=True)
mints = psi4.core.MintsHelper(wfn.basisset())
S_mat = np.asarray(mints.ao_overlap())
# verify overlap matrix in my program to check AO ordering, this always works
# s.check_overlap_mat(S_mat,"psi4")
nstates= 10
psi4.set_options({"opdm":True,"num_roots":nstates,"tdm":True,"basis":"aug-cc-pvqz"})
ci_energy, ci_wfn = psi4.energy('FCI', return_wfn=True)
mo_coeff = ci_wfn.Ca_subset('AO', 'ALL').to_array(dense=True)
for i in range(0,nstates):
    for j in range(i,nstates):
        opdm_mo = ci_wfn.get_opdm(i, j, "SUM", True).to_array(dense=True)
        opdm_ao = np.einsum('mi,ij,nj->mn', mo_coeff, opdm_mo, mo_coeff, optimize=True)

I checked that the above code runs, but didn’t check the results, so proceed with caution if you choose to use any of it!

6

This topic was automatically closed 60 days after the last reply. New replies are no longer allowed.