LR-TDDFT of excited state


Just Got Here
Dear all,
I am trying to calculate the absorption spectrum of excited CHBr3 with an extended C-Br bond. CHBr3 is in excited T2, where one HOMO-1 electron is excited to LUMO. I can use the following tag to do it, however, it doesn't necessarily converge to T2 state. Sometimes, it only converges to T1 state.

vectors swap beta 55 56
MAX_OVL

Is it possible to run a first LR-TDDFT calculation with tag CIVECS to export the T2 state and perform the second LR-TDDFT based on the CI vectors from the first step?

If we have to revise the code, what subroutines I should revise?

Thanks a lot!

Forum Regular
No, you cannot perform a second LR-TDDFT calculation with a reference built from the CI vectors of an excited state. LR-TDDFT is formulated for a single determinant reference state and your T2 solution from the first LR-TDDFT would be a linear combination of determinants.

Quadratic response TDDFT could give you the absorption spectrum for an excited state; however, it is not implemented in NWChem (the only implementation I am aware of is in the code Dalton).

Within NWChem, you could attempt our RT-TDDFT based approach to excited state absorption ("Excited state absorption from real-time time-dependent density functional theory" JCTC vol. 11 4294-4303 [2015] https://pubs.acs.org/doi/abs/10.1021/acs.jctc.5b00473). While not documented, the functionality for this approach is in version 6.8 of the code. You would need to first run a LR-TDDFT gradient calculation on the excited state of interest in order to generate the excited state density matrix (or matrices in the case of an open shell calculation). You can then use the excited state density as your starting point by adding the appropriate load command to your RT-TDDFT input block:

load density <file name>.dmat

for closed shell or

load density <file name>.dmatA <file name>.dmatB

for open shell, where you would obviously insert the file name for the density matrix (matrices) generated from the LR-TDDFT gradient calculation. From there you should be able to follow the method outlined in the above cited paper to generate the absorption spectrum.

Just Got Here
Hi..in my case i was trying to do excited state geometry optimization using the analytical gradient of LR-TDDFT. I did not find the relevant part of the manual that describes the gradient calculation of LR-TDDFT. so, I am wondering if this function available in dalton? if yes, how to use it? Thank you for your help!

https://www.7pcb.com/

Forum Regular
Sorry I misunderstood you. Geometry optimization on the excited state is straightforward within NWChem. See the third example input in the documentation https://github.com/nwchemgit/nwchem/wiki/Excited-State-Calculations#sample-input

Clicked A Few Times
Quote:Sean Dec 7th 4:50 am
No, you cannot perform a second LR-TDDFT calculation with a reference built from the CI vectors of an excited state. LR-TDDFT is formulated for a single determinant reference state and your T2 solution from the first LR-TDDFT would be a linear combination of determinants.

Quadratic response TDDFT could give you the absorption spectrum for an excited state; however, it is not implemented in NWChem (the only implementation I am aware of is in the code Dalton).

Within NWChem, you could attempt our RT-TDDFT based approach to excited state absorption ("Excited state absorption from real-time time-dependent density functional theory" JCTC vol. 11 4294-4303 [2015] https://pubs.acs.org/doi/abs/10.1021/acs.jctc.5b00473). While not documented, the functionality for this approach is in version 6.8 of the code. You would need to first run a LR-TDDFT gradient calculation on the excited state of interest in order to generate the excited state density matrix (or matrices in the case of an open shell calculation). You can then use the excited state density as your starting point by adding the appropriate load command to your RT-TDDFT input block:

load density <file name>.dmat

for closed shell or

load density <file name>.dmatA <file name>.dmatB

for open shell, where you would obviously insert the file name for the density matrix (matrices) generated from the LR-TDDFT gradient calculation. From there you should be able to follow the method outlined in the above cited paper to generate the absorption spectrum.


Hi Sean,
Could you please elaborate a bit on the RT-TDDFT with an example input?
Also once the ESA spectrum is obtained, can you give an example input on how to plot the RT-TDHF transition densities for each of the features of the spectrum?
Thanks in advance

Forum Regular
Below is an example for calculating the excited state response with RT-TDDFT for the 5th excited state of a water molecule. Some items to note:
-A LR-TDDFT gradient calculation needs to be performed first in order to generate the excited state density matrix
-That density matrix needs to be loaded within the RT-TDDFT block of the input
-In addition to the RT-TDDFT runs with field applied (here I have only done the excitation along the x-axis), an additional run needs to be done without any field applied; this extra run is used to correct the dipole moments from the runs with an applied field

Obtaining transition densities to plot is significantly more involved than just the spectrum. You need to add a visualization block to each RT-TDDFT input block. This will generate a density matrix file for each step in each RT-TDDFT run. Then as with the dipole moment, the run without a field is used to correct the runs with a field. In the examples that we showed in the paper, the probed transitions were polarized along a single axis so we just used the density matrices from the run with a field along that axis. After a series of corrected density matrices are generated, then you need to take the Fourier transform of the density matrices at the transition frequency. The resulting density matrix can then be fed into DPLOT to generate a cube file for visualization.

echo
start water

geometry "system" units angstroms nocenter noautoz noautosym
 O     0.00000000     0.00000000     0.15104872
 H    -0.75808206     0.00000000    -0.48477224
 H     0.75808206     0.00000000    -0.48477224
end

set geometry "system"

basis
 * library STO-3G
end

dft
  xc hfexch 1.0
end

tddft
 nroots 10
 notriplet
 civecs
 grad
  root 5
 end
end
task tddft gradient

unset rt_tddft:*
rt_tddft
  tmax 200.0
  dt 0.2
  load density water.dmat
  tag "ref"
  print dipole
 end
task dft rt_tddft

unset rt_tddft:*
rt_tddft
  tmax 200.0
  dt 0.2
  load density water.dmat
  tag "kick_x"
  print dipole
  field "kick"
    type delta
    polarization x
    max 0.0001
  end
  excite "system" with "kick"
 end
task dft rt_tddft

Clicked A Few Times
Thanks for your prompt reply Sean,
Do I need to add this additional no-field step at every polarization x, y and z?

Would the method work after a tddft optimization such as below?
tddft
nroots 10
notriplet
civecs
target 5
targetsym a
grad
root 5
end
end
task tddft optimize

What is the difference between the two calculations? Does the gradient calculation produce a spectrum near the vertical excitation while the other at the equilibrated state?

As for the visualization, would something like this work:

unset rt_tddft:*
rt_tddft
 tmax 200.0
dt 0.2
load density water.dmat
tag "ref"
print dipole
end
task dft rt_tddft

unset rt_tddft:*
rt_tddft
 tmax 200.0
dt 0.2
load density water.dmat
tag "kick_x"
print dipole
field "kick"
type delta
polarization x
max 0.0001
end
excite "system" with "kick"
end
visualization
 tstart 0.0
tend 200.0
treference 0.0
dplot
end
task dft rt_tddft

What FFT utility would you use for the fourier transform?

Thanks for all your help

Forum Regular
You only need to run the calculation without a field once.

Running an optimization first would be fine. My example input assumed you already had the geometry of interest. The reason it was set as a call to tddft gradient is that the excited state density matrix is not computed if you are only calculating the excitation energies; it is only calculated if you need an excited state gradient. If you are doing a geometry optimization on the excited state, then you will be calculating excited state gradients and will be generating the excited state density matrix. The caveat to this is that the XC functional must have analytical 2nd derivatives implemented; otherwise, the code is calculating gradients through numerical differentiation and an excited state density matrix is not calculated. For those functionals (e.g. M06), we showed that one possible workaround is to use a different functional to calculated the excited state density matrix and then use the functional of interest for the RT-TDDFT propagation (see, https://pubs.acs.org/doi/abs/10.1021/acs.jpclett.6b00282)

You need a visualization block for each propagation, including the one without a field. Also, don't include the treference keyword or the dplot keyword.

Use whatever FFT utility works for you. In the past I wrote analysis code in Fortran around FFTPACK from netlib, but most programing languages have FFT libraries that can be linked.

Clicked A Few Times
Once the calculaton is finished, how would you amend these commands to include the reference without the field?

Extracting the x-dipole moment for the x-kick:
nw_rtparse.py -xdipole -px -tkick_x input.nwo > x.dat

FT for x polarization:
fft1d -d50 -z -p50000 < x.dat | rotate_fft > xw.dat

Forum Regular
One possibility is

run:
nw_rtparse.py -xdipole -px -tkick_x input.nwo > x.dat

then:
nw_rtparse.py -xdipole -px -tref input.nwo > refx.dat

then:
paste x.dat refx.dat | tail -n <number of data points> | awk '{print $1,$2-$4}' > corrected_x.dat


The top of the x.dat/refx.dat files should tell you how many data points you have. You should then be able to feed the corrected_x.dat file into the fft1d utility. You can repeat the same for the other polarizations, e.g.

run:
nw_rtparse.py -xdipole -py -tkick_y input.nwo > y.dat

then:
nw_rtparse.py -xdipole -py -tref input.nwo > refy.dat

then:
paste y.dat refy.dat | tail -n <number of data points> | awk '{print $1,$2-$4}' > corrected_y.dat


The point is just that you want to subtract the no-field propagated dipole from the field propagated dipole at each time step before doing the FT.

Clicked A Few Times
Thank you Sean for the clear explanation. Meanwhile I have some troubles with the rotate_fft.py script:
I cannot use it together with fft1d.m as suggested in the H20 example. The command below,
octave -q fft1d.m < x.dat | rotate_fft.py > xw.dat
produces the error:
error: fft1d.m: A(I): index out of bounds; value 1 out of bound 0
error: called from
   fft1d.m at line 41 column 7

While octave -q fft1d.m x.dat > xw.dat works fine, and trying to do the rotate on the resulting file xw,dat also does not work
rotate_fft.py xw.dat > xw_rot.dat

Forum Regular
Let me begin by saying that the scripts in $NWCHEM_TOP/contrib/parsers are not mine, but I think I have managed to get them to work.

As you have already discovered, the fft1d.m utility takes the input as a command line argument rather than a redirect from stdin. Also note that the documentation for RT-TDDFT says that there are a couple parameters hardcoded into that utility, specifically, the time constant for the exponential damping of the signal and the amount of zero-padding. You may want to considering playing around with those numbers to improve the appearance of your output.

To use the rotate_fft.py script on the output of the fft1d.m utility, you need to first delete the header of that xw.dat file (i.e. delete all the text at the top so that the file only contains the 4 columns of numbers. Then rotate_fft.py takes its input from a redirect rather than the command line, i.e. rotate_fft.py < xw.dat (alternatively, python rotate_fft.py < xw.dat, if you haven't made rotate_fft.py executable).

Clicked A Few Times
Thanks Sean rotate_fft works now with the redirect.

I now encountered a problem when using damping values >5.0 and for the first points of my data. rotate_fft.py fails with the exception:
 ("abs not equal to sqrt(re^2 + im^2)") 
In order for rotate_fft to run through all the data, I have to either discard the first ~300 points or loosen the difference between Abs and the sum of square roots of Im and Re to < 1e-5 instead of the default < 1e-6 in rotate_fft.py. By adding more damping the error is increased and (Abs - sqrt(Im^2+ Re^2) has to be loosened further.
In your experience, do you think this is normal? What does this say about the quality of the data?

Forum Regular
This is not normal. That part of the rotate_fft script is just checking the data that you are inputing, so if you are getting that error message, something has gone wrong with the previous steps. Have you looked at your time-dependent dipole moment to check that it looks reasonable? Have you looked at the output of fft1d.m to see if that looks reasonable? For the example input I provided above, I am able to run both fft1d.m and rotate_fft.py without issue, have you been able to do the same?


Forum >> NWChem's corner >> General Topics