Analytical derivative coupling for multistate CASPT2 theory

TL;DR

This paper develops an analytical method for derivative couplings in multistate CASPT2 theory, enabling efficient and accurate simulations of photochemical processes involving non-radiative transitions.

Contribution

It introduces an algorithm for analytical evaluation of derivative couplings in MS-CASPT2 and XMS-CASPT2, matching the computational cost of energy gradient calculations.

Findings

Geometries and energies at MECIs agree with multireference CI results.

The method successfully optimized MECIs for stilbene and GFP chromophore.

Computational efficiency is comparable to nuclear energy gradient calculations.

Abstract

The probability of non-radiative transitions in photochemical dynamics is determined by the derivative couplings, the couplings between different electronic states through the nuclear degrees of freedom. Efficient and accurate evaluation of the derivative couplings is, therefore, of central importance to realize reliable computer simulations of photochemical reactions. In this work, the derivative couplings for multistate multireference second-order perturbation theory (MS-CASPT2) and its 'extended' variant (XMS-CASPT2) are studied, in which we present an algorithm for their analytical evaluation. The computational costs for evaluating the derivative couplings are essentially the same as those for calculating the nuclear energy gradients. The geometries and energies calculated with XMS-CASPT2 for small molecules at minimum energy conical intersections (MECIs) are in good agreement with those computed by multireference configuration interaction. As numerical examples, MECIs are optimized using XMS-CASPT2 for stilbene and a GFP model chromophore (the 4-para-hydroxybenzylidene-1,2-dimethyl-imidazolin-5-one anion).

Tables

Table 1. Table 1: Energies of the ground and excited states of ethylene at the S 0 subscript S 0 \mathrm{S}_{0} minimum geometry and the MECIs relative to the S 0 subscript S 0 \mathrm{S}_{0} minimum energy (eV). The aug-cc-pVDZ basis set and CAS(2 e ,2 o ) are used.

	XMS-CASPT2			unc-MRCI+Qa
	$S_0$	$S_1$	$S_2$	$S_0$	$S_1$	$S_2$
$S_0$	0.00	7.70	13.04	0.00	7.79	…
Py	4.42	4.42	9.05	4.40	4.46	9.24
Et	4.63	4.63	5.79	4.56	4.56	5.59
$C_{3v}$-Et	4.71	4.71	5.78	4.69	4.69	5.59
Hm	5.06	5.06	8.25	5.15	5.20	8.48
To	2.99	5.43	5.43	2.92	5.30	5.42

Equations

$$|L⟩=\sum _I{c}_{I,L}|I⟩,$$

$$|\tilde{M}⟩=\sum _L|L⟩U_{LM},$$

$$\sum _{KL}{U}_{KM}⟨K|\hat{f}|L⟩U_{LN}=E_M^{(0)}{\delta }_{MN}.$$

$${\tilde{c}}_{I,M}=\sum _L{c}_{I,L}{U}_{LM}.$$

$$|{\Phi }_M⟩=|\tilde{M}⟩+\sum _{N\Omega }{\hat{E}}_{\Omega }|\tilde{N}⟩T_{\Omega ,LN},$$

$${\hat{T}}_{MN}=\sum _{\Omega }{\hat{E}}_{\Omega }{T}_{\Omega ,MN},$$

$$|{\Phi }_M^{(1)}⟩=\sum _N{\hat{T}}_{MN}|\tilde{N}⟩.$$

$$⟨\tilde{M}|{\hat{E}}_{\Omega }^{†}(\hat{f}-E_N^{(0)}+E_{\text{shift}})|{\Phi }_N^{(1)}⟩+⟨\tilde{M}|{\hat{E}}_{\Omega }^{†}\hat{H}|\tilde{N}⟩=0,$$

$$H_{LM}^{\text{eff}}=$$

$$-$$

$$\sum _M{H}_{LM}^{\text{eff}}{R}_{MP}=R_{LP}{E}_P^{\text{XMS}},$$

$$|{\Psi }_P⟩=\sum _M|{\Phi }_M⟩R_{MP}.$$

$$|P⟩=\sum _M|\tilde{M}⟩R_{MP},$$

$$\mathcal{L}=⟨{\Psi }_{\mathrm{tot}}|\hat{\mathcal{H}}|{\Psi }_{\mathrm{tot}}⟩-\lambda \left[⟨{\Psi }_{\mathrm{tot}}|{\Psi }_{\mathrm{tot}}⟩-1\right],$$

$$\hat{\mathcal{H}}\left(\mathbf{r},\mathbf{X}\right)={\hat{T}}_{\mathrm{nuc}}\left(\mathbf{X}\right)+{\hat{H}}_{\mathrm{el}}\left(\mathbf{r},\mathbf{X}\right),$$

$${\hat{T}}_{\mathrm{nuc}}=-\sum _A\frac{1}{2M_A}\frac{d^2}{d{\mathbf{X}}_A^2},$$

$$\mathcal{E}=⟨{\Psi }_{\mathrm{tot}}|\hat{\mathcal{H}}|{\Psi }_{\mathrm{tot}}⟩.$$

$${\mathcal{L}}_{\mathrm{BO},\mathbf{X}}$$

$$-\lambda \left[⟨{\Psi }_{\mathrm{el}}(\mathbf{r};\mathbf{X})|{\Psi }_{\mathrm{el}}(\mathbf{r};\mathbf{X})⟩-1\right],$$

$$\frac{\partial {\mathcal{L}}_{\mathrm{BO},\mathbf{X}}}{\partial {\Psi }_{\mathrm{el}}(\mathbf{r};\mathbf{X})}=0,$$

$$|{\Psi }_{\mathrm{tot}}(\mathbf{r},\mathbf{X})⟩=\sum _P|{\Psi }_{\mathrm{el},P}(\mathbf{r};\mathbf{X})⟩{\chi }_P(\mathbf{X}),$$

