Clarifying TDDFT Output

sbembenek · October 18, 2021, 6:37am

Looking to clarify some of the output from the TDDFT calculation.

For my input:

molecule mol {
0 1
O 0.00000000 0.00000000 0.00000000
H 0.00000000 -0.79069000 0.61221700
H 0.00000000 0.79069000 0.61221700

symmetry c1
}

set scf_type df
set basis aug-cc-pVTZ
set reference rks
set s_tolerance 1e-9
set_num_threads(62)

set save_jk true

import numpy as np
import psi4
from psi4.driver.procrouting.response.scf_response import tdscf_excitations
from psi4.driver.p4util import spectrum

e, wfn = energy('b3lyp-d3bj', return_wfn=True)
res = tdscf_excitations(wfn, states=20)

Question #1: I added ‘symmetry c1’ per an example I saw – is this the standard trick to get the breakdown of the excitations & de-excitations for a symmetric molecule?

So, I did get the breakdown of the excitations & de-excitations. For example:

Excited State    1 (1 A):   0.24630 au   185.00 nm f = 0.0434
  Sums of squares: Xssq =  1.000490e+00; Yssq =  4.902097e-04; Xssq - Yssq =  1.000000e+00
     5 ->  6    0.991925 (98.391%)
     5 ->  8    0.117839 ( 1.389%)

I’m reading this as, for Excited State 1 (with wavelength 185nm and oscillator strength of 0.0434) you have two contributions:

A dominant transition (i.e., 98.391%) from MOs 5–>6
A ‘minor’ (1.389%) transition from MOs 5–>8

Question #2: Is this correct?

Question #3: Can you clarify ‘Sums of squares’ to me?

Question #4: In addition to the oscillator strengths, can the respective transition dipole moment also be printed? I understand their relationship – just would like to see them both at the same time.

laz · October 18, 2021, 7:06pm

It was me who pulled together the code that outputs a lot of the TDDFT output that you are asking about, mostly because I discovered the information I wanted was there but not in a very nice format.

Note that I am no way an expert in this but I worked it out from reading the TDSCF page, a bit of literature, and then comparing what I got with the output from other software. The following is my understanding of the theory, which may be a bit off.

(I was hoping that the documentation would pull most things from the source code comments but it seems not. I need to add to it! There are a couple of undocumented options here that I mention below.)

I think you have two options to run the TDSCF / TDDFT code. Either calling energy() followed by tdscf_excitations() (as you are), or you can simply use something like energy(‘td-b3lyp’) (see the top of the manual page).

Q1: yes: this code only currently works if you force C1 symmetry. At the time, I couldn’t work out how to extract the relevant information from the different submatrices corresponding to the different representations when higher symmetry is used. It would be nice to “mend” that but I don’t have the time at present.

Q2: pretty much that. By default, only coefficients >0.1 are printed. If you set TDSCF_COEFF_CUTOFF to <0.1 then you will get all coefficients larger than the cutoff printed. (Alternatively, I think you can add something like “coeff_cutoff = 0.01” to a call to tdscf_excitations() to print coefficients larger than 0.01.) Note that the output becomes very big very quickly if you lower this by much.

Q3. The X matrix contains excitations for the excited state and the Y matrix contains de-excitations (or nothing if using TDA). If you have a large amounts of de-excitations (large coefficients in the Y matrix) then I think it means that the ground state isn’t very well represented by that level of theory.

The matrices are normalised so that Xssq - Yssq = 1.0. Hence, Xssq will generally be approx 1 and Yssq will generally be a lot smaller. (I’m not sure if these values are output by other software or whether I simply added it for my own piece of mind when generating X and Y from the left and right eigenvectors.)

