Excited States & Transition Dipole with TD-DFT

sbembenek · October 9, 2021, 3:33am

I’m trying to get excited states and the transition dipole moment using TD-DFT

Here’s an example of what I am doing:

molecule mol {
0 1
O 0.00000000 0.00000000 0.00000000
H -0.75725300 -0.00000000 -0.58843500
H 0.75725300 0.00000000 -0.58843500

}

set scf_type df
set reference rks
set s_tolerance 1e-9

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

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

mol.print_out()
print_variables()

The key output is then:

                         TDSCF excitation energies                
                 by Andrew M. James and Daniel G. A. Smith        
         ---------------------------------------------------------
********************************************************************************
**********                          WARNING                           **********
**********      TDSCF is experimental results may be inaccurate       **********
********************************************************************************

  ==> Requested Excitations <==

      0 states with A1 symmetry
      0 states with A2 symmetry
      0 states with B1 symmetry
      0 states with B2 symmetry

  ==> Options <==

     etol      :              1e-06
     rtol      :              1e-08
     guess_type:              denominators
     restricted:              True
     triplet   :              False
     ptype     :              rpa



Final Energetic Summary:
                                    Excitation Energy         Total Energy  
     #   Sym: GS->ES (Trans)       [au]            [eV]            (au)

In other words, no result is found.

The similar calculation with EOM-CC2 as:

0 1
O 0.00000000 0.00000000 0.00000000
H -0.75725300 -0.00000000 -0.58843500
H 0.75725300 0.00000000 -0.58843500

}

set scf_type pk
set basis aug-cc-pVTZ
set reference rhf


e, wfn = energy('cc2', return_wfn=True)

set roots_per_irrep [0, 0, 2, 2]

properties('eom-cc2', properties=['oscillator_strength'])

mol.print_out()
print_variables()

Gives ‘desirable’ results, for example:


	                   Ground State -> Excited State Transitions

	                   Excitation Energy          OS       RS        RS     Einstein A
	State   (eV)    (cm^-1)    (nm)     (au)              (l,au)   (v,au)     (s^-1)
	 1 B1  11.207   90389.9   110.6   0.411847   0.0142   0.0000   0.0000  7.731354E+07
	 2 B1  13.017  104987.9    95.2   0.478360   0.0832   0.0000   0.0000  6.116557E+08
	 1 B2   7.244   58424.9   171.2   0.266203   0.0578   0.0000   0.0000  1.316965E+08
	 2 B2  10.376   83686.2   119.5   0.381302   0.0006   0.0000   0.0000  2.774351E+06

Doing transition

	*********************************************************************
	*********************************************************************
	************                                             ************
	************ Excited State-Excited State Transition Data ************
	************                                             ************
	*********************************************************************
	*********************************************************************


	*** Computing <1B1|X{pq}}|2B1> (LEFT) Transition Density ***

		*** LTD Setup complete.
		***...density has been built...
		***...density has been sorted...

	*** Computing <2B1|X{pq}}|1B1> (RIGHT) Transition Density ***

		*** RTD Setup complete.
		***...density has been built...
		***...density has been sorted...

What needs to be different in the TD-DFT setup to get the transitions and corresponding transition dipole moments?

thanks!

jmisiewicz · October 9, 2021, 11:51am

It looks like you’re using a pre-1.4 version of Psi. Upgrade to 1.4 and follow the example in the docs.

sbembenek · October 10, 2021, 7:24pm

thanks.

I’ve upgraded to the 1.4.1.

Here’s the input file:

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

}

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


mol.update_geometry()
mol.symmetrize(1e-3)

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

e, wfn = optimize('b3lyp-d3bj', return_wfn=True)
res = tdscf_excitations(wfn, states=8)

mol.print_out()
print_variables()

However, I now get:

        ---------------------------------------------------------
                         TDSCF excitation energies                
                 by Andrew M. James and Daniel G. A. Smith        
         ---------------------------------------------------------

  ==> Options <==

     Residual threshold  : 1.0000e-04
     Initial guess       : denominators
     Reference           : RHF
     Solver type         : RPA (Hamiltonian)


  ==> Requested Excitations <==

      2 singlet states with A1 symmetry
      2 singlet states with A2 symmetry
      2 singlet states with B1 symmetry
      2 singlet states with B2 symmetry


  ==> Seeking the lowest 2 singlet states with A1 symmetry

                        Generalized Hamiltonian Solver                       
                              By Andrew M. James                             

  ==> Options <==

    Max number of iterations        = 60   
    Eigenvector tolerance           = 1.0000e-04
    Max number of expansion vectors = 400  

  => Iterations <=
                              Max[D[value]]     Max[|R|]   # vectors