$$\mathcal{E}=\sum _P{\chi }_P^2\left(T_{\mathrm{nuc},P}+E_P\right)+\sum _{Q\ne P}{\chi }_Q⟨{\Psi }_{\mathrm{el},Q}|{\hat{T}}_{\mathrm{nuc}}|{\Psi }_{\mathrm{el},P}⟩{\chi }_P,$$

$$⟨{\Psi }_{\mathrm{el},Q}|{\hat{T}}_{\mathrm{nuc}}|{\Psi }_{\mathrm{el},P}⟩$$

$$-\frac{1}{2}\sum _A\frac{1}{M_A}\text{innerproduct}{\Psi }_{\mathrm{el},Q}\frac{d^2{\Psi }_{\mathrm{el},P}}{d{\mathbf{X}}_A^2},$$

$${\mathbf{d}}^{QP}=\text{innerproduct}{\Psi }_{\mathrm{el},Q}\frac{d{\Psi }_{\mathrm{el},P}}{d\mathbf{X}}.$$

$${\mathcal{L}}^{(0)}=⟨{\Psi }_{\mathrm{tot}}^{(0)}|\hat{\mathcal{H}}|{\Psi }_{\mathrm{tot}}^{(0)}⟩-\lambda \left[⟨{\Psi }_{\mathrm{tot}}^{(0)}|{\Psi }_{\mathrm{tot}}^{(0)}⟩-1\right],$$

$${\mathcal{L}}^{(2)}=⟨{\Psi }_{\mathrm{tot}}^{(1)}|{\hat{\mathcal{H}}}^{(0)}-{\mathcal{E}}^{(0)}|{\Psi }_{\mathrm{tot}}^{(1)}⟩+2⟨{\Psi }_{\mathrm{tot}}^{(1)}|\hat{\mathcal{H}}|{\Psi }_{\mathrm{tot}}^{(0)}⟩,$$

$${\hat{\mathcal{H}}}^{(0)}={\hat{T}}_{\mathrm{nuc}}+{\hat{H}}_{\mathrm{el}}^{(0)},$$

$${\mathcal{E}}^{(0)}=⟨{\Psi }_{\mathrm{tot}}^{(0)}|{\hat{\mathcal{H}}}^{(0)}|{\Psi }_{\mathrm{tot}}^{(0)}⟩.$$

$${\mathcal{E}}^{PT2}=⟨{\Psi }_{\mathrm{tot}}^{(0)}|\hat{\mathcal{H}}|{\Psi }_{\mathrm{tot}}^{(0)}⟩+⟨{\Psi }_{\mathrm{tot}}^{(1)}|\hat{\mathcal{H}}|{\Psi }_{\mathrm{tot}}^{(0)}⟩.$$

$${\mathcal{L}}_{\mathrm{BO},\mathbf{X}}^{(0)}$$

$$-\lambda \left[⟨{\Psi }_{\mathrm{el}}^{(0)}(\mathbf{r};\mathbf{X})|{\Psi }_{\mathrm{el}}^{(0)}(\mathbf{r};\mathbf{X})⟩-1\right],$$

$${\mathcal{L}}_{\mathrm{BO},\mathbf{X}}^{(2)}$$

$$+2⟨{\Psi }_{\mathrm{el}}^{(1)}(\mathbf{r};\mathbf{X})|{\hat{H}}_{\mathrm{el}}|{\Psi }_{\mathrm{el}}^{(0)}(\mathbf{r};\mathbf{X})⟩,$$

$$\frac{\partial {\mathcal{L}}_{\mathrm{BO},\mathbf{X}}^{(0)}}{\partial {\Psi }_{\mathrm{el}}^{(0)}(\mathbf{r};\mathbf{X})}=0,$$

Full text

I Introduction

Understanding the interaction between molecules and light is an important challenge, not only in basic science but also for technological developments, because it could lead to the efficient utilization of light in photo-functional materials. When molecules are irradiated by photons, the molecules undergo various photochemical processes to relax from their electronic excited states.Valeur (2002) Non-radiative deactivation is an example of a process that plays a vital role in photo-induced structural changes of molecules that are used as photochromic and photomechanical materials.Levine and Martínez (2007); Erbas-Cakmak et al. (2015); Matsika and Krause (2011); Matsika et al. (2014) Non-radiative transitions also act as a competitive deactivation pathway in light emission devices,Olsen et al. (2010) reducing the quantum yield of emission.

These non-radiative transitions are induced by the derivative couplings [also often referred to as nonadiabatic coupling matrix elements (NACMEs)],Lengsfield et al. (1984); Lengsfield and Yarkony (1992); Galvan et al. (2016); Lischka et al. (2004); Barbatti and Lischka (2008); Tajti and Szalay (2009); Christiansen (1999); Ichino et al. (2009); Fatehi et al. (2011); Chernyak and Mukamel (2000); Send and Furche (2010); Zhang and Herbert (2014, 2015); Herbert et al. (2016) which are the couplings between the electronic and nuclear degrees of freedom.111The terms ’derivative coupling’ and ’nonadiabatic coupling’ are conventionally used interchangeably, except for some recent reports, such as Ref. 9 We will mathematically define the derivative couplings in the following. Efficient computation of derivative couplings together with nuclear energy gradients enables on-the-fly dynamics simulations of photochemical processes.Tully (1990); Levine and Martínez (2007); Nelson et al. (2014); Tavernelli (2015) It also allows for locating conical intersections between potential energy surfaces,Domcke et al. (2004); Baer (2006); Bearpark et al. (1994); Manaa and Yarkony (1993); Yarkony (1996); Matsika and Krause (2011) which are the set of geometries where two or more potential energy surfaces intersect with each other. Since the computational costs of these applications are strongly dominated by the underlying computation of the derivative couplings and nuclear energy gradients, development of quantum chemical approaches for their efficient and accurate evaluation has the potential to significantly advance what is considered to be the state of the art method in computational photophysics and photochemistry.

