Can a spin-flip TD-DFT singlet ground state wavefunction be written as a sum over MOs with occupancy information?

Question

Using collinear spin-flip(SF) TD-DFT with the Tamm-Dancoff approximation, BHanHLYP and a simple test molecule (ethene, with a singlet ground state).

It's my understanding, using ethene as an example, that spin-flip DFT starts with a unrestricted Kohn-Sham single-reference triplet reference state, then generates the other TD-DFT 'excited' states using the MOs and occupancies from that reference state and a calculated excitation vector. The relative energy of the ground state can be calculated, and, at least in ORCA, the molecular geometry can be optimized on this ground-state PES. So far, so good.

My question is: can the wavefunction of the resulting singlet SF ground state, which appears as an effective negative-energy 'excited' state, have its electron density distribution represented (and written in machine-readable form), as another set of only MOs and occupancies?

It doesn't matter if this output set is single-reference or multi-reference, for my purposes, just so long as it can be output as only MOs and occupations.

If this is implementation-dependent, which codes implementing spin-flip TD-DFT can do this? Does the Tamm-Dancoff approximation need to be dropped in order to achieve this?

Answer 1

What you are looking for are the natural orbitals of the spin-flip ground state wavefunction. As a matter of fact, the one-electron density matrix from any electronic wavefunction, no matter whether it is multireferential or whether it is an excited state wavefunction, can be written as a linear combination of contributions of natural orbitals (NOs):

$$ D = \sum_i n_i |\psi_i\rangle \langle\psi_i| \tag{1} $$

Given the one-electron density matrix in an orthogonal basis, one can find the NOs and the corresponding occupation numbers by diagonalizing the density matrix. If the density matrix is only available in a non-orthogonal basis (such as the AO basis), one can transform it to an orthogonal AO basis (via either Schmidt orthogonalization or symmetric orthogonalization. EDIT: other more robust orthogonalization methods also exist, see Susi's comment below), diagonalize the density matrix, and transform the NO coefficients back to the AO basis. The density matrix does not contain all the information in the wavefunction, but it contains sufficient information for determining the electronic density.

Any decent program that supports SF-TDDFT or SF-TDA can generate the natural transition orbitals (NTOs) or difference density matrices of the spin-flip states (both the ground state and the excited states), or both; the former can be transformed to the latter by plugging the NTOs and their occupation numbers into Eq. (1). To get the NOs of a spin-flip state, one must add the difference density matrix to the density matrix of the reference state to get the density matrix of the spin-flip state. Diagonalizing this density matrix yields the NOs and their occupation numbers. Depending on the computational chemistry software that you use, there is a good chance that some or all of these post-processing steps can be done by Multiwfn.