Q-Chem 5.1 User’s Manual

10.6 Nonadiabatic Couplings and Optimization of Minimum-Energy Crossing Points

10.6.1 Nonadiabatic Couplings

Conical intersections are degeneracies between Born-Oppenheimer potential energy surfaces that facilitate non-adiabatic transitions between excited states, i.e., internal conversion and intersystem crossing processes, both of which represent a breakdown of the Born-Oppenheimer approximation.[Matsika and Krause(2011), Herbert et al.(2016)Herbert, Zhang, Morrison, and Liu] Although simultaneous intersections between more than two electronic states are possible,[Matsika and Krause(2011)] consider for convenience the two-state case, and let

  \begin{equation} \label{eq:H_2x2} \mathbf{H} = \left( \begin{array}{cc} H^{}_{JJ}(\mathbf{R}) & H^{}_{JK}(\mathbf{R}) \\ H_{KJ}^\ast (\mathbf{R}) & H^{}_{KK}(\mathbf{R}) \\ \end{array} \right) \; . \end{equation}   (10.2)

denote the matrix representation of the vibronic (vibrational + electronic) Hamiltonian [Eq. eq:Hamiltonian] in a basis of two electronic states, $J$ and $K$. (Electronic degrees of freedom have been integrated out of this expression, and $\mathbf{R}$ represents the remaining, nuclear coordinates.) By definition, the Born-Oppenheimer states are the ones that diagonalize $\mathbf{H}$ at a particular molecular geometry $\mathbf{R}$, and thus two conditions must be satisfied in order to obtain degeneracy in the Born-Oppenheimer representation: $H^{}_{JJ} = H^{}_{KK}$ and $H^{}_{JK}=0$. As such, degeneracies between two Born-Oppenheimer potential energy surfaces exist in subspaces of dimension $N_{\rm int}-2$, where $N_{\rm int} = 3 N_{\rm atoms}-6$ is the number of internal (vibrational) degrees of freedom (assuming the molecule is non-linear). This ($N_{\rm int}-2$)-dimensional subspace is known as the seam space because the two states are degenerate everywhere within this space. In the remaining two degrees of freedom, known as the branching space, the degeneracy between Born-Oppenheimer surfaces is lifted by an infinitesimal displacement, which in a three-dimensional plot resembles a double cone about the point of intersection, hence the name conical intersection.

The branching space is defined by the span of a pair of vectors $\mathbf{g}^{}_{JK}$ and $\mathbf{h}^{}_{JK}$. The former is simply the difference in the gradient vectors of the two states in question,

  \begin{equation} \label{eq:gJK} \mathbf{g}^{}_{JK} = \frac{\partial E_ J}{\partial \textbf{R}} - \frac{\partial E_ K}{\partial \textbf{R}} \; , \end{equation}   (10.3)

and is readily evaluated at any level of theory for which analytic energy gradients are available (or less-readily, via finite difference, if they are not!). The definition of the non-adiabatic coupling vector $\mathbf{h}^{}_{JK}$, on the other hand, is more involved and not directly amenable to finite-difference calculations:

  \begin{equation} \label{eq:hJK} \mathbf{h}^{}_{JK} = \left\langle \Psi _ J \left| \left( \partial \hat{H}/\partial \mathbf{R} \right)\right|\Psi _ K\right\rangle \; . \end{equation}   (10.4)

This is closely related to the derivative coupling vector

  \begin{equation} \label{eq:dJK} \mathbf{d}^{}_{JK} = \bigl \langle \Psi _ J \bigl | \bigl ( \partial /\partial \mathbf{R} \bigr )\bigr |\Psi _ K\bigr \rangle = \frac{\mathbf{h}^{}_{JK}}{E_ J-E_ K} \; . \end{equation}   (10.5)

The latter expression for $\mathbf{d}^{}_{JK}$ demonstrates that the coupling between states becomes large in regions of the potential surface where the two states are nearly degenerate. The relative orientation and magnitudes of the vectors $\mathbf{g}_{JK}^{}$ and $\mathbf{h}^{}_{JK}$ determined the topography around the intersection, i.e., whether the intersection is “peaked” or “sloped”; see Ref. Herbert:2016b for a pedagogical overview.