To achieve efficient evaluation of derivative couplings, analytical differentiation techniques have been explored in the last few decades. Historically, analytical derivative couplings for multi-configurational methods were first studied [such as the state-averaged complete active space self-consistent field (SA-CASSCF)Lengsfield et al. (1984); Lengsfield and Yarkony (1992); Galvan et al. (2016) and uncontracted multireference configuration interaction (unc-MRCI) methods.Lengsfield et al. (1984); Lischka et al. (2004)] These two models nevertheless have disadvantages: SA-CASSCF is essentially a mean-field model with static correlation treatment, and does not describe dynamical correlation; unc-MRCI is computational very demanding and is often used without double excitations (i.e., MRCIS that does not describe dynamical correlation).Barbatti and Lischka (2008)

More recently, the analytical evaluation of derivative couplings for single-reference theories has been extensively investigated, including those based on equation-of-motion coupled-cluster theory (EOM-CC),Tajti and Szalay (2009); Christiansen (1999); Ichino et al. (2009) configuration interactions singles (CIS),Fatehi et al. (2011) time-dependent density functional theories (TDDFT),Chernyak and Mukamel (2000); Send and Furche (2010) and their spin-flip variants.Zhang and Herbert (2014, 2015); Herbert et al. (2016) Standard single-reference methods are, however, known to incorrectly predict the dimensionality of the conical intersection spaces between the ground and excited states, because they do not compute the states on an equal footing.Levine et al. (2006); Huix-Rotllant et al. (2013); Li et al. (2014); Tuna et al. (2015); Herbert et al. (2016) There have been attempts to resolve this problem.Li et al. (2014); Huix-Rotllant et al. (2013) In addition, spin-flip single-reference methods that treat the ground and excited states equally have been successfully applied to the computation of the electronic structure around conical intersections,Epifanovsky and Krylov (2007); Shepler et al. (2008); Minezawa and Gordon (2009); Zhang and Herbert (2015); Herbert et al. (2016) yielding conical intersection structures that are in agreement with those obtained with CASSCF or MRCI. There have nevertheless been discrepancies in the energies of the conical intersections between SF-TDDFT and CASSCF computations.Herbert et al. (2016) To assess these models, especially for large molecular systems, the development of efficient multi-reference electron correlation methods for calculating derivative couplings is warranted.

We therefore turn our attention to one of the most successful multi-reference models, complete active space second order perturbation theory (CASPT2).Andersson et al. (1990, 1992) The CASPT2 method is a post-CASSCF method that describes dynamical correlation up to the second order. It uses so-called fully internally contracted basis functions to expand the first-order wave functions for efficiency. The electronic structure around conical intersections can be accurately described by its multistate variant (MS-CASPT2),Finley et al. (1998) which diagonalizes an effective Hamiltonian formed from the state-specific CASPT2 wave functions. The MS-CASPT2 method has subsequently been improved by the ‘extended’ variant (XMS-CASPT2)Shiozaki et al. (2011) to remedy the erratic behavior of MS-CASPT2 potential energy surfaces when the mixing is strong.Serrano-Andrés et al. (2005) Note that this extension was first proposed by Granovsky for uncontracted multireference perturbation theory.Granovsky (2011) Very recently, One of the authors and co-workers developed an analytical nuclear gradient program for CASPT2,MacLeod and Shiozaki (2015); Vlaisavljevich and Shiozaki (2016) which forms the basis for the work presented herein.

In this work, we report the derivation and implementation of the analytical MS-CASPT2 and XMS-CASPT2 derivative couplings. The computer program has been implemented as an extension of the aforementioned nuclear gradient program. We note in passing that the interstate couplings studied for MS-CASPT2 (with partial internal contractionCelani and Werner (2003); Werner et al. (2011)) by Mori and KatoMori and Kato (2009) are part of the derivative coupling, which is what we call the MS-mixing term (see below). In the following, we first present the definition of the derivative couplings for MS-CASPT2 and XMS-CASPT2 wave functions, followed by the working equations for analytically evaluating the (X)MS-CASPT2 derivative couplings. We compare the minimum energy conical intersections (MECIs) of ethylene optimized by XMS-CASPT2 with those of unc-MRCI. In addition, the shapes of the potential energy surfaces near the MECIs of a cationic retinal model chromophore (PSB3) are investigated. We then present the MECIs of stilbene and an anionic GFP model chromophore (p-HBDI) optimized by XMS-CASPT2.

II Theoretical Background

In this section, we briefly review the previous work relevant to the development of the analytical (X)MS-CASPT2 derivative couplings. We first present the XMS-CASPT2 theory that is the basis of this work. The definition of the derivative coupling is then presented, followed by an algorithm for the analytical evaluation of the derivative couplings for SA-CASSCF.

II.1 XMS-CASPT2 Wave Functions

XMS-CASPT2 is a quasi-degenerate second-order perturbation theory on the basis of the CASSCF reference functions. The CASSCF wave functions are a linear combination of the Slater determinants,

[TABLE]

where $c_{I,L}$ are the configuration-interaction (CI) coefficients. In the following, $I$ and $J$ label Slater determinants, and $K$, $L$, $M$, and $N$ label reference functions. In XMS-CASPT2, the rotated reference functions are formed by diagonalizing the state-averaged Fock operator, $\hat{f}$, within the reference space,Granovsky (2011); Shiozaki et al. (2011); Vlaisavljevich and Shiozaki (2016)

[TABLE]

where $U_{LM}$ is chosen such that it satisfies

[TABLE]

$E^{(0)}_{M}$ is the zeroth-order energy of the state $|\tilde{M}\rangle$. For latter convenience, we introduce the rotated reference coefficients,

[TABLE]

The state-specific CASPT2 wave function is the sum of the reference function and the first-order corrections expanded in the internally contracted basis, i.e.,

[TABLE]

