Hi all,

I’m attempting to reproduce the the process described here [https://doi.org/10.1016/S0009-2614(01)00828-4] by Mayer involving Mayer Bond Indices. If I understand the paper correctly (which is questionable), I am attempting to 1) calculate the orbitals/density matrix at scf, 2) drop an electron from the HOMO, 3) calculate a one-electron property without recalculating the orbitals. Is this possible with PSI4?

Thanks for any assistance

The bad news is that the code can’t do this right now. The good news is that it’s very easy to write the code in an input. Here’s my hastily hacked together example:

```
import psi4
import numpy as np
psi4.geometry("""
O 0.0 0.0 0.0
H 1.0 0.0 0.0
H 0.0 1.0 0.0
symmetry c1
""")
psi4.set_options({'basis': 'cc-pvdz',
'reference' : 'uhf'})
e, wfn = psi4.energy('scf', return_wfn=True)
nalpha = wfn.nalpha()
oeprop = psi4.core.OEProp(wfn)
oeprop.add('MAYER_INDICES')
# Compute from native psi4 quantities
oeprop.compute()
# Recompute the beta density, and rerun as a sanity check
Ca = wfn.Ca().np
nbeta = nalpha
wfn.Db().np[:] = np.einsum('pi,qi->pq', Ca[:,:nbeta], Ca[:,:nbeta])
oeprop.compute()
# Same again, but remove a beta electron
nbeta = nalpha-1
wfn.Db().np[:] = np.einsum('pi,qi->pq', Ca[:,:nbeta], Ca[:,:nbeta])
oeprop.compute()
```

This works by breaking into the wavefunction’s beta density memory and making changes. The second call is just to verify that I computed the density matrix correctly. Good luck!

Works great! Thanks Andy.