Algorithms to compute the non-adiabatic couplings $\mathbf{d}^{}_{JK}$ are not widely available in quantum chemistry codes, but thanks to the efforts of the Herbert and Subotnik groups, they are available in Q-Chem when the wave functions $\Psi _ J$ and $\Psi _ K$, and corresponding electronic energies $E_ J$ and $E_ K$, are computed at the CIS or TDDFT level,[Fatehi et al.(2011)Fatehi, Alguire, Shao, and Subotnik, Zhang and Herbert(2014b), Zhang and Herbert(2015a), Ou et al.(2015)Ou, Bellchambers, Furche, and Subotnik] or at the corresponding spin-flip (SF) levels of theory (SF-CIS or SF-TDDFT). The spin-flip implementation[Zhang and Herbert(2014b)] is particularly significant, because only that approach—and not traditional spin-conserving CIS or TDDFT—affords correct topology around conical intersections that involve the ground state. To understand why, suppose that $J$ in Eq. eq:H_2x2 represents the ground state; call it $J=0$ for definiteness. In linear response theory (TDDFT) or in CIS (by virtue of Brillouin’s theorem), the coupling matrix elements between the reference (ground) state and all of the excited states vanish identically, hence $H_{0K}(\mathbf{R})\equiv 0$. This means that there is only one condition to satisfy in order to obtain degeneracy, hence the branching space is one- rather than two-dimensional, for any conical intersection that involves the ground state.[Levine et al.(2006)Levine, Ko, Quenneville, and Mart\'{\i }nez] (For intersections between two excited states, the topology should be correct.) In the spin-flip approach, however, the reference state has a different spin multiplicity than the target states; if the latter have spin quantum number $S$, then the reference state has spin $S+1$. This has the effect that the ground state of interest (spin $S$) is treated as an excitation, and thus on a more equal footing with excited states of the same spin, and it rigorously fixes the topology problem around conical intersections.[Zhang and Herbert(2014b)]

Nonadiabatic (derivative) couplings are available for both CIS and TDDFT. The CIS non-adiabatic couplings can be obtained from direct differentiation of the wave functions with respect to nuclear positions.[Fatehi et al.(2011)Fatehi, Alguire, Shao, and Subotnik, Zhang and Herbert(2014b)] For TDDFT, the same procedure can be carried out to calculate the approximate non-adiabatic couplings, in what has been termed the “pseudo-wave function” approach.[Ou et al.(2014)Ou, Fatehi, Alguire, Shao, and Subotnik, Zhang and Herbert(2014b)] Formally more rigorous TDDFT non-adiabatic couplings derived from quadratic response theory are also available, although they are subject to certain undesirable, accidental singularities if for the two states $J$ and $K$ in Eq. eq:hJK, the energy difference $|E_ J-E_ K|$ is quasi-degenerate with the excitation energy $\omega _ I^{} = E_ I - E_0$ for some third state, $I$.[Zhang and Herbert(2015a), Ou et al.(2015)Ou, Bellchambers, Furche, and Subotnik] As such, the pseudo-wave function method is the recommended approach for computing non-adiabatic couplings with TDDFT, although in the spin flip case the pseudo-wave function approach is rigorously equivalent to the pseudo-wave function approach, and is free of singularities.[Zhang and Herbert(2015a)]

Finally, we note that there is some evidence that SF-TDDFT calculations are most accurate when used with functionals containing $\sim $50% Hartree-Fock exchange,[Shao et al.(2003)Shao, Head-Gordon, and Krylov, Huix-Rotllant et al.(2010)Huix-Rotllant, Natarajan, Ipatov, Wawire, Deutsch, and Casida] and many studies with this method (see Ref. Herbert:2016b for a survey) have used the BH&HLYP functional, in which LYP correlation is combined with Becke’s “half and half” (BH&H) exchange functional, consisting of 50% Hartree-Fock exchange and 50% Becke88 exchange (EXCHANGE = BHHLYP in Q-Chem.)

10.6.2 Job Control and Examples

In order to perform non-adiabatic coupling calculations, the $derivative_coupling section must be given:

$derivative_coupling
one line comment
$i,j,k,...$
$end
Nonadiabatic couplings will then be computed between all pairs of the states $i, j, k, \ldots $; use “0” to request the HF or DFT reference state, “1” for the first excited state, etc.

CALC_NAC

Determines whether we are calculating non-adiabatic couplings.


TYPE:

LOGICAL


DEFAULT:

FALSE


OPTIONS:

TRUE

Calculate non-adiabatic couplings.

FALSE

Do not calculate non-adiabatic couplings.


RECOMMENDATION:

None.


CIS_DER_NUMSTATE

Determines among how many states we calculate non-adiabatic couplings. These states must be specified in the $derivative_coupling section.


TYPE:

INTEGER


DEFAULT:

0


OPTIONS:

0

Do not calculate non-adiabatic couplings.

$n$

Calculate $n(n-1)/2$ pairs of non-adiabatic couplings.


RECOMMENDATION:

None.


SET_QUADRATIC

Determines whether to include full quadratic response contributions for TDDFT.


TYPE:

LOGICAL


DEFAULT:

FALSE


OPTIONS:

TRUE

Include full quadratic response contributions for TDDFT.

FALSE

Use pseudo-wave function approach.


RECOMMENDATION:

The pseudo-wave function approach is usually accurate enough and is free of accidental singularities. Consult Refs. Zhang:2015a and Ou:2015 for additional guidance.


Example 10.222  Nonadiabatic couplings among the lowest five singlet states of ethylene, computed at the TD-B3LYP level using the pseudo-wave function approach.

$molecule
   0 1
   C    1.85082356   -1.78953123    0.00000000
   H    2.38603593   -2.71605577    0.00000000
   H    0.78082359   -1.78977646    0.00000000
   C    2.52815456   -0.61573833    0.00000000
   H    1.99294220    0.31078621    0.00000000
   H    3.59815453   -0.61549310    0.00000000
$end

$rem
   CIS_N_ROOTS        4
   CIS_TRIPLETS       false
   SET_ITER           50
   CIS_DER_NUMSTATE   5
   CALC_NAC           true
   EXCHANGE           b3lyp
   BASIS              6-31G*
$end

$derivative_coupling
   0 is the reference state
   0 1 2 3 4
$end


Example 10.223  Nonadiabatic couplings between $S_0$ and $S_1$ states of ethylene using BH&HLYP and spin-flip TDDFT.

$molecule
   0 3
   C    1.85082356   -1.78953123    0.00000000
   H    2.38603593   -2.71605577    0.00000000
   H    0.78082359   -1.78977646    0.00000000
   C    2.52815456   -0.61573833    0.00000000
   H    1.99294220    0.31078621    0.00000000
   H    3.59815453   -0.61549310    0.00000000
$end

$rem
   SPIN_FLIP          true
   UNRESTRICTED       true
   CIS_N_ROOTS        4
   CIS_TRIPLETS       false
   SET_ITER           50
   CIS_DER_NUMSTATE   2
   CALC_NAC           true
   EXCHANGE           bhhlyp
   BASIS              6-31G*
$end

$derivative_coupling
   comment
   1 3
$end

Example 10.224  Nonadiabatic couplings between $S_1$ and $S_2$ states of ethylene computed via quadratic response theory at the TD-B3LYP level.

$molecule
   0 1
   C    1.85082356   -1.78953123    0.00000000
   H    2.38603593   -2.71605577    0.00000000
   H    0.78082359   -1.78977646    0.00000000
   C    2.52815456   -0.61573833    0.00000000
   H    1.99294220    0.31078621    0.00000000
   H    3.59815453   -0.61549310    0.00000000
$end

$rem
   CIS_N_ROOTS       4
   CIS_TRIPLETS      false
   RPA               true
   SET_ITER          50
   CIS_DER_NUMSTATE  2
   CALC_NAC          true
   EXCHANGE          b3lyp
   BASIS             6-31G*
   SET_QUADRATIC     true #include full quadratic response
$end

$derivative_coupling
   comment
   1 2
$end

10.6.3 Minimum-Energy Crossing Points

The seam space of a conical intersection is really a (hyper)surface of dimension $N_{\rm int}-2$, and while the two electronic states in question are degenerate at every point within this space, the electronic energy varies from one point to the next. To provide a simple picture of photochemical reaction pathways, it is often convenient to locate the minimum-energy crossing point (MECP) within this $(N_{\rm int}-2)$-dimensional seam. Two separate minimum-energy pathway searches, one on the excited state starting from the ground-state geometry and terminating at the MECP, and the other on the ground state starting from the MECP and terminating at the ground-state geometry, then affords a photochemical mechanism. (See Ref. Zhang:2014a for a simple example.) In some sense, then, the MECP is to photochemistry what the transition state is to reactions that occur on a single Born-Oppenheimer potential energy surface. One should be wary of pushing this analogy too far, because whereas a transition state reasonably be considered to be a bottleneck point on the reaction pathway, the path through a conical intersection may be downhill and perhaps therefore more likely to proceed from one surface to the other at a point “near" the intersection, and in addition there can be multiple conical intersections between the same pair of states so more than one photochemical mechanism may be at play. Such complexity could be explored, albeit at significantly increased cost, using non-adiabatic “surface hopping" ab initio molecular dynamics, as described in Section 10.7.6. Here we describe the computationally-simpler procedure of locating an MECP along a conical seam.

Recall that the branching space around a conical intersection between electronic states $J$ and $K$ is spanned by two vectors, $\mathbf{g}^{}_{JK}$ [Eq. eq:gJK] and $\mathbf{h}^{}_{JK}$ [Eq. eq:hJK]. While the former is readily available in analytic form for any electronic structure method that has analytic excited-state gradients, the non-adiabatic coupling vector $\mathbf{h}^{}_{JK}$ is not available for most methods. For this reason, several algorithms have been developed to optimize MECPs without the need to evaluate $\mathbf{h}^{}_{JK}$, and three such algorithms are available in Q-Chem.

Mart\'{\i }nez and coworkers[Levine et al.(2008)Levine, Coe, and Martinez] developed a penalty-constrained MECP optimization algorithm that consists of minimizing the objective function

  \begin{equation} F_\sigma (\mathbf{R}) = \tfrac {1}{2}\bigl [E_ I(\mathbf{R}) + E_ J(\mathbf{R})\bigr ] + \sigma \left(\frac{\bigl [E_ I(\mathbf{R}) - E_ J(\mathbf{R})\bigr ]^2}{E_ I(\mathbf{R}) - E_ J(\mathbf{R}) + \alpha }\right) \; , \end{equation}   (10.6)

where $\alpha $ is a fixed parameter to avoid singularities and $\sigma $ is a Lagrange multiplier for a penalty function meant to drive the energy gap to zero. Minimization of $F_\sigma $ is performed iteratively for increasingly large values $\sigma $.

A second MECP optimization algorithm is a simplification of the penalty-constrained approach that we call the “direct” method. Here, the gradient of the objective function is

  \begin{equation} \label{eq:direct_ MECP} \mathbf{G}=\mathbf{P}\mathbf{G}_{\rm mean}+2(E_ K-E_ J)\mathbf{G}^{}_{\rm diff} \; , \end{equation}   (10.7)

where

  \begin{equation} \mathbf{G}^{}_{\rm mean}=\tfrac {1}{2}(\mathbf{G}_ J+\mathbf{G}_ K) \end{equation}   (10.8)

is the mean energy gradient, with $\mathbf{G}_ i = \partial E_ i/\partial \mathbf{R}$ being the nuclear gradient for state $i$, and

  \begin{equation} \mathbf{G}^{}_{\rm diff} = \frac{\mathbf{G}_ K-\mathbf{G}_ J}{||\mathbf{G}_ K-\mathbf{G}_ J||} \end{equation}   (10.9)

is the normalized difference gradient. Finally,

  \begin{equation} \label{eq:MECP_ projector1} \mathbf{P} = \bm@general \boldmath \m@ne \mv@bold \bm@command {1}-\mathbf{G}^{}_{\rm diff}\mathbf{G}_{\rm diff}^\top \end{equation}   (10.10)

projects the gradient difference direction out of the mean energy gradient in Eq. eq:direct_MECP. The algorithm then consists in minimizing along the gradient $\mathbf{G}$, with for the iterative cycle over a Lagrange multiplier, which can sometimes be slow to converge.

The third and final MECP optimization algorithm that is available in Q-Chem is the branching-plane updating method developed by Morokuma and coworkers[Maeda et al.(2010)Maeda, Ohno, and Morokuma] and implemented in Q-Chem by Zhang and Herbert.[Zhang and Herbert(2014a)] This algorithm uses a gradient that is similar to that in Eq. eq:direct_MECP but projects out not just $\mathbf{G}_{\rm diff}$ in Eq. eq:MECP_projector1 but also a second vector that is orthogonal to it, representing an iteratively-updated approximation to the branching space.

None of these three methods requires evaluation of non-adiabatic couplings, and all three can be used to optimize MECPs at the CIS, SF-CIS, TDDFT, SF-TDDFT, and SOS-CIS(D0) levels. The direct algorithm can also be used for EOM-XX-CCSD methods (XX = EE, IP, or EA). It should be noted that since EOM-XX-CCSD is a linear response method, it suffers from the same topology problem around conical intersections involving the ground state that was described in regards to TDDFT in Section 10.6.1. With spin-flip approaches, correct topology is obtained.[Zhang and Herbert(2014b)]

Analytic derivative couplings are available for (SF-)CIS and (SF-)TDDFT, so for these methods one can alternatively employ an optimization algorithm that makes use of both $\mathbf{g}^{}_{JK}$ and $\mathbf{h}^{}_{JK}$. Such an algorithm, due to Schlegel and coworkers,[Bearpark et al.(1994)Bearpark, Robb, and Schlegel] is available in Q-Chem and consists of optimization along the gradient in Eq. eq:direct_MECP but with a projector

  \begin{equation} \label{eq:MECP_ projector2} \mathbf{P} = \bm@general \boldmath \m@ne \mv@bold \bm@command {1}-\mathbf{G}^{}_{\rm diff} \mathbf{G}_{\rm diff}^\top - \mathbf{y}\mathbf{y}^\top \end{equation}   (10.11)

where

  \begin{equation} \mathbf{y} = \frac{ (\bm@general \boldmath \m@ne \mv@bold \bm@command {1}-\mathbf{x}\mathbf{x}^\top )\mathbf{h}^{}_{JK} }{ || (\bm@general \boldmath \m@ne \mv@bold \bm@command {1}-\mathbf{x}\mathbf{x}^\top )\mathbf{h}^{}_{JK} || } \; , \end{equation}   (10.12)

in place of the projector in Eq. eq:MECP_projector1. Equation eq:MECP_projector2 has the effect of projecting the span of $\mathbf{g}^{}_{JK}$ and $\mathbf{h}^{}_{JK}$ (i.e., the branching space) out of state-averaged gradient in Eq. eq:direct_MECP. The tends to reduce the number of iterations necessary to converge the MECP, and since calculation of the (optional) $\mathbf{h}^{}_{JK}$ vector represents only a slight amount of overhead on top of the (required) $\mathbf{g}^{}_{JK}$ vector, this last algorithm tends to yield significant speed-ups relative to the other three.[Zhang and Herbert(2014b)] As such, it is the recommended choice for (SF-)CIS and (SF-)TDDFT.

It should be noted that while the spin-flip methods cure the topology problem around conical intersections that involve the ground state, this method tends to exacerbate spin contamination relative to the corresponding spin-conserving approaches.[Zhang and Herbert(2015b)] While spin contamination is certainly present in traditional, spin-conserving CIS and TDDFT, it presents the following unique challenge in spin-flip methods. Suppose, for definiteness, that one is interested in singlet excited states. Then the reference state for the spin-flip methods should be the high-spin triplet. A spin-flipping excitation will then generate S$_0$, S$_1$, S$_2,\ldots $ but will also generate the $M_ S=0$ component of the triplet reference state, which therefore appears in what is ostensibly the singlet manifold. Q-Chem attempts to identify this automatically, based on a threshold for $\langle \Hat {S}^2\rangle $, but severe spin contamination can sometimes defeat this algorithm,[Zhang and Herbert(2014a)] hampering Q-Chem’s ability to distinguish singlets from triplets (in this particular example). An alternative might be the state-tracking procedure that is described in Section 10.6.5.

10.6.4 Job Control and Examples

For MECP optimization, set MECP_OPT = TRUE in the $rem section, and note that the $derivative_coupling input section discussed in Section 10.6.2 is not necessary in this case.

MECP_OPT

Determines whether we are doing MECP optimizations.


TYPE:

LOGICAL


DEFAULT:

FALSE


OPTIONS:

TRUE

Do MECP optimization.

FALSE

Do not do MECP optimization.


RECOMMENDATION:

None.


MECP_METHODS

Determines which method to be used.


TYPE:

STRING


DEFAULT:

BRANCHING_PLANE


OPTIONS:

BRANCHING_PLANE

Use the branching-plane updating method.

MECP_DIRECT

Use the direct method.

PENALTY_FUNCTION

Use the penalty-constrained method.


RECOMMENDATION:

The direct method is stable for small molecules or molecules with high symmetry. The branching-plane updating method is more efficient for larger molecules but does not work if the two states have different symmetries. If using the branching-plane updating method, GEOM_OPT_COORDS must be set to 0 in the $rem section, as this algorithm is available in Cartesian coordinates only. The penalty-constrained method converges slowly and is suggested only if other methods fail.


MECP_STATE1

Sets the first Born-Oppenheimer state for MECP optimization.


TYPE:

INTEGER/INTEGER ARRAY


DEFAULT:

None


OPTIONS:

[$i$,$j$]

Find the $j$th excited state with the total spin $i$; $j=0$ means the SCF ground state.


RECOMMENDATION:

$i$ is ignored for restricted calculations; for unrestricted calculations, $i$ can only be 0 or 1.


MECP_STATE2

Sets the second Born-Oppenheimer state for MECP optimization.


TYPE:

INTEGER/INTEGER ARRAY


DEFAULT:

None


OPTIONS:

[$i$,$j$]

Find the $j$th excited state with the total spin $i$; $j=0$ means the SCF ground state.


RECOMMENDATION:

$i$ is ignored for restricted calculations; for unrestricted calculations, $i$ can only be 0 or 1.


CIS_S2_THRESH

Determines whether a state is a singlet or triplet in unrestricted calculations.


TYPE:

INTEGER


DEFAULT:

120


OPTIONS:

$n$

Sets the $\langle \Hat {S}^2\rangle $ threshold to $n/100$


RECOMMENDATION:

For the default case, states with $\langle \Hat {S}^2\rangle > 1.2$ are treated as triplet states and other states are treated as singlets.


MECP_PROJ_HESS

Determines whether to project out the coupling vector from the Hessian when using branching plane updating method.


TYPE:

LOGICAL


DEFAULT:

TRUE


OPTIONS:

TRUE

FALSE


RECOMMENDATION:

Use the default.


Example 10.225  MECP optimization of an intersection between the S$_2$ and S$_3$ states of NO$_2^-$, using the direct method at the SOS-CIS(D0) level.

$molecule
   -1 1
   N1
   O2  N1  RNO
   O3  N1  RNO  O2  AONO
   
   RNO  = 1.50
   AONO = 100
$end

$rem
   JOBTYPE         = opt
   METHOD          = soscis(d0)
   BASIS           = aug-cc-pVDZ
   AUX_BASIS       = rimp2-aug-cc-pVDZ
   PURECART        = 1111
   CIS_N_ROOTS     = 4
   CIS_TRIPLETS    = false
   CIS_SINGLETS    = true
   MEM_STATIC      = 900
   MEM_TOTAL       = 1950
   MECP_OPT        = true
   MECP_STATE1     = [0,2]
   MECP_STATE2     = [0,3]
   MECP_METHODS    = mecp_direct
$end


Example 10.226  Optimization of the ethylidene MECP between S$_0$ and S$_1$ in $\rm C_2H_2$, at the SF-TDDFT level using the branching-plane updating method.

$molecule
   0 3
   C   0.044626  -0.2419240   0.357157
   C   0.008905   0.6727548   1.460500
   H   0.928425  -0.1459163  -0.272095
   H  -0.831032  -0.1926895  -0.288529
   H  -0.009238   0.9611331   2.479936
   H   0.068314  -1.2533580   0.778847
$end

$rem
   JOBTYPE           opt
   MECP_OPT          true
   MECP_METHODS      branching_plane
   MECP_PROJ_HESS    true ! project out y vector from the hessian
   GEOM_OPT_COORDS   0  ! currently only works for Cartesian coordinate
   METHOD            bhhlyp
   SPIN_FLIP         true
   UNRESTRICTED      true
   BASIS             6-31G(d,p)
   CIS_N_ROOTS       4
   MECP_STATE1       [0,1]
   MECP_STATE2       [0,2]
   CIS_S2_THRESH     120
$end

Example 10.227  Optimization of the twisted-pyramidalized ethylene MECP between S$_0$ and S$_1$ in $\rm C_2H_2$ using SF-TDDFT.

$molecule
   0 3
   C  -0.015889   0.073532  -0.059559
   C   0.012427  -0.002468   1.315694
   H   0.857876   0.147014  -0.710529
   H  -0.936470  -0.011696  -0.626761
   H   0.764557   0.663381   1.762573
   H   0.740773  -0.869764   1.328583
$end

$rem
   JOBTYPE         opt
   MECP_OPT        true
   MECP_METHODS    penalty_function
   METHOD          bhhlyp
   SPIN_FLIP       true
   UNRESTRICTED    true
   BASIS           6-31G(d,p)
   CIS_N_ROOTS     4
   MECP_STATE1     [0,1]
   MECP_STATE2     [0,2]
   CIS_S2_THRESH   120
$end


Example 10.228  Optimization of the $\tilde{B}^{1}A_2$ and $\tilde{A}^{1}B_2$ states of N$_3^+$ using the direct method at the EOM-EE-CCSD level.

$molecule
   1 1
   N1
   N2 N1 rNN
   N3 N2 rNN N1 aNNN
   
   rNN  =  1.54
   aNNN = 50.0
$end

$rem
   JOBTYPE                opt
   MECP_OPT               true
   MECP_METHODS           mecp_direct
   METHOD                 eom-ccsd
   BASIS                  6-31g
   EE_SINGLETS            [0,2,0,2]
   XOPT_STATE_1           [0,2,2]
   XOPT_STATE_2           [0,4,1]
   CCMAN2                 false
   GEOM_OPT_TOL_GRADIENT  30
$end

Example 10.229  Optimization of the ethylidene MECP between S$_0$ and S$_1$, using BH&HLYP spin-flip TDDFT with analytic derivative couplings.

$molecule
   0 3
   C   0.044626  -0.241924   0.357157
   C   0.008905   0.672754   1.460500
   H   0.928425  -0.145916  -0.272095
   H  -0.831032  -0.192689  -0.288529
   H  -0.009238   0.961133   2.479936
   H   0.068314  -1.253358   0.778847
$end

$rem
   JOBTYPE            opt
   MECP_OPT           true
   MECP_METHODS       branching_plane
   MECP_PROJ_HESS     true
   GEOM_OPT_COORDS    0
   MECP_STATE1        [0,1]
   MECP_STATE2        [0,2]
   UNRESTRICTED       true
   SPIN_FLIP          true
   CIS_N_ROOTS        4
   CALC_NAC           true
   CIS_DER_NUMSTATE   2
   SET_ITER           50
   EXCHANGE           bhhlyp
   BASIS              6-31G(d,p)
   SYMMETRY_IGNORE    true
$end

10.6.5 State-Tracking Algorithm

For optimizing excited-state geometries and other applications, it can be important to find and follow electronically excited states of a particular character as the geometry changes. Various state-tracking procedures have been proposed for such cases.[Harabuchi et al.(2014)Harabuchi, Keipert, Zahariev, Taketsugu, and Gordon, Zhang and Herbert(2015b)] An excited-state, state-tracking algorithm available in Q-Chem is based on the overlap of the attachment/detachment densities at successive steps (Section 7.12.1).[Closser et al.(2014)Closser, Gessner, and Head-Gordon] Using the densities avoids any issues that may be introduced by sign changes in the orbitals or configuration-interaction coefficients.

Two parameters are used to influence the choice of the electronic surface. One ($\gamma ^{}_ E$) controls the energy window for states included in the search, and the other ($\gamma ^{}_ S$) controls how well the states must overlap in order to be considered of the same character. These can be set by the user or generated automatically based on the magnitude of the nuclear displacement. The energy window is defined relative to the estimated energy for the current step (i.e., $E_{\rm est}\pm \gamma ^{}_ E$), which in turn is based on the energy, gradient and nuclear displacement of previous steps. This estimated energy is specific to the type of calculation (e.g., geometry optimization).

The similarity metric for the overlap is defined as

  \begin{equation} \label{eq:similarity_ metric} \mathcal S = 1-\tfrac {1}{2}\bigl ( ||\Delta \mathbf A|| +|| \Delta \mathbf D||\bigr ) \end{equation}   (10.13)

where $\Delta \mathbf A =\mathbf A_{t+1}-\mathbf A_ t$ is the difference in attachment density matrices (Eq. eq:attach_density_matrix) and $\Delta \mathbf D = \mathbf D_{t+1}-\mathbf D_ t$ is the difference in detachment density matrices (Eq. eq:detach_density_matrix), at successive steps. Equation eq:similarity_metric uses the matrix spectral norm,

  \begin{equation} ||\mathbf{M}|| = \left( \lambda _{\rm max}\mathbf{M}^\dagger \mathbf{M} \right)^{1/2} \end{equation}   (10.14)

where $\lambda _{\rm max}$ is the largest eigenvalue of $\mathbf{M}$.

The selected state always satisfies one of the following

  1. It is the only state in the window defined by $\gamma ^{}_ E$.

  2. It is the state with the largest overlap, provided at least one state has $\mathcal{S} \geq \gamma ^{}_ S$.

  3. It is the nearest state energetically if all states in the window have $\mathcal{S} < \gamma ^{}_ S$, or if there are no states in the energy window.

State-following can currently be used with CIS or TDDFT excited states and is initiated with the $rem variable STATE_FOLLOW. It can be used with geometry optimization, ab initio molecular dynamics,[Closser et al.(2014)Closser, Gessner, and Head-Gordon] or with the freezing/growing-string method. The desired state is specified using SET_STATE_DERIV for optimization or dynamics, or using SET_STATE_REACTANT and SET_STATE_PRODUCT for the freezing- or growing-string methods. The results for geometry optimizations can be affected by the step size (GEOM_OPT_DMAX), and using a step size smaller than the default value can provide better results. Also, it is often challenging to converge the strings in freezing/growing-string calculations.

STATE_FOLLOW

Turns on state following.


TYPE:

LOGICAL


DEFAULT:

FALSE


OPTIONS:

FALSE

Do not use state-following.

TRUE

Use state-following.


RECOMMENDATION:

None.


FOLLOW_ENERGY

Adjusts the energy window for near states


TYPE:

INTEGER


DEFAULT:

0


OPTIONS:

0

Use dynamic thresholds, based on energy difference between steps.

$n$

Search over selected state $E_{\rm est}\pm n\times 10^{-6}\; E_ h$.


RECOMMENDATION:

Use a wider energy window to follow a state diabatically, smaller window to remain on the adiabatic state most of the time.


FOLLOW_OVERLAP

Adjusts the threshold for states of similar character.


TYPE:

INTEGER


DEFAULT:

0


OPTIONS:

0

Use dynamic thresholds, based on energy difference between steps.

$n$

Percentage overlap for previous step and current step.


RECOMMENDATION:

Use a higher value to require states have higher degree of similarity to be considered the same (more often selected based on energy).