Examples: first-order nonadiabatic couplings

First-order nonadiabatic couplings

The calculations of first-order nonadiabatic couplings (NAC) between ground and excited-states (<0|Dx|Sn>), and those between the excited-states (<Sm|Dx|Sn> or <Tm|Dx|Tn>) at the TD-DFT/TDA level can be achieved by generalizing the standard linear and quadratic response theories, for details, see

1. Zhendong Li and Wenjian Liu, "First-order nonadiabatic coupling matrix elements between excited states: A Lagrangian formulation at the CIS, RPA, TD-HF, and TD-DFT levels", J. Chem. Phys. 141, 014110 (2014).

2. Zhendong Li, Bingbing Suo, and Wenjian Liu, "First-order nonadiabatic coupling matrix elements between excited states: II. Implementation and applications at the TD-DFT and pp-RPA levels", J. Chem. Phys. 141, 244105 (2014).

For convenience, however, in the input they are both specified by the QUAD keyword with single and double, respectively. Either analytic derivative or finite difference approach can be used. The latter is only allowed for C(1) symmetry, and for molecules without orbital degeneracy!

Compared with Hellman-Feynman results

By using the turn over rule, the NAC between two adiabatic states can be written as a integration of the gradients of V_nuc and the transition density scaled by the energy difference. We may call it Hellman-Feynman like expression.

However, for real application, neither the state is exact nor the basis is complete, so this form will not yield accurate results !!! As a simple illustration, we consider the NAC between two triplet states of MgH2 calculated at the pp-TDA level using several different basis sets.

pp-tda for MgH2:

$COMPASS
Title
 nh3
Basis
 aug-cc-pvqz
uncontracted
# nac-mg
#primitives
Geometry
 Mg 0. 0. 5.0
 H 0. 1.25 0. 
 H 0. -1.25 0.
End geometry
units
bohr
skeleton
$END

$xuanyuan
direct
schwarz
$end

$scf
RHF
charge
2
spin
1
THRESHCONV
1.d-12 1.d-10
OPTSCR
1
$end

$tddft
imethod
4
isf
1
nexit
0 0 2 0
itda
1
ipprpa
1
idiag
1
istore
1
crit_e
1.d-12
crit_vec
1.d-10
$end

$resp
iprt
2
QUAD
FNAC
norder
1
method
2
nfiles
1
$end

BASIS SET 1: nac-mgh2 (only containing s/p)

  No. Pair   ExSym   ExEnergies      f     D<S^2>          Dominant Excitations             IPA   Ova     En-E1

    1  B1    1  B1  -18.6728 eV   0.0000   0.0000  92.8%  VV(1):  B1(   2 )->  A1(   5 ) -23.956 0.000    0.0000
    2  B1    2  B1  -16.3909 eV   0.0000   0.0000  88.4%  VV(1):  B1(   3 )->  A1(   5 ) -21.065 0.000    2.2819

...

  Gradient contribution from Final-NAC(R)-HFey
     1       -0.0000000000       -0.0000000000       -0.7243310867
     2        0.3451224889        0.0000000000        0.0561331731
     3       -0.3451224889       -0.0000000000        0.0561331731
  Sum of gradient contribution from Final-NAC(R)-HFey
             -0.0000000000       -0.0000000000       -0.6120647404

...

  Gradient contribution from Final-NAC(R)-Escaled
     1       -0.0000000000       -0.0000000000        0.3549353015
     2        0.7989929157        0.0000000000       -0.1473168767
     3       -0.7989929157       -0.0000000000       -0.1473168774
  Sum of gradient contribution from Final-NAC(R)-Escaled
             -0.0000000000       -0.0000000000        0.0603015474

BASIS SET 2: aug-cc-pVTZ

  No. Pair   ExSym   ExEnergies      f     D<S^2>          Dominant Excitations             IPA   Ova     En-E1

    1  B1    1  B1  -18.6736 eV   0.0000   0.0000  90.6%  VV(1):  B1(   2 )->  A1(   5 ) -23.766 0.000    0.0000
    2  B1    2  B1  -16.9243 eV   0.0000   0.0000  87.8%  VV(1):  B1(   3 )->  A1(   5 ) -21.638 0.000    1.7493