Traceback (most recent call last):
  File "/data/Apps/anaconda3/envs/psi4env/bin/psi4", line 333, in <module>
    exec(content)
  File "<string>", line 36, in <module>
  File "/data/Apps/anaconda3/envs/psi4env/lib//python3.7/site-packages/psi4/driver/procrouting/response/scf_response.py", line 742, in tdscf_excitations
    res_1 = _solve_loop(wfn, ptype, solve_function, singlets_per_irrep, maxiter, restricted, "singlet")
  File "/data/Apps/anaconda3/envs/psi4env/lib//python3.7/site-packages/psi4/driver/procrouting/response/scf_response.py", line 348, in _solve_loop
    ret = solve_function(engine, nstates, guess_, maxiter)
  File "/data/Apps/anaconda3/envs/psi4env/lib//python3.7/site-packages/psi4/driver/procrouting/response/scf_response.py", line 717, in rpa_solver
    verbose=verbose)
  File "/data/Apps/anaconda3/envs/psi4env/lib//python3.7/site-packages/psi4/driver/p4util/solvers.py", line 961, in hamiltonian_solver
    H1x, H2x, nprod = engine.compute_products(vecs)
  File "/data/Apps/anaconda3/envs/psi4env/lib//python3.7/site-packages/psi4/driver/procrouting/response/scf_products.py", line 260, in compute_products
    twoel = self.wfn.twoel_Hx(compute_vectors, False, "SO")

RuntimeError: 
Fatal Error: RHF::twoel_Hx: JK object is not initialized, please set option SAVE_JK to True.
Error occurred in file: /scratch/psilocaluser/conda-builds/psi4-multiout_1633639682895/work/psi4/src/psi4/libscf_solver/rhf.cc on line: 456
The most recent 5 function calls were:




Printing out the relevant lines from the Psithon --> Python processed input file:
    set_num_threads(62)
    mol.update_geometry()
    mol.symmetrize(1e-3)
    from psi4.driver.procrouting.response.scf_response import tdscf_excitations
    e, wfn = optimize('b3lyp-d3bj', return_wfn=True)
--> res = tdscf_excitations(wfn, states=8)
    mol.print_out()
    print_variables()

!----------------------------------------------------------------------------------!
!                                                                                  !
! Fatal Error: RHF::twoel_Hx: JK object is not initialized, please set option      !
!     SAVE_JK to True.                                                             !
! Error occurred in file: /scratch/psilocaluser/conda-                             !
!     builds/psi4-multiout_1633639682895/work/psi4/src/psi4/libscf_solver/rhf.cc   !
!     on line: 456                                                                 !
! The most recent 5 function calls were:                                           !
!                                                                                  !
!----------------------------------------------------------------------------------!

ajs99778 · October 10, 2021, 8:04pm

Using set save_jk true should get this to run. Here is a similar input file that works for me:

#asdf.xyz
set_num_threads(6)
memory 6 GB
from psi4.driver.procrouting.response.scf_response import tdscf_excitations
basis {
    assign    cc-pVDZ
}

molecule {
     0 1
     symmetry c1
     O        0.00000        0.69500       -0.09249
     O       -0.00000       -0.69500       -0.09249
     H       -0.38814        0.89525        0.73989
     H        0.38814       -0.89525        0.73989
}

set {
    save_jk                 true
    scf_type                df
    reference               rks
    s_tolerance             1e-9
}

mol = get_active_molecule()

nrg, wfn = optimize('B3LYP', return_wfn=True)
res = tdscf_excitations(wfn, states=10)

mol.print_out()
print_variables()

sbembenek · October 11, 2021, 5:12am

thanks – it worked.

Also, do I need to add something particular if I want transitions other than from the ground state? Or is it simply a matter of adding more ‘states’? In particular, I’m especially interested in transitions around the HOMO-LUMO MOs.