Q4. Various types of TDM vectors are calculated but not usually output. They are assigned to variables with pretty horrendous names. I added another option variable: TDSCF_TDM_PRINT (or call tdscf_excitations() with a tdm_print parameter). This is a list variable that can contain one or more of “E_TDM_LEN”, “E_TDM_VEL”, and “M_TDM” to print out Electric Transition Dipole Moments (length), Electric Transition Dipole Moments (velocity), and Magnetic Transition Dipole Moments, respectively. By default none are printed.

Caveat: it was a while back that I did this and some of the above is from memory.

sbembenek · October 18, 2021, 10:14pm

Thanks! this help a lot.

I tried a few formats on trying to get to the TDM info, but was unsuccessful.

Can you give me an explicit example of how I would get the (for example) E_TDM_LEN”, “E_TDM_VEL to print out. And does it print out as a column along with →


                                    Excitation Energy         Total Energy        Oscillator Strength             Rotatory Strength       
     #   Sym: GS->ES (Trans)        au              eV              au          au (length)    au (velocity)    au (length)    au (velocity)

laz · October 19, 2021, 5:32pm

This is with current dev version (through conda) but I don’t think the TDSCF code has changed since the 1.4 release. Here’s a basic example for benzene:

memory 1 gb
set_num_threads(4)

molecule benzene {
    C         -1.394814023431    -0.000000000000     0.000000000000
    C         -0.697441089778     1.207970230209     0.000000000000 
    C         -0.697441089778    -1.207970230209     0.000000000000
    C          0.697441089778     1.207970230209     0.000000000000 
    C          0.697441089778    -1.207970230209     0.000000000000  
    C          1.394814023431     0.000000000000     0.000000000000 
    H         -2.480916576762    -0.000000000000     0.000000000000 
    H         -1.240492371820     2.148551413522    -0.000000000000 
    H         -1.240492371820    -2.148551413522     0.000000000000 
    H          1.240492371820     2.148551413522    -0.000000000000  
    H          1.240492371820    -2.148551413522    -0.000000000000 
    H          2.480916576762     0.000000000000     0.000000000000 
symmetry c1
}

set {
  basis cc-pVDZ
  tdscf_states 10
  tdscf_coeff_cutoff 0.01
  tdscf_tdm_print ['E_TDM_LEN', 'E_TDM_VEL']
}

energy('td-b3lyp/cc-pvdz')

This prints excitations and de-excitations larger than 0.01 and also outputs tables of TDM vectors in length and velocity representations. E.g. (just after the table of transitions):

Electric Transition Dipole Moments (Length) (au):
State      X          Y          Z
    1  0.000000   0.000011  -0.000000 
    2  0.000745  -0.000000   0.000000 
    3  1.793550   0.000000   0.000000 
    4 -0.000000   1.793646   0.000000 
    5 -0.000000  -0.000000   0.000000 
    6 -0.000000   0.000000  -0.000000 
    7 -0.000000   0.000000   0.000000 
    8 -0.000000   0.000000   0.000000 
    9 -0.000000  -0.000001  -0.002783 
   10 -0.000000  -0.000000   0.206327 

Electric Transition Dipole Moments (Velocity) (au):
State      X          Y          Z
    1 -0.000000  -0.000001   0.000000 
    2 -0.000173   0.000000  -0.000000 
    3 -0.478996  -0.000000   0.000000 
    4  0.000000  -0.479028   0.000000 
    5 -0.000000   0.000000  -0.000000 
    6  0.000000  -0.000000   0.000000 
    7 -0.000000  -0.000000  -0.000000 
    8  0.000000   0.000000  -0.000000 
    9  0.000000  -0.000000   0.000738 
   10  0.000000   0.000000  -0.054969 


Contributing excitations and de-excitations
Only contributions with coefficients > 1.00e-02 will be printed:

Excited State    1 (1 A):   0.20205 au   225.50 nm f = 0.0000
  Sums of squares: Xssq =  1.000886e+00; Yssq =  8.863692e-04; Xssq - Yssq =  1.000000e+00
    10 -> 33    0.010751 ( 0.012%)
    11 -> 34   -0.010750 ( 0.012%)
    14 -> 24   -0.017165 ( 0.029%)
    15 -> 31   -0.010687 ( 0.011%)
    16 -> 32   -0.010684 ( 0.011%)
    20 -> 22   -0.706812 (49.958%)
    21 -> 23   -0.706784 (49.954%)
    20 <- 22   -0.011917 ( 0.014%)
    21 <- 23   -0.011922 ( 0.014%)

etc.

(I’ve not checked the quality of the output, though!)

sbembenek · October 19, 2021, 11:02pm

Thanks – this did work for me.

What do you mean ‘quality’ of the output?

laz · October 19, 2021, 11:09pm

I was just meaning the level of theory and any match with experiment, etc.

sbembenek · October 19, 2021, 11:23pm

Yes, understood Thanks again.

sbembenek · October 20, 2021, 11:15pm

One thing I noticed is that while this (‘your syntax’) works for me, such that I get the the TDMs and I am able to control the coefficient cutoff:

molecule mol {
0 1
O 0.00000000 0.00000000 0.00000000
H 0.00000000 -0.79069000 0.61221700
H 0.00000000 0.79069000 0.61221700

 symmetry c1

}

set scf_type df
set reference rks
set s_tolerance 1e-9
set_num_threads(64)

set save_jk true

set {
  tdscf_states 20
  tdscf_coeff_cutoff 0.01
  tdscf_tdm_print ['E_TDM_LEN', 'E_TDM_VEL']
}

e, wfn = energy('td-b3lyp-d3bj/aug-cc-pVTZ', return_wfn=True)

This syntax does not allow me to get the TDMs nor control coefficient cutoff:

molecule mol {
0 1
O 0.00000000 0.00000000 0.00000000
H 0.00000000 -0.79069000 0.61221700
H 0.00000000 0.79069000 0.61221700

 symmetry c1

}

set scf_type df
set reference rks
set s_tolerance 1e-9
set_num_threads(64)

set save_jk true

import psi4

from psi4.driver.procrouting.response.scf_response import tdscf_excitations

set {
 
  tdscf_coeff_cutoff 0.01
  tdscf_tdm_print ['E_TDM_LEN', 'E_TDM_VEL']
}

e, wfn = energy('b3lyp-d3bj/aug-cc-pVTZ', return_wfn=True)
res = tdscf_excitations(wfn, states=20)

The latter is desirable because I am then able to use the ‘spectrum’ function.

In short, I have not been able to figure out how to bring together getting all these results in a single TDDFT. I appreciate we’ve gone a ‘bit’ past our original discussion but wanted to hit you up first before I posted as new topic.

laz · October 21, 2021, 1:55pm

Ahhhh…I see what you mean. Something like this should do the trick:

e, wfn = energy('b3lyp-d3bj/aug-cc-pVTZ', return_wfn=True)
res = tdscf_excitations(wfn, states=20, tdm_print=['E_TDM_LEN', 'E_TDM_VEL'], coeff_cutoff=0.01)

(Modifying my benzene example above in that way gave me the same output as above with the “td-b3lyp” sort of usage.)

sbembenek · November 4, 2021, 8:04pm

Hey,

Can you clarify why the TDMs in the length versus velocity form are different for a given excited state? This isn’t readily obvious to me. So, are the units then different for these two forms? Thanks.

hokru · November 5, 2021, 8:27am

Due to how they are formulated they are different for approximate wave functions.
The length form depends on the coordinate system (it’s gauge dependent) while the velocity is not.
However, the length forms is less sensitive to the quality of the wave function, which is why for TD-DFT the length form is more often used.

Opinions what to use may differ from group to group though.

sbembenek · November 10, 2021, 12:07am

Thanks.

I picked up some papers to get more up to speed (do you have papers you’d suggest?)

Thanks again for getting me started.