where $\Omega$ denotes all possible double-excitation manifold (see Refs. Andersson et al., 1990, 1992; MacLeod and Shiozaki, 2015; Vlaisavljevich and Shiozaki, 2016). For brevity, we introduce the following short-hand notations:

[TABLE]

The perturbation amplitudes $T_{\Omega,LN}$ are obtained by solving the amplitude equation,

[TABLE]

where $E_{\mathrm{shift}}$ is the level shift to circumvent intruder state problems.Roos and Andersson (1995) Once the amplitudes are determined, the effective Hamiltonian $H^{\text{eff}}$ is constructed as

[TABLE]

This effective Hamiltonian is then diagonalized to obtain the XMS-CASPT2 energies and wave functions,

[TABLE]

Hereafter $P$ and $Q$ label the physical states of interest. We define the XMS-rotated reference function,

[TABLE]

which will be used later.

II.2 Definition of Derivative Coupling in CASPT2

The Lagrangian for the total energy of a molecule, including the nuclear and electronic kinetic energy, resulting in the total wave function $\Psi_{\mathrm{tot}}$ is

[TABLE]

and solving the stationary condition for this Lagrangian, $\partial\mathcal{L}/\partial\Psi_{\mathrm{tot}}=0$, leads to the molecular Schrödinger equation. Here, the molecular Hamiltonian $\hat{\mathcal{H}}$ is written as a sum of the nuclear kinetic energy operator $\hat{T}_{\mathrm{nuc}}$ and the electronic Hamiltonian $\hat{H}_{\mathrm{el}}$,

[TABLE]

where $\mathbf{r}$ and $\mathbf{X}$ are the electronic and nuclear coordinates. The nuclear kinetic energy operator is explicitly written as

[TABLE]

where $A$ labels the nuclei, and $M_{A}$ is the nuclear mass. The expectation value for the total energy $\mathcal{E}$ with $\Psi_{\mathrm{tot}}$ is

[TABLE]

Within the Born–Oppenheimer approximation, one ignores the $\hat{T}_{\mathrm{nuc}}$ term by assuming that the kinetic energies of the nuclei are much smaller than the kinetic energies of the electrons. Then, the Lagrangian for the electronic Hamiltonian becomes

[TABLE]

where $|\Psi_{\mathrm{el}}\rangle$ is an electronic wave function. The stationary condition for this Lagrangian,

[TABLE]

leads to the electronic Schrödinger equation. The set of the solutions for this equation yields the adiabatic basis set (for example, the full CI eigenvectors and their corresponding energies). We will denote the individual adiabatic state as $|\Psi_{\mathrm{el},P}(\mathbf{r};\mathbf{X})\rangle$. To calculate the total energy with the motions of the electrons and nuclei using the electronic wave functions, we expand the total wave function in the adiabatic basis set

[TABLE]

where the expansion coefficients $\chi$ are normalized. Substituting this into the expression for the total energy $\mathcal{E}$ [Eq. (14)] yields the expression

[TABLE]

The off-diagonal elements of the nuclear kinetic energy term can be explicitly written as

[TABLE]

where $A$ labels the nuclei, and $M_{A}$ is the nuclear mass. The matrix element in the first term is the derivative coupling between the adiabatic states $Q$ and $P$, which we hereafter denote as $\mathbf{d}^{QP}$,

[TABLE]

The second term on the right hand side of Eq. (19) is the scalar coupling.Lengsfield et al. (1984); Tully (1990); Lengsfield and Yarkony (1992); Domcke et al. (2004); Baer (2006)

The derivative coupling formulas for (X)MS-CASPT2 can be obtained in a similar fashion. The Lagrangians for the zeroth-order energy and the energy corrected to second order (the Hylleraas functional) are

[TABLE]

where the zeroth-order Hamiltonian and the zeroth-order total energy are defined as

[TABLE]

$\hat{H}_{\mathrm{el}}^{(0)}$ is the zeroth-order Hamiltonian in XMS-CASPT2. The expressions for the energy corrected to second order $\mathcal{E}^{\mathrm{PT2}}$ with optimized $\Psi_{\mathrm{tot}}^{(0)}$ and $\Psi_{\mathrm{tot}}^{(1)}$ are

[TABLE]

Again, within the Born–Oppenheimer approximation, one ignores the $\hat{T}_{\mathrm{nuc}}$ term, and the Lagrangians become

[TABLE]

where $E_{\mathrm{el}}^{(0)}$ is the zeroth-order electronic energy [see Eq. (3)]. The stationary conditions for $\mathcal{L}^{(0)}_{\mathrm{BO},\mathbf{X}}$, and $\mathcal{L}^{(1)}_{\mathrm{BO},\mathbf{X}}$ with fixed $\Psi_{\mathrm{el}}^{(0)}(\mathbf{r};\mathbf{X})$,

[TABLE]

give us a set of the zeroth- and first-order wave functions [the latter corresponds to Eq. (7)]. Analogous to Eq. (17), we expand the total wave functions in the basis of the electronic wave functions as

[TABLE]

By substituting these into the total second-order energy $\mathcal{E}^{\mathrm{PT2}}$ [Eq. (23)], we obtain

[TABLE]

Note that this substitution is an approximation if the total wave functions are not fully optimized. In Eq. (27), the electronic wave functions are rotated to make $\hat{H}_{\mathrm{el}}$ diagonal [which is Eq. (9a)],

[TABLE]

The off-diagonal elements of the $\hat{T}_{\mathrm{nuc}}$ term can therefore be explicitly written as

[TABLE]

which leads to the following expression for the derivative coupling in (X)MS-CASPT2,

[TABLE]