...

  Gradient contribution from Final-NAC(R)-HFey
     1       -0.0000000000       -0.0000000000       -0.2135694894
     2        0.3992751267       -0.0000000000        0.0624241280
     3       -0.3992751267       -0.0000000000        0.0624241280
  Sum of gradient contribution from Final-NAC(R)-HFey
             -0.0000000000       -0.0000000000       -0.0887212334

...

  Gradient contribution from Final-NAC(R)-Escaled
     1        0.0000000000       -0.0000000000        0.3849500968
     2        0.8010940332        0.0000000000       -0.1247569146
     3       -0.8010940332       -0.0000000000       -0.1247570212
  Sum of gradient contribution from Final-NAC(R)-Escaled
              0.0000000000       -0.0000000000        0.1354361609

BASIS SET 3: uncontracted aug-cc-pVQZ

  No. Pair   ExSym   ExEnergies      f     D<S^2>          Dominant Excitations             IPA   Ova     En-E1
 
    1  B1    1  B1  -18.6693 eV   0.0000   0.0000  90.6%  VV(1):  B1(   2 )->  A1(   5 ) -23.760 0.000    0.0000
    2  B1    2  B1  -16.9256 eV   0.0000   0.0000  87.5%  VV(1):  B1(   3 )->  A1(   5 ) -21.636 0.000    1.7438

 ...

  Gradient contribution from Final-NAC(R)-HFey
     1       -0.0000000000        0.0000000000        0.1886204721
     2       -0.4020811753        0.0000000000       -0.0612611462
     3        0.4020811753        0.0000000000       -0.0612611462
  Sum of gradient contribution from Final-NAC(R)-HFey
             -0.0000000000        0.0000000000        0.0660981797

...

  Gradient contribution from Final-NAC(R)-Escaled
     1       -0.0000000000        0.0000000000       -0.3856366870
     2       -0.8006137635        0.0000000000        0.1251752431
     3        0.8006137635       -0.0000000000        0.1251750107
  Sum of gradient contribution from Final-NAC(R)-Escaled
              0.0000000000        0.0000000000       -0.1352864331

It is seen that the results [NAC-HFey] are quite differenct, while those [NAC(R)-Escaled] consider the finite basis set effects (Pulay terms) are similar. While for small molecules like H2, such difference seems to be not significant, however, for large molecules the difference can be very large. (This might be due to the fact that for small systems, the density is small such that the contraction with integrals is also not large).

The same conclusion can also be found for TD-DFT, as already demonstrated by the paper of Send and Furche (see TABLE IV), where u-aug-TZVPP basis set have been used to get qualitatively correct results for the coupling between the ground and first singlet excited state of Cinnoline (C8H6N2) using the PBE functional !

Analytic derivative approach

NAC between two 1A2 states of CH2O:

$COMPASS
Title
 ch2o
Basis
 6-31GP
Geometry
 C                  0.00000000   -0.00000000   -0.53964037
 O                  0.00000000    0.00000000    0.68767663
 H                  0.00000000    0.93940400   -1.13178537
 H                  0.00000000   -0.93940400   -1.13178537
End geometry
skeleton
$END

$xuanyuan
direct
schwarz
$end

$scf
RHF
charge
0
spin
1
THRESHCONV
1.d-10 1.d-8
OPTSCR
1
$end

$tddft
imethod
1
isf
0
nexit
0 2 0 0
itda
0
idiag
1
istore
1
crit_e
1.d-10
crit_vec
1.d-8
lefteig
DirectGrid
$end

$resp
iprt
1
QUAD
FNAC
double
pairs
1
1 2 1 1 2 2
norder
1
method
2
nfiles
1
ignore
1
noresp
$end

Finite difference approach

