I’ve been digging around in the TDSCF code as well as various bits of the nwchem source, and also various literature sources to try to get my head around what is what under the different conditions. I can reproduce the excitation coefficients that I was initially seeking fairly easily by appending a chunk of code to an input file but it would be lot cleaner if it were within PSI4 itself (and symmetry, etc. would become a lot easier to deal with).

Before I attempt to back-engineer what I have into the TDSCF code, I though it would be wise to check that I’m doing this the “correct” way because it seems to be very easy to get similar values using a variety of approaches!

My understanding is this:

The left and right eigenvectors are left = |Xn - Yn> and right = |Xn + Yn>, where Xn contains the excitation coefficients and Yn contains the de-excitation coefficients.

The left and right eigenvectors are both nocc rows by nvirt columns.

If TDA is used, de-excitations are ignored so Yn are zero, and so the left and right eigenvectors are identical: left = right = |Xn>. In this case, elements of either eigenvector are the excitation coefficients that I am after. Typically, these are normalised so that the sum of the squares of the coefficients = 0.5 (presumably, so that the alpha + beta coefficients are together normalised to 1).

If TDA isn’t used, then (I think!) I can simply use X = 1/2(right + left) and Y = 1/2(right - left), and then use X for the excitation coefficients as above and Y for the (usually smaller) de-excitation coefficients. Typically, these are normalised so that the (sum of the squares of the excitation coefficients) - (sum of the squares of the de-excitation coefficients) = 0.5.

(Normalisation of these coefficients is not always well documents and may not be performed by some software.)

Assuming the above is correct, it will be fairly easy for me to add the relevant bits to scf_response.py. I *think* that the TDSCF code currently only handles restricted formalism and singlets for TD-DFT and so I wouldn’t have to worry about alpha and beta coefficients, at least initially (although I think triplets are currently possible for TD-HF so I should probably explore that, too!).

I’d also add a new parameter (TDSCF_COEFF?) as a threshold for which coefficients to print in the output.

A couple of (fairly) basic questions that I can’t readily find an answers to that might improve the output I’m generating…

Is there a simple way to obtain MO symmetries from a wavefunction object? (As printed by HF::print_orbitals() in libscf_solver/hf.cc: I can’t see if the labels make it to the python interface or not.)

Is there a direct way to get the number of virtual orbitals or do I need to use something like: wfn.nmo()-wfn.nalpha()? (Although, I can easily get these from the dimensions of the left or right eigenvectors.)