The validity of the expression for $\mathbf{d}^{\mathrm{XMS},QP}$ has been confirmed by calculating the line integral around conical intersections (see below).Domcke et al. (2004); Baer (2006) This quantity can be readily used in many formulations of non-adiabatic dynamics, for instance in fewest-switch surface-hopping (FSSH) non-adiabatic dynamics.Tully (1990); Fabiano et al. (2008) Using the derivative coupling defined in Eq. (30), the equation for XMS-CASPT2 FSSH dynamics that propagates $\chi_{Q}$ is

[TABLE]

where we introduce the velocity of the trajectory as $\mathbf{v}=d\mathbf{X}/dt$.

II.3 Analytical SA-CASSCF Derivative Coupling

The derivative coupling between SA-CASSCF wave functions can be written as

[TABLE]

We call $\mathbf{d}^{QP}_{\mathrm{CI}}$ and $\mathbf{d}^{QP}_{\mathrm{det}}$ the CI and determinant terms, respectively.Lischka et al. (2004); Galvan et al. (2016); Lengsfield et al. (1984); Domcke et al. (2004) The CI term can be evaluated using the following relationship,

[TABLE]

where $E^{\mathrm{CAS}}_{P}$ and $E^{\mathrm{CAS}}_{Q}$ are the CASSCF energies of states $P$ and $Q$. Note that this term corresponds to the so-called interstate couplings,Domcke et al. (2004); Bearpark et al. (1994); Manaa and Yarkony (1993); Mori and Kato (2009); Matsika and Krause (2011) and are analogous to the Hellmann–Feynman forces in nuclear gradient theory.Domcke et al. (2004); Lengsfield and Yarkony (1992); Matsika and Krause (2011) This formula follows from the fact that the off-diagonal elements of the Hamiltonian in the basis of the SA-CASSCF states are zero, $\matrixelement{Q}{\hat{H}}{P}=0$. Since the SA-CASSCF wave functions are optimized with respect to both the CI and orbital coefficients, we can write the LagrangianCelani and Werner (2003); Koch et al. (1990) for the CI term multiplied by $E_{P}^{\mathrm{CAS}}-E_{Q}^{\mathrm{CAS}}$ as

[TABLE]

The first term is the off-diagonal element of the CI Hamiltonian, and the remaining terms represent the SA-CASSCF convergence conditions, where $\mathbf{Z}$, $\mathbf{V}$, $\mathbf{z}_{N}$, and $x_{N}$ are their Lagrange multipliers. $\mathbf{A}$ is the derivative of the SA-CASSCF energy with respect to the orbital rotations, $\mathbf{C}$ and $\mathbf{S}$ are the orbital coefficients and overlap matrix elements, and $W_{N}$ and $\mathbf{c}_{N}$ are the weight and CI coefficients of state $N$. See details in Refs. Shiozaki et al., 2011 and Celani and Werner, 2003. The so-called $Z$-vector equationHandy and Schaefer (1984); Lengsfield and Yarkony (1992) is obtained by differentiating the Lagrangian with respect to the orbital rotation parameters and CI coefficients; the source terms of the $Z$-vector equation are

[TABLE]

The multiplier-dependent terms in the $Z$-vector equation are identical to those in nuclear gradient theory.Celani and Werner (2003); Shiozaki et al. (2011); Győrffy et al. (2013) After solving the $Z$-vector equation for $\mathbf{Z}$ and $\mathbf{z}$, we compute the effective density matrices ($\gamma^{\mathrm{eff}}$ and $\Gamma^{\mathrm{eff}}$) and the Lagrange multipliers $\mathbf{V}$, whose explicit forms are given in Refs. Shiozaki et al., 2011 and Celani and Werner, 2003. These matrices are then contracted with the derivative integrals to yield the CI term:

[TABLE]

where $\mu$, $\nu$, $\lambda$ and $\sigma$ label atomic orbitals, and the superscript ${\mathbf{X}}$ denotes the integral derivatives with respect to $\mathbf{X}$. The formula for the determinant term can be easily derived if one rewrites the operator $d/d\mathbf{X}$ as a one-electron operator,Lengsfield et al. (1984)

[TABLE]

in which $\phi_{\mu}$ are atomic orbitals, $d\kappa_{rs}/d\mathbf{X}$ are the orbital response parameters, and $\hat{E}_{rs}$ is the spin-adapted one-particle excitation operator. This leads to a compact form of the determinant term,

[TABLE]

where $\gamma^{QP}_{rs}=\langle Q|\hat{E}_{rs}|P\rangle$. In practice, the evaluation of the derivative of $\kappa_{rs}$ can be avoided in one of two ways: One approach is to include it in the $Z$-vector algorithm by adding the following $Y^{C}_{rs}$ to $Y_{rs}$,

[TABLE]

The other approach is to use the following equivalent expressionLischka et al. (2004); Galvan et al. (2016) that is written in terms of $\sigma_{rs}^{X}$ alone,

[TABLE]

The latter approach does not require modification of the $Z$-vector equation. A similar algorithm to the former is used in the evaluation of the XMS-CASPT2 derivative couplings described below.

III Analytical XMS-CASPT2 Derivative Coupling

The XMS-CASPT2 derivative coupling can be formally written as

[TABLE]

where we introduced a short-hand notation for the products of the MS mixing matrices,

[TABLE]

The expressions in Eq. (42) are manifestly invariant with respect to any rotations among active orbitals. The orbital invariance has also been confirmed numerically. In the following, we present the working equations for its analytical evaluation.

Since the XMS-CASPT2 states are obtained from the eigenvalue equation, Eq. (9), the MS-mixing term $\mathbf{d}^{\mathrm{XMS},QP}_{\mathrm{mix}}$ can be evaluated using a strategy similar to the one for the CI term [Eq. (33)] in the evaluation of the SA-CASSCF derivative coupling,

[TABLE]

in which we define $\Delta E_{PQ}^{\mathrm{XMS}}=E_{P}^{\text{XMS}}-E_{Q}^{\text{XMS}}$ for brevity. This MS-mixing term corresponds to the interstate coupling and can be evaluated using the methodologies described in previous workMori and Kato (2009); Vlaisavljevich and Shiozaki (2016) using the PT2 Lagrangian,