To use the finite difference approach, nosym must be used to avoid the rotation of molecules. Currently, only C(1) group is permitted. For the example considered above, it turns out that the first two 1A2 excited states are the 1st and 4th excited states. However, if only use iexit=4 with iterative diagonalization, the initialization based on orbital energy difference (IPA) will miss the 2rd excited states, as its IPA is very large, viz.,

  No. Pair   ExSym   ExEnergies      f     D<S^2>          Dominant Excitations             IPA   Ova     En-E1
    1   A    2   A    4.3504 eV   0.0000   0.0000  88.9%  CV(0):   A(   8 )->   A(   9 )  15.558 0.476    0.0000
    2   A    3   A    9.3394 eV   0.0008   0.0000  92.4%  CV(0):   A(   6 )->   A(   9 )  21.045 0.546    4.9890
    3   A    4   A    9.3480 eV   0.1782   0.0000  85.8%  CV(0):   A(   7 )->   A(   9 )  17.814 0.822    4.9975
    4   A    5   A   11.3372 eV   0.0000   0.0000  92.3%  CV(0):   A(   5 )->   A(   9 )  22.197 0.535    6.9868
    5   A    6   A   11.6654 eV   0.3266   0.0000  92.9%  CV(0):   A(   8 )->   A(  10 )  18.527 0.461    7.3150

Thus, iexit=5 is used in the following inputs:

$COMPASS
Title
 ch2o
Basis
 6-31GP
Geometry
 C                  0.00000000   -0.00000000   -0.53964037
 O                  0.00000000    0.00000000    0.68767663
 H                  0.00000000    0.93940400   -1.13178537
 H                  0.00000000   -0.93940400   -1.13178537
End geometry
skeleton
group
c(1)
nosym
$END

$xuanyuan
direct
schwarz
$end

$scf
RHF
charge
0
spin
1
THRESHCONV
1.d-10 1.d-8
OPTSCR
1
$end

$tddft
imethod
1
isf
0
iexit
5
itda
0
idiag
1
istore
1
crit_e
1.d-10
crit_vec
1.d-8
lefteig
DirectGrid
$end

$resp
iprt
1
QUAD
FNAC
double
pairs
1
1 1 1 1 1 4
norder
1
method
2
nfiles
1
ignore
1
noresp
$end

The corresponding input for the finite difference approach is to add two keywords: FDIF for specification and step followed by a real number for step size in finite difference. If the unit is bohr in the COMPASS part, the 'BOHR' should be added in the RESP part also.

$COMPASS
Title
 ch2o
Basis
 6-31GP
Geometry
 C                  0.00000000   -0.00000000   -0.53964037
 O                  0.00000000    0.00000000    0.68767663
 H                  0.00000000    0.93940400   -1.13178537
 H                  0.00000000   -0.93940400   -1.13178537
End geometry
skeleton
group
c(1)
nosym
$END

$xuanyuan
direct
schwarz
$end

$scf
RHF
charge
0
spin
1
THRESHCONV
1.d-10 1.d-8
OPTSCR
1
iaufbau
0
$end

$tddft
imethod
1
isf
0
iexit
5
itda
0
idiag
1
istore
1
crit_e
1.d-10
crit_vec
1.d-8
lefteig
AOKXC
DirectGrid
$end

$resp
iprt
1
QUAD
FNAC
#single
#states
#1
#1 1 2
double
pairs
1
1 1 1 1 1 4
norder
1
method
2
nfiles
1
FDIF
step
0.001
ignore
1
noresp
$end

To use finite-difference, a script fdiff.py along with fdmol.py in bdf-pkg/source/tools/fdiff should be used as

./fbdiff.py run.sh xxx.inp > log

After the calculation is done, an output file xxx.out will present in the current directory. The log file saves the information during the calculations.

Illustration

Visualisation of the results can be achieved via the tools in bdf-pkg/source/tools/fdiff/NACplot.nb.

Static: t1_nac.PNG

Dynamic: t1_nac.GIF

Examples: first-order nonadiabatic couplings (last edited 2018-12-11 13:55:27 by lzd)