Understanding MRSF TDDFT - Open-Quantum-Platform/openqp GitHub Wiki
- Excited State Derivation: TDDFT allows us to derive excited states as variations or responses of the ground state, eliminating the need for additional methods to obtain excited states. In essence, any excited state can be obtained as a "response state" of the ground state.
- Ground State as a Response State: TDDFT extends its utility beyond excited states; it also enables the ground state to be obtained as a response state using spin-flip techniques.
However, despite its popularity, the Linear-Response (LR)-TDDFT method has some well-documented limitations:
- Non-Local Properties: LR-TDDFT may fail to capture non-local properties, particularly for long-range charge transfer excited states.
- Double Excitation Characters: It may provide a poor description of excited states with double excitation character.
- Static Correlation: LR-TDDFT may struggle to accurately describe static correlation in closed-shell reference states that undergo bond-breaking.
- Lack of Coupling: It lacks the capability to model the coupling between ground and excited states, which is essential for conical intersections (CI) and avoided crossings.
Researchers and scientists should be aware of these limitations when applying LR-TDDFT and consider alternative methods or approaches for cases where these limitations may impact the accuracy of their calculations.
By incorporating these "red missing" configurations, SF-TDDFT aims to provide a more complete and accurate description of excited states, mitigating the spin contamination issue associated with the conventional approach.
Yet another important but not much-appreciated source of spin contamination is from the mismatched contributions of L and R of OO TYPE (the blues). This is because they are coming from different orbitals' spin-flip transitions as shown below. The former and latter come from the spin-flip of O1 and O2 orbitals, respectively. Sometimes, the mismatched contributions introduce a major spin contamination.
The red missing configurations can be added into the response space by the (M_S = -1) component of ROHF. A hypothetical single reference combining (M_S = +1) and (M_S = -1) components of the ROHF triplet can be constructed by a spinor-like transformation. See more here. The resulting mixed-reference SF-TDDFT (MRSF-TDDFT) eliminates the spin contamination of SF-TDDFT, allowing automatic identification of the electronic states as singlets and triplets. It should be emphasized that MRSF-TDDFT produces not only excited but also ground electronic states. Therefore, open-shell singlets such as diradicals, which cannot be studied by Kohn-Sham DFT, can be naturally described by MRSF-TDDFT.
There are multiple advantages of MRSF-TDDFT over TDDFT. Here, we list just some of them.
In the case of C and N+ atoms, MRSF-TDDFT produces 8 states of ( ^3 P, ^1 D, ^1 S ) as:
On the other hand, only 4 states can be represented by TDDFT as:
This demonstrates that the conventional TDDFT misses many electronic states. Especially, the doubly excited configuration of ( ^1 S ) is completely missing in TDDFT.
The (3 \Sigma_g^{+}) state (in green), entirely composed of a doubly excited configuration, is missing in TDDFT (LR) but present in MRSF-TDDFT as shown below in the case of (H_2) dissociation.
The doubly excited configuration makes a greater contribution to the ground electronic state in the much-stretched region, eventually leading to the correct dissociation limit shown by the dashed line. On the other hand, the corresponding ground state of the LR-PBE0 curve (TDDFT) has to be obtained by DFT, which does not have the correct asymptotic behavior. See more here.
In linear all-trans polyenes, internal conversion (IC) between the (1^1B_u^+) (optically bright at the Franck-Condon (FC) geometry, the red curve below) and the (2^1A_g^-) (optically dark, the blue curve) states has long been debated. However, it hasn't been proven until recently.
One of the difficulties in describing the (1^1B_u^+) and the (2^1A_g^-) states is their radically different nature. The former state comprises a one-electron HOMO (\rightarrow) LUMO transition and displays pronounced ionic characteristics, while the latter is dominated by HOMO (\rightarrow) LUMO double excitations, requiring a balanced theory with both dynamic and nondynamic electron correlation.
As seen above, while both MRSF-TDDFT and REKS(4,4) produce the correct CI as seen in (\delta)-CR-EOMCC(2,3), the TDDFT does not. See more here.
A Conical Intersection (CI) is a molecular geometry at which two (or more) adiabatic electronic states become degenerate. As the intersecting electronic states at a CI are coupled by non-adiabatic coupling, CIs provide efficient funnels for state-to-state population transfer mediated by nuclear motion. The degeneracy of the intersecting states at a CI is lifted along two directions in the space of internal molecular coordinates Q, which are defined by the gradient difference and derivative coupling vectors (GDV and DCV, respectively):
This is for the case of a crossing between the ground ((S_0)) and the lowest excited (S_1) singlet states. The degeneracy is lifted linearly along the GDV and DCV directions, which lends the potential energy surfaces (PESs) of the intersecting states to the topology
of a double cone below, hence the name. The remaining 3N-8 internal coordinates leave the degeneracy intact, thus defining the crossing seam (or the intersection space) of the CI.
The popular TDDFT fails to yield the correct dimensionality of the CI seam and predicts a linear crossing (3N-7) as shown in the right figure above. The problem is simply rooted in the absence of coupling between the ground (described by reference DFT) and first excited states (described by the response of TDDFT).
On the other hand, both ground and first excited states of MRSF-TDDFT are obtained by the same response, producing the correct double cone topology. See more here.
MRSF-TDDFT currently calculates singlet, triplet, and quintet states using the ROHF triplet reference. As a result, one may be confused by the two different triplets of reference and response. A triplet has three states with (M_s = +1, 0, -1). Since they are energetically degenerate, nearly all quantum chemistry programs calculate the simplest possible state of (M_s = +1) as a representative triplet state of ROHF. This particular triplet is the high-spin triplet. On the other hand, the triplet states generated by MRSF-TDDFT are (M_s = 0) triplets. Both the _reference_ and lowest _response_ triplets are supposed to describe the same lowest triplet states with just different (M_s). However, their energies are not degenerate. This is because the _reference_ triplet is variationally obtained, while the _response_ triplet is calculated by linear response theory. When you use MRSF-TDDFT, it is always better not to utilize the _reference_ triplet in your analysis.
This simple example can explain a great deal of MRSF-TDDFT calculations. The input example of MRSF-TDDFT with BHHLYP/6-31G basis set for OpenQP is:
[input]
system=
4 0.000000000 0.000000000 0.000000000
charge=0
runtype=energy
basis=6-31g
functional=bhhlyp
method=tdhf
[guess]
type=huckel
[scf]
multiplicity=3
type=rohf
[tdhf]
type=mrsf
multiplicity=1
nstate=5
The important keywords are:
-
[tdhf], type=mrsf: Specifying MRSF-TDDFT method -
[tdhf], nstate=5: Total number of excited states requested -
[scf], multiplicity=3and[tdhf], multiplicity=1:multiplicityspecifies spin multiplicity. Singlet and triplet are 1 and 3, respectively. As you can see, there are 2 different places ofmultiplicity. The one in[scf]group specifies the spin multiplicity of the reference wavefunction, in this case, the ROHF reference state. Themultiplicityin[tdhf]group specifies the spin multiplicity of response states, which are the ones generated by MRSF-TDDFT theory.
This input will print out the SCF procedures as:
Direct SCF iterations begin.
=============================================================================================
Iter Energy Delta E Int Skip DIIS Error Shift Method
=============================================================================================
1 -14.2146699972 -14.2146699972 0 0.24795652 0.000 SD
2 -14.5596466705 -0.3449766733 0 0.04973230 0.000 C-DIIS
3 -14.5664546039 -0.0068079334 0 0.00556032 0.000 C-DIIS
4 -14.5665968308 -0.0001422269 0 0.00046311 0.000 C-DIIS
5 -14.5665974671 -0.0000006363 0 0.00005130 0.000 C-DIIS
6 -14.5665974743 -0.0000000072 0 0.00000247 0.000 C-DIIS
7 -14.5665972410 0.0000002333 0 0.00029265 0.000 C-DIIS
8 -14.5665974717 -0.0000002307 0 0.00003074 0.000 C-DIIS
9 -14.5665974740 -0.0000000023 0 0.00000669 0.000 C-DIIS
10 -14.5665974741 -0.0000000001 0 0.00000001 0.000 C-DIIS
----------------------------------------------------------------
SCF convergence achieved ....
Final ROHF energy is -14.5665974741 after 10 iterations
DFT: XC energy = -1.4315120718
DFT: total electron density = 3.9999999987
DFT: number of electrons = 4
As you can see above, the final ROHF-DFT with BHHLYP functional energy is -14.5665974741 Hartree, which corresponds to the (1s^2 2s^1 2p^1) electron configuration. (The ground-state electron configuration of the Be atom is (1s^2 2s^2)). There are 3 (2p) orbitals, where one can put an electron. Since ROHF chooses one of the three possible (2p) orbitals, it inevitably breaks the symmetry of them and makes one of them more stable than the other two, which is exactly seen in the final orbitals below.
===============================
Molecular Orbitals and Energies
===============================
-------------- Alpha Orbitals -------------
1 2 3 4 5
-4.3406033266 -0.1902242226 -0.0441863943 0.0130839299 0.0130839455
1 Be 1 S 0.9967308542 -0.2281711167 -0.0000000000 -0.0000000000 0.0000000000
2 Be 1 S 0.0276583420 0.3642811412 0.0000000000 -0.0000000000 0.0000000000
3 Be 1 X 0.0000000000 0.0000000000 -0.5629800381 0.0000161123 -0.0000000000
4 Be 1 Y -0.0000000000 -0.0000000000 -0.0000000000 -0.0000009912 -0.4107251178
5 Be 1 Z 0.0000000000 0.0000000000 0.0000221305 0.4107251715 -0.0000009912
6 Be 1 S -0.0105052228 0.6891136408 0.0000000000 0.0000000000 -0.0000000000
7 Be 1 X -0.0000000000 0.0000000000 -0.5315765683 0.0000265969 -0.0000000001
8 Be 1 Y 0.0000000000 0.0000000000 -0.0000000000 -0.0000016348 -0.6774590577
9 Be 1 Z -0.0000000000 0.0000000000 0.0000208319 0.6774590081 -0.0000016348
6 7 8 9
0.3350537379 0.3379393730 0.3379393752 0.3602164892
1 Be 1 S -0.0000000000 -0.0000000000 0.0000000000 0.0070794669
2 Be 1 S -0.0000000000 -0.0000000000 0.0000000000 2.0048584067
3 Be 1 X -1.2221506776 -0.0000565098 0.0000000005 -0.0000000000
4 Be 1 Y 0.0000000000 0.0000101635 1.2813678952 -0.0000000000
5 Be 1 Z 0.0000544684 -1.2813678766 0.0000101635 -0.0000000000
6 Be 1 S 0.0000000000 0.0000000000 -0.0000000000 -1.9255515375
7 Be 1 X 1.2361331096 0.0000518900 -0.0000000004 0.0000000000
8 Be 1 Y -0.0000000000 -0.0000092215 -1.1626039485 0.0000000000
9 Be 1 Z -0.0000544214 1.1626039760 -0.0000092215 0.0000000000
Even though the three orbitals (3, 4, and 5 above) are degenerate, the orbital energy of 3 is -0.0441863943, while those of 4 and 5 are 0.0130839455 as a result of the ROHF calculation above.
The major transition contributions of each state by MRSF-TDDFT are seen from:
===================================
Spin-adapted spin-flip excitations
===================================
State # 1 Energy = -2.336035 eV
<S^2> = 0.0000
DRF Coeff OCC VIR
--- -------- ------ ------
3 -0.992115 3 -> 2
5 -0.124395 2 -> 3
State # 2 Energy = 2.340563 eV
<S^2> = 0.0000
DRF Coeff OCC VIR
--- -------- ------ ------
9 -0.221863 3 -> 4
12 0.972646 3 -> 5
21 -0.067094 3 -> 8
State # 3 Energy = 2.340563 eV
<S^2> = 0.0000
DRF Coeff OCC VIR
--- -------- ------ ------
9 -0.972646 3 -> 4
12 -0.221863 3 -> 5
18 0.067094 3 -> 7
State # 4 Energy = 2.571426 eV
<S^2> = 0.0000
DRF Coeff OCC VIR
--- -------- ------ ------
2 -0.998412 2 -> 2
15 -0.056300 3 -> 6
State # 5 Energy = 5.042455 eV
<S^2> = 0.0000
DRF Coeff OCC VIR
--- -------- ------ ------
8 -0.992697 2 -> 4
17 0.120630 2 -> 7
As seen above, the first two response states are listed as STATE #1 and #2 with the coefficients of major contributions. For example, the major transition of STATE #1 is 3 (OCC) (\rightarrow) 2 (VIR) spin-flip transition, which means that the (\alpha) spin electron in orbital #3 goes to (\beta) spin electron in orbital #2 with a coefficient of -0.992115. As a result of this transition, orbital #2 now becomes doubly occupied with (\alpha) and (\beta) electrons. One can also see the state energy of -2.336035 eV for STATE #1, which is the relative energy with respect to the ROHF reference. The program also lists the summary as:
State Energy Excitation Excitation(eV) <S^2> Transition dipole moment, a.u. Oscillator
Hartree eV rel. GS X Y Z Abs. strength
1 -14.6524451714 -2.336035 0.000000 0.000 0.0000 0.0000 0.0000 0.0000 0.0000
0 -14.5665974741 0.000000 2.336035 (ROHF/UHF Reference state)
2 -14.4805833619 2.340563 4.676598 0.000 -0.0000 -2.0494 -0.4675 2.1021 0.5063
3 -14.4805833619 2.340563 4.676598 0.000 -0.0001 0.4675 -2.0494 2.1021 0.5063
4 -14.4720993070 2.571426 4.907461 0.000 -1.8525 -0.0000 0.0001 1.8525 0.4126
5 -14.3812906885 5.042455 7.378489 0.000 -0.0000 -0.0000 0.0000 0.0000 0.0000
Transition Excitation Transition dipole, a.u. Oscillator
eV x y z Abs. strength
1 -> 2 4.676598 -0.0000 -2.0494 -0.4675 2.1021 0.5063
1 -> 3 4.676598 -0.0001 0.4675 -2.0494 2.1021 0.5063
1 -> 4 4.907461 -1.8525 -0.0000 0.0001 1.8525 0.4126
1 -> 5 7.378489 -0.0000 -0.0000 0.0000 0.0000 0.0000
2 -> 3 0.000000 0.0000 -0.0000 -0.0000 0.0000 0.0000
2 -> 4 0.230863 -0.0000 -0.0000 -0.0000 0.0000 0.0000
2 -> 5 2.701891 -0.3273 0.0000 0.0000 0.3273 0.0071
3 -> 4 0.230863 0.0000 0.0000 -0.0000 0.0000 0.0000
3 -> 5 2.701891 -1.4347 -0.0000 0.0001 1.4347 0.1363
4 -> 5 2.471028 -0.0001 0.0000 -1.4704 1.4704 0.1309
STATE #0 is the ROHF reference. The response states are from STATE #1. In the above case, STATE #1 has a negative energy of -2.336035 as compared to the reference ROHF energy. The first response state, STATE #1 corresponds to the ground singlet state, which is typically lower than the first triplet state. Therefore, the negative energy is not unusual. When you report the relative vertical excitation energy (VEE), you should use the lowest singlet state (STATE #1) as your reference state, not the reference ROHF (STATE #0). The program also reports the corresponding (\langle S^2 \rangle), transition dipole, and oscillator strengths of each response state. The (\langle S^2 \rangle) value indicates the spin state (0 = singlet, 2 = triplet states), while the oscillator strengths indicate the strength of absorption. As you can see, the SELECTING EXCITED STATE IROOT= 1 indicates that STATE #1, which is the ground singlet state, is chosen for further calculations such as geometry optimizations.
The results of singlets and triplets by various methods are summarized in the table below. The MRSF-TDDFT corresponds to the current example. The ( ^1 S ) is STATE #1 above, which serves as the reference energy for all the other states. The ( ^3 P_z, ^3 P_x, ^3 P_y
) should be degenerate. However, they are not exactly degenerate, since the orbital optimization step by reference ROHF breaks the symmetry among them by selectively choosing one particular orbital out of three (p) orbitals. As a result, MRSF-TDDFT produces 2.900 and 2.667 eV as compared to the ( ^1 S ) ground singlet state. Likewise, MRSF-TDDFT produces 4.913 and 4.690 eV VEEs of ( ^1 P_z, ^1 P_x, ^1 P_y ) singlet states. The (\langle S^2 \rangle) values are presented in parentheses. As compared to MRSF-TDDFT, SF-TDDFT is missing one degenerate state of ( ^1 P_{x,y} ). In fact, SF-TDDFT mixes it with ( ^3 P_{x,y} ), producing a half-half mixture of singlet and triplet states with the averaged singlet (2.667 eV of MRSF-TDDFT) and triplet (4.690 eV of MRSF-TDDFT) energy of 3.688 eV. This half-half mixture can be easily seen from its (\langle S^2 \rangle) value of 1.00 in parentheses. The (\langle S^2 \rangle = 1.00) represents neither a singlet nor a triplet state.