[TABLE]

The so-called $\lambda$-equation is obtained by making the PT2 Lagrangian stationary with respect to the perturbation amplitudes, $T_{\Omega,KN}$:

[TABLE]

The explicit expression for the total Lagrangian is

[TABLE]

The second, third, and fourth terms on the right-hand side account for the SA-CASSCF stationary conditions and also appeared in Eq. (34). The last two terms account for the frozen-core approximation and the XMS rotation, respectively. Using this Lagrangian, the MS-mixing term is simply written as

[TABLE]

where the subscript $R$ indicates that the MS rotation matrix elements in the first term of the Lagrangian are held fixed when taking the derivative. The CAS and PT2 terms are

[TABLE]

To avoid the evaluation of the derivatives of $U_{MK}$, $c_{I,N}$, and $\kappa_{rs}$ with respect to $\mathbf{X}$, we simultaneously evaluate Eqs. (48) and (49) using the $Z$-vector equation. First, the Lagrange multiplier $w_{MN}$ is evaluated as

[TABLE]

The $Z$-vector algorithm to calculate $\mathbf{Z}$, $\mathbf{z}_{N}$, $\mathbf{V}$, and $x_{N}$ is analogous to the one used in the nuclear gradient algorithm;Celani and Werner (2003); Shiozaki et al. (2011); Vlaisavljevich and Shiozaki (2016) the source terms for the $Z$-vector equation are

[TABLE]

The contribution from the second term in the square bracket in Eq. (51b) vanishes in the case of XMS-CASPT2 (but not in MS-CASPT2), because $y_{I,N}$ is taken to be orthogonal to the reference space when XMS reference functions are used. The $Z$-vector equation is solved iteratively. Having determined $w_{MN}$, $\mathbf{Z}$, $\mathbf{z}_{N}$, and $\mathbf{V}$, the derivative couplings $\mathbf{d}^{\mathrm{XMS},QP}$ can be computed as

[TABLE]

The first term is computed as a contraction of the effective density matrices with the derivative integrals as in the nuclear gradient algorithms. The algorithm for MS-CASPT2 can be obtained by setting $U_{MN}=\delta_{MN}$ and neglecting its derivative. If desired, translational invariance can be achieved by setting the second term to zero.Fatehi et al. (2011); Fatehi and Subotnik (2012) Note that density fitting is used when evaluating the above expressions in our implementation. The additional computational cost associated with the CAS and PT2 terms is negligible compared to the costs of computing the MS-mixing term.

IV Numerical Examples

In this section, we first present the MECIs of ethylene computed by XMS-CASPT2 and compare their structures and energetics with those computed by unc-MRCI. Second, we investigate the shapes of the potential energy surfaces of a model retinal chromophore (the penta-2,4-dieniminium cation or PSB3) near the MECIs. Finally, we report the computation of the MECIs for the large molecules [stilbene and an anionic GFP model chromophore (4-para-hydroxybenzylidene-1,2-dimethyl-imidazolin-5-one or pHBDI)]. The geometries were optimized using XMS-CASPT2 with the cc-pVDZ basis set and its corresponding JKFIT basis set for density fitting unless mentioned otherwise. The so-called SS-SR contraction schemeVlaisavljevich and Shiozaki (2016) with vertical shift (0.2 $E_{\mathrm{h}}$) was used unless otherwise specified. We searched for the MECI by the gradient projection method,Bearpark et al. (1994) that uses the gradient difference and interstate coupling vectors (instead of the full derivative coupling vector). The flowchart method was used for updating model Hessians for quasi-Newton steps.Birkholz and Schlegel (2016) For comparison, the MECIs were also optimized using SA-CASSCF. All of the calculations were performed using the program package bagel.bag

IV.1 Conical Intersections of Ethylene

First, to assess the accuracy of XMS-CASPT2 against unc-MRCI,Barbatti et al. (2004); Zhang and Herbert (2014) we optimized the MECIs of ethylene. We adopted the reference SA-CASSCF wave functions from Ref. Barbatti et al., 2004 with the (2e,2o) active space using the aug-cc-pVDZ basis set. The three lowest states were averaged in SA-CASSCF. The XMS-CASPT2 calculations were performed using a vertical shift of 0.5 $E_{\mathrm{h}}$. All of the MECIs reported using unc-MRCIBarbatti et al. (2004); Zhang and Herbert (2014) were also found by XMS-CASPT2: the twisted-pyramid (Py), ethylidene (Et), $C_{3v}$ ethylidene ($C_{3v}$-Et), H-migration (Hm) ($\mathrm{S}_{0}$/$\mathrm{S}_{1}$), and twisted-orthogonal (To) ($\mathrm{S}_{1}$/$\mathrm{S}_{2}$) MECIs.

The geometries of the MECIs and the energies relative to the $\mathrm{S}_{0}$ minimum are shown and compared with the unc-MRCI results in Fig. 1 and in Table 1. They are in good agreement with those computed by unc-MRCI with the Davidson(+Q) correction and the additional restricted active space in the reference function. This suggests that XMS-CASPT2 can yield results comparable to unc-MRCI+Q in finding the geometries and energies of MECIs, as in the case for calculating the spectroscopic constants of a range of small molecules.Azizi et al. (2006)

IV.2 Conical Intersection of PSB3

Next, we show the potential energy surfaces of cationic PSB3 near the MECIs. PSB3 is a minimal model for the photoisomerization of protonated Schiff bases,Gozem et al. (2012, 2014) whose conical intersections have been characterized by various electronic structure methods including MS-CASPT2 (using numerical nuclear gradients).Gozem et al. (2012) In this work, we optimized the geometries of the MECIs using SA-CASSCF, MS-CASPT2, and XMS-CASPT2. The two lowest states were averaged in SA-CASSCF. We used the (6e,6o) active space consisting of three $\pi$ and three $\pi^{\ast}$ orbitals. The vertical shift of 0.5 $E_{\mathrm{h}}$ was applied. The energies at the MECI geometries relative to the trans structure in the ground state were found to be 2.36 eV with both MS-CASPT2 and XMS-CASPT2, which is lower by 0.4 eV than that computed by CASSCF (2.74 eV). This attests to the importance of dynamical correlation in predicting energies at MECIs.

On the other hand, the shapes of the potential energy surfaces around the MECIs computed by MS-CASPT2 and XMS-CASPT2 were found to differ significantly. The difference is schematically shown in Fig. 2(a) where the energy gaps are presented along the loop on the branching plane [Fig. 2(b)] centered at the MECI, as in Ref. Gozem et al., 2014. Note that the branching plane is the plane defined by the gradient difference and interstate coupling vectors ($\mathbf{g}$ and $\mathbf{h}$, respectively). The radius of the loop was set to 0.002 Å. We normalized the vectors $\mathbf{g}$ and $\mathbf{h}$ when generating the loop. The MS-CASPT2 energy gaps showed spurious oscillatory behavior along the loop, as reported in Ref. Gozem et al., 2014, whose amplitude was as large as 0.02 eV. The gaps computed by XMS-CASPT2 were completely smooth. The amplitude of the oscillation (0.84 kcal/mol) was comparable to that computed by unc-MRCI ($\sim$1 kcal/mol).Gozem et al. (2014) The line integral of the XMS-CASPT2 derivative coupling along the same loopDomcke et al. (2004); Baer (2006) was found to be very close to $\pi$ (3.1412). The potential energy surfaces near the MECIs are also presented in Fig. 2(c). This example stresses the importance of using XMS-CASPT2 (instead of MS-CASPT2) when simulating the electronic structure around conical intersections.

IV.3 Conical Intersections of Stilbene

Stilbene has been extensively studied computationally as a model compound for photoisomerization around a C=C double bond,Quenneville and Martínez (2003); Levine and Martínez (2007); Minezawa and Gordon (2011); Ioffe and Granovsky (2013); Lei et al. (2014); Harabuchi et al. (2014) partly because reliable experimental results are available.Waldeck (1991) The mechanisms for the photodynamics and the locations of the conical intersections have been well characterized. There are three types of low-lying conical intersections: one-bond flip (OBF), hula-twist (HT), and 4a,4b-dihydrophenanthrene(DHP)-like conical intersections. The OBF and HT conical intersections are related to the cis–trans photoisomerization, and the DHP-like conical intersections are relevant in the formation of DHP, which is a minor product. We adopted the CASSCF MECI geometries reported in Refs. Ioffe and Granovsky, 2013 and Lei et al., 2014 as the starting geometries for the optimizations. The reference CASSCF wave functions were calculated using two-state averaging with the (6e,6o) active space.

The optimized structures at the $\mathrm{S}_{0}$ minima and MECIs are shown in Figure 3 along with their energies. With XMS-CASPT2, two of the MECIs are located below the trans Franck–Condon point, and the other two are above that point. This is in contrast with SA-CASSCF results in which all of the MECIs are about 1 eV below the trans Franck–Condon point. That is, the MECIs are less accessible from the Franck–Condon point when dynamical correlation is considered. This result is consistent with the previous study based on CASPT2//CASSCF/6-31G.Lei et al. (2014)

Next, we computed the relative energies of the MECIs using SA-CASSCF and XMS-CASPT2. The OBF1 MECIs are 0.21 eV and 0.02 eV lower than the HT1 MECIs with SA-CASSCF using (6e,6o) and (10e,10o) active spaces, which is different from the XMS-CASPT2 results (0.24 eV higher). This is consistent with previous SF-TDDFTMinezawa and Gordon (2011) and XMCQDPT2Ioffe and Granovsky (2013) results. The DHP-like MECIs are also stabilized by dynamical correlation. These MECIs have higher energies than the other MECIs when dynamical correlation is not included in the calculations; they are 0.50 eV and 1.18 eV higher than the OBF1 MECIs with SA-CASSCF using (6e,6o) and (10e,10o) active spaces, respectively. Note that the DHP-like conical intersections become unfavorable when the active space is extended,Ioffe and Granovsky (2013) and therefore, further investigations using XMS-CASPT2 and larger active spaces are required to unravel the importance of the DHP-like conical intersections in photodynamics. The structures of the HT1 MECI of stilbene located using SA-CASSCF and XMS-CASPT2 are shown in Figure 4 along with the branching planes. The root mean square deviation of the MECI structures using SA-CASSCF and XMS-CASPT2 was 0.36 Å after aligning the two to minimize the deviation, which is the accuracy that one would achieve when the CASPT2//CASSCF approach is used.

IV.4 Conical Intersections of pHBDI

Next, we show the MECIs of the 4-para-hydroxybenzylidene-1,2-dimethyl-imidazolin-5-one (pHBDI) anion. Anionic pHBDI is considered to be an emitting species of the green fluorescent protein (GFP) and its variants.Tsien (1998); Weber et al. (1999); Martin et al. (2004); Olsen and Smith (2008); Park and Rhee (2016); van Thor (2009); Olsen et al. (2010); Polyakov et al. (2010); Acharya et al. (2017); Morozov and Groenhof (2015); Minezawa and Gordon (2012); Zhang et al. (2014) While GFP exhibits strong fluorescence with a lifetime on the order of nanoseconds,Tsien (1998); Martin et al. (2004); Olsen and Smith (2008); van Thor (2009); Olsen et al. (2010); Polyakov et al. (2010); Park and Rhee (2016); Acharya et al. (2017) the nonadiabatic transition is known to occur in about a few picoseconds when the chromophore is not embedded in the protein environment.Mandal et al. (2004) As a resonant monomethine dye, it is widely accepted that the anionic GFP chromophore undergoes nonadiabatic transitions when the chromophore is twisted along the bridge.Olsen and Smith (2008); Olsen et al. (2010); Olsen (2010); Morozov and Groenhof (2015) There are two available bridge channels in this molecule, which are the imidazolinone (I) and phenolate (P) channels, named after the moiety connected to the bridge bond that twists.Olsen et al. (2010)