thanks again!

ajs99778 · October 11, 2021, 6:16pm

I do not know how to choose the initial state with Psi4’s TD-DFT. As far as I know, it is always the ground state.

sbembenek · October 11, 2021, 8:13pm

thanks for your help!

jmisiewicz · October 11, 2021, 11:31pm

While I’m not an expert on excited states, I’m fairly sure this is one of the inherent limitations of TD-DFT and response theory for excited states: it only describes excitations from the state you start with. In principle, you could converge the excited state density. Try MOM for that.

Consult an expert if you want to pursue this further.

sbembenek · October 12, 2021, 2:32am

thanks for the help.

But you do get excited state → excited state from EOM-CC - yes?

If you have an ‘expert’ in mind, please point me in the right direction.

thanks again.

robertodr · October 12, 2021, 1:48pm

@jmisiewicz is right: you cannot compute transition properties between excited states using linear response. You would need the quadratic response function and its double residues to access those quantities. You can try to use a different ground state, i.e. run a DFT calculation for an excited state and then run TDDFT on top. However, I suspect this kind of calculation would fail: the iterative subspace diagonalization cannot cope with non-positive definite matrices and a non-ground state reference is bound to be a saddle point for the electronic Hessian.

The EOM-CC methods can compute excited state-excited state transition properties because they do build actual representations of the excited states while solving for the requested amount of excitation energies. See also Figure 2 in this review The review offers a very thorough overview of response theory, though it’s a bit math-heavy.

sbembenek · October 12, 2021, 7:39pm

thanks for the clarification and I’ll check out the review.

A followup question:

So it sounds like I cannot get the excited-state geometries form TDDFT – it this correct?

Can I get the molecular configuration such as to know which electronic transition occurs, or do you only get the molecular state to state transition?

jmisiewicz · October 12, 2021, 11:49pm

Analytic gradients are not supported, no. They’re theoretically possible, but just not implemented. Getting an excited-state geometry by a finite difference optimization is possible, but if that’s what you want, I would need to check how you actually do that. It may require a very ugly input file.
I don’t understand the question. By “molecular configuration,” do you mean a geometry? An MO picture of the transition? How does “molecular configuration” differ from “molecular state”?

sbembenek · October 13, 2021, 12:33am

Thanks – I’ll put this on the back burner for now.
I’m interested in knowing what specific electronic transition occurred (e.g., HOMO → LUMO). I don’t know if this possible with TDDFT.

For example:

                                    Excitation Energy         Total Energy        Oscillator Strength             Rotatory Strength       
     #   Sym: GS->ES (Trans)        au              eV              au          au (length)    au (velocity)    au (length)    au (velocity) 
    ---- -------------------- --------------- --------------- --------------- --------------- --------------- --------------- ---------------
     1      A1->B1 (1 B1)      0.25391         6.90912        -76.21287        0.0471          0.0473          0.0000         -0.0000

I’m reading this as: the molecule with an overall (state) symmetry of A1 has a certain probability of making a single-electron transition to another state with overall (state) symmetry of B1. By molecular configuration, I simply meant the actual electronic transition that occurred or the MOs that the electrons are occupying in the transition state.

Thanks again for the help – I really appreciate it.

jmisiewicz · October 13, 2021, 1:02am

Ah, I see.

In the TD-DFT picture, an excited state doesn’t unambiguously correspond to a single excitation of the MOs. Instead, all excitations have some “weight”. There are two approaches from here:

Stay in the DFT MO picture, and interpret your excitation in terms of the excitations with largest weights.
Get a picture in terms of very high weight excitations of orbitals called the natural transition orbitals. These are not the DFT MOs and are computed by applying an SVD to the weight vector.

(2) isn’t implemented in Psi. (1) is implemented, but only with C1 symmetry, and it’s only done automatically in the TDA approximation for reasons I don’t understand.

I should make a “TD-DFT expansion wishlist” issue…

sbembenek · October 13, 2021, 1:36am

Thanks for clarifying – helps a lot!

I’m good with option #1 – but do add #2 to the wish list

Oh and for EOM-CC? No go with getting the particular electronic transition out there?

system · December 12, 2021, 1:37am

This topic was automatically closed 60 days after the last reply. New replies are no longer allowed.