We optimized the planar equilibrium geometry of the ground state, the geometries of the I- and P-twisted minima of the first excited state, and the MECIs between these states near the twisted geometries. The reference CASSCF wave functions were optimized using three-state averaging with the (4e,3o) active space, which provides a reliable description of twisting of the bridge bond.Olsen and Smith (2008); Olsen and McKenzie (2009)

The structures and the energy diagram for these geometries are shown in Figure 5. The twisted geometries optimized on the first excited state were found to be lower in energy than the Franck–Condon point by 0.38 and 0.14 eV for the I- and P-twisted minima, respectively. However, the conical intersections associated with these twisted geometries are located much higher in energy than the minima; the I-twisted MECI has about the same energy as the Franck–Condon point, and the P-twisted MECI lies 0.44 eV above the Franck–Condon point. This result is in stark contrast to those obtained using SA-CASSCF: the SA-CASSCF energy at the I-twisted and P-twisted MECIs are 1.24 and 0.64 eV lower than $\mathrm{S}_{1}$ energy at the Franck–Condon point. This suggests that the thermal accessibility of the $\mathrm{S}_{0}$/$\mathrm{S}_{1}$ P-twisted CI predicted by SA-CASSCF is an artifact, because the P-twisted MECI lies ca. 10 kcal/mol above the Franck–Condon point on the XMS-CASPT2 surface. SA-CASSCF results remain qualitatively the same even when the active space is expanded to (12e,11o); the MECIs computed using CAS(12e,11o) were 0.91 and 0.30 eV lower than the Franck–Condon point (see also Ref. Polyakov et al., 2010). It is notable that this differs qualitatively from the previous CAS(12e,11o) result with two-state averaging,Martin et al. (2004) which suggests that the choice of the states to be averaged is as important as the choice of the active space.

The structures of the P-twisted MECI optimized using SA-CASSCF and XMS-CASPT2 are shown in Figure 6, in which the full derivative coupling vectors are also shown. The MECI geometries differ primarily on the bridge: The hydrogen atom on the bridge comes out of the plane of the imidazolinone ring more severely in XMS-CASPT2 than in SA-CASSCF. The torsional angle around the imidazolinone bridge (H1-C2-C3-C4 torsional angle, see Figure 6) was $-45.1^{\circ}$ and $-32.0^{\circ}$ for XMS-CASPT2 and SA-CASSCF, respectively. Furthermore, the bond length for the imidazolinone bridge (C2-C3 bond) was found to be 1.49 Å (XMS-CASPT2) and 1.45 Å (SA-CASSCF), indicating that the bridging carbon atom has slightly more $sp^{3}$ character than $sp^{2}$ at the P-twisted MECI when dynamical correlation is taken into account. We also note that the hula twist (HT) conical intersection reported in Refs. Weber et al., 1999 and Martin et al., 2004 was not found in our optimization (the optimization converged to the I-twisted conical intersection). Overall, our results imply that including dynamical correlation quantitatively (or, even qualitatively) affects the photochemical dynamics.

V Conclusions

In this work, we have derived the working equations for analytically evaluating the derivative couplings using (X)MS-CASPT2. The equations have been translated into an efficient computer program as an extension of previously reported analytical gradient programs for (X)MS-CASPT2, which is included in the bagel packagebag and publicly available under the GNU General Public License. The fully internally contracted wave functions were used together with the density fitting approximation. The computational cost for calculating the derivative couplings was found to be essentially the same as that for computing the nuclear energy gradients for one state. The energetics at the MECIs were substantially influenced by dynamical correlation. The optimization of MECIs for ethylene, PSB3, stilbene, and pHBDI was presented to demonstrate the versatility of our program. This finding encourages us to develop a methodology for large-scale direct dynamics in complex systems using a model that takes into account dynamical correlation contributions such as XMS-CASPT2. To achieve this goal, efforts to improve our algorithms and implementations (especially for large active spaces) are underway and will be reported in the near future.

VI Acknowledgments

The authors thank Dr. Bess Vlaisavljevich for her constructive comments on the manuscript. The debugging of the SA-CASSCF derivative coupling code in bagel, on which this work is based, was facilitated by the existing implementation in molpro.Werner et al. (2011) This work has been supported by the Air Force Office of Scientific Research Young Investigator Program (Grant No. FA9550-15-1-0031). The development of the program infrastructure has been in part supported by National Science Foundation [ACI-1550481 (JWP) and CHE-1351598 (TS)]. T.S. is an Alfred P. Sloan Research Fellow.

Appendix A Numerical XMS-CASPT2 Derivative Coupling

In part to validate our new analytical derivative coupling program, we have also implemented the code that numerically evaluates (X)MS-CASPT2 derivative coupling. In this case, $\mathbf{d}^{\mathrm{XMS},QP}_{\mathrm{mix}}$ and $\mathbf{d}^{\mathrm{XMS},QP}_{\mathrm{CAS}}$ can be combined to the form that resembles SA-CASSCF derivative coupling, using ${\tilde{r}}_{I,Q}=\sum_{M}\tilde{c}_{I,M}R_{MQ}$, as

[TABLE]

and is evaluated numerically asGalloy and Lorquet (1977)

[TABLE]

The PT2 term is calculated using the following expression:

[TABLE]

The derivatives ($\Delta{\tilde{r}}_{I,P}/\Delta X$ and $\Delta C_{\mu r}/\Delta X$) were calculated by means of the finite difference formula. The numerical derivative couplings agreed with those analytically evaluated.