Analytical derivative coupling for multistate CASPT2 theory

Analytical derivative coupling for multistate CASPT2 theory

Jae Woo Park    Toru Shiozaki Department of Chemistry, Northwestern University, 2145 Sheridan Rd., Evanston, IL 60208, USA.
August 22, 2019

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).

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,


where are the configuration-interaction (CI) coefficients. In the following, and label Slater determinants, and , , , and label reference functions. In XMS-CASPT2, the rotated reference functions are formed by diagonalizing the state-averaged Fock operator, , within the reference space,Granovsky (2011); Shiozaki et al. (2011); Vlaisavljevich and Shiozaki (2016)


where is chosen such that it satisfies


is the zeroth-order energy of the state . For latter convenience, we introduce the rotated reference coefficients,


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.,


where 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:


The perturbation amplitudes are obtained by solving the amplitude equation,


where is the level shift to circumvent intruder state problems.Roos and Andersson (1995) Once the amplitudes are determined, the effective Hamiltonian is constructed as


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


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


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 is


and solving the stationary condition for this Lagrangian, , leads to the molecular Schrödinger equation. Here, the molecular Hamiltonian is written as a sum of the nuclear kinetic energy operator and the electronic Hamiltonian ,


where and are the electronic and nuclear coordinates. The nuclear kinetic energy operator is explicitly written as


where labels the nuclei, and is the nuclear mass. The expectation value for the total energy with is


Within the Born–Oppenheimer approximation, one ignores the 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


where is an electronic wave function. The stationary condition for this Lagrangian,


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 . 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


where the expansion coefficients are normalized. Substituting this into the expression for the total energy [Eq. (14)] yields the expression


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


where labels the nuclei, and is the nuclear mass. The matrix element in the first term is the derivative coupling between the adiabatic states and , which we hereafter denote as ,


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


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


is the zeroth-order Hamiltonian in XMS-CASPT2. The expressions for the energy corrected to second order with optimized and are


Again, within the Born–Oppenheimer approximation, one ignores the term, and the Lagrangians become


where is the zeroth-order electronic energy [see Eq. (3)]. The stationary conditions for , and with fixed ,


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


By substituting these into the total second-order energy [Eq. (23)], we obtain


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 diagonal [which is Eq. (9a)],


The off-diagonal elements of the term can therefore be explicitly written as


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


The validity of the expression for 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 is


where we introduce the velocity of the trajectory as .

ii.3 Analytical SA-CASSCF Derivative Coupling

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


We call and 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,


where and are the CASSCF energies of states and . 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, . 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 as


The first term is the off-diagonal element of the CI Hamiltonian, and the remaining terms represent the SA-CASSCF convergence conditions, where , , , and are their Lagrange multipliers. is the derivative of the SA-CASSCF energy with respect to the orbital rotations, and are the orbital coefficients and overlap matrix elements, and and are the weight and CI coefficients of state . See details in Refs. Shiozaki et al., 2011 and Celani and Werner, 2003. The so-called -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 -vector equation are


The multiplier-dependent terms in the -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 -vector equation for and , we compute the effective density matrices ( and ) and the Lagrange multipliers , 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:


where , , and label atomic orbitals, and the superscript denotes the integral derivatives with respect to . The formula for the determinant term can be easily derived if one rewrites the operator as a one-electron operator,Lengsfield et al. (1984)


in which are atomic orbitals, are the orbital response parameters, and is the spin-adapted one-particle excitation operator. This leads to a compact form of the determinant term,


where . In practice, the evaluation of the derivative of can be avoided in one of two ways: One approach is to include it in the -vector algorithm by adding the following to ,


The other approach is to use the following equivalent expressionLischka et al. (2004); Galvan et al. (2016) that is written in terms of alone,


The latter approach does not require modification of the -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


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


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 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,


in which we define 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,


The so-called -equation is obtained by making the PT2 Lagrangian stationary with respect to the perturbation amplitudes, :


The explicit expression for the total Lagrangian is


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


where the subscript 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


To avoid the evaluation of the derivatives of , , and with respect to , we simultaneously evaluate Eqs. (48) and (49) using the -vector equation. First, the Lagrange multiplier is evaluated as


The -vector algorithm to calculate , , , and 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 -vector equation are


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 is taken to be orthogonal to the reference space when XMS reference functions are used. The -vector equation is solved iteratively. Having determined , , , and , the derivative couplings can be computed as


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 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 ) 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 . 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), ethylidene (-Et), H-migration (Hm) (/), and twisted-orthogonal (To) (/) MECIs.

Figure 1: The geometries of the MECIs of ethylene optimized by XMS-CASPT2. The selected geometrical parameters are shown with the unc-MRCI geometrical parametersBarbatti et al. (2004) in the parentheses (bondlengths are shown in Å).
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
-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

a Taken from Ref. 56.

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

The geometries of the MECIs and the energies relative to the 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 and three orbitals. The vertical shift of 0.5  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.

Figure 2: The potential energy surfaces near the MECIs of PSB3. (a) The and energy gaps around the 0.002 Å radius loop centered at the MECIs, with MS-CASPT2 (dashed) and XMS-CASPT2 (solid). The loop is on the branching plane. (b) The branching plane at the MECI computed by MS-CASPT2 (left) and XMS-CASPT2 (right). The vectors and are plotted with blue and red arrows, respectively. (c) The and potential energies on the branching plane, with MS-CASPT2 (left) and XMS-CASPT2 (right). The energy at the MECI is set to zero.

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 ( and , respectively). The radius of the loop was set to 0.002 Å. We normalized the vectors and 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 (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 (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 cistrans 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.

Figure 3: XMS-CASPT2 optimized , geometries and MECIs of stilbene. Energies (in eV) are reported relative to the minimum energy.

The optimized structures at the 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.

Figure 4: (a) Overlay of the HT1 MECI structures of stilbene optimized by SA-CASSCF (pink) and XMS-CASPT2 (red); (b) The branching plane at the MECI computed by SA-CASSCF (left) and XMS-CASPT2 (right). The vectors and are plotted with blue and red arrows, respectively.

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)

Figure 5: XMS-CASPT2 optimized , geometries and conical intersections of pHBDI. Energies (in eV) are reported relative to the minimum energy.

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 energy at the Franck–Condon point. This suggests that the thermal accessibility of the / 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.

Figure 6: (a) Overlay of the P-twisted MECI structures of pHBDI optimized by SA-CASSCF (pink) and XMS-CASPT2 (red); (b) The branching plane at the MECI computed by SA-CASSCF (left) and XMS-CASPT2 (right). The vectors and are plotted with blue and red arrows, respectively.

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 and 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 character than 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, and can be combined to the form that resembles SA-CASSCF derivative coupling, using , as


and is evaluated numerically asGalloy and Lorquet (1977)


The PT2 term is calculated using the following expression:


The derivatives ( and ) were calculated by means of the finite difference formula. The numerical derivative couplings agreed with those analytically evaluated.


  • Valeur (2002) B. Valeur, Molecular Fluorescence: Principles and Applications (Wiley-VCH, Weinheim, Germany, 2002) pp. 20–124.
  • Levine and Martínez (2007) B. G. Levine and T. J. Martínez, Annu. Rev. Phys. Chem. 58, 613 (2007).
  • Erbas-Cakmak et al. (2015) S. Erbas-Cakmak, D. A. Leigh, C. T. McTernan,  and A. L. Nussbaumer, Chem. Rev. 115, 10081 (2015).
  • Matsika and Krause (2011) S. Matsika and P. Krause, Annu. Rev. Phys. Chem 62, 621 (2011).
  • Matsika et al. (2014) S. Matsika, X. Feng, A. V. Luzanov,  and A. I. Krylov, J. Phys. Chem. A 118, 11943 (2014).
  • Olsen et al. (2010) S. Olsen, K. Lamothe,  and T. J. Martínez, J. Am. Chem. Soc. 132, 1192 (2010).
  • Lengsfield et al. (1984) B. H. Lengsfield, III, P. Saxe,  and D. R. Yarkony, J. Chem. Phys. 81, 4549 (1984).
  • Lengsfield and Yarkony (1992) B. H. Lengsfield, III and D. R. Yarkony, Adv. Chem. Phys. 82, 1 (1992).
  • Galvan et al. (2016) I. F. Galvan, M. G. Delcey, T. B. Pedersen, F. Aquilante,  and R. Lindh, J. Chem. Theory Comput. 12, 3636 (2016).
  • Lischka et al. (2004) H. Lischka, M. Dallos, P. G. Szalay, D. R. Yarkony,  and R. Shepard, J. Chem. Phys. 120, 7322 (2004).
  • Barbatti and Lischka (2008) M. Barbatti and H. Lischka, J. Am. Chem. Soc. 130, 6831 (2008).
  • Tajti and Szalay (2009) A. Tajti and P. G. Szalay, J. Chem. Phys. 131, 124104 (2009).
  • Christiansen (1999) O. Christiansen, J. Chem. Phys. 110, 711 (1999).
  • Ichino et al. (2009) T. Ichino, J. Gauss,  and J. F. Stanton, J. Chem. Phys. 130, 174105 (2009).
  • Fatehi et al. (2011) S. Fatehi, E. Alguire, Y. Shao,  and J. E. Subotnik, J. Chem. Phys. 135, 234105 (2011).
  • Chernyak and Mukamel (2000) V. Chernyak and S. Mukamel, J. Chem. Phys. 112, 3572 (2000).
  • Send and Furche (2010) R. Send and F. Furche, J. Chem. Phys. 132, 044107 (2010).
  • Zhang and Herbert (2014) X. Zhang and J. M. Herbert, J. Chem. Phys. 141, 064104 (2014).
  • Zhang and Herbert (2015) X. Zhang and J. M. Herbert, J. Chem. Phys. 143, 234107 (2015).
  • Herbert et al. (2016) J. M. Herbert, X. Zhang, A. F. Morrison,  and J. Liu, Acc. Chem. Res. 49, 931 (2016).
  • (21) The terms ’derivative coupling’ and ’nonadiabatic coupling’ are conventionally used interchangeably, except for some recent reports, such as Ref. 9.
  • Tully (1990) J. C. Tully, J. Chem. Phys. 93, 1061 (1990).
  • Nelson et al. (2014) T. Nelson, S. Fernandez-Alberti, A. E. Roitberg,  and S. Tretiak, Acc. Chem. Res. 47, 1155 (2014).
  • Tavernelli (2015) I. Tavernelli, Acc. Chem. Res. 48, 792 (2015).
  • Domcke et al. (2004) W. Domcke, D. R. Yarkony,  and H. Köppel, eds., Conical Intersections: Theory, Computation and Experiment (World Scientific, Singapore, 2004) pp. 3–174.
  • Baer (2006) M. Baer, Beyond Born–Oppenheimer: Electronic Nonadiabatic Coupling Terms and Conical Intersections (John Wiley & Sons, Hoboken, NJ, 2006) pp. 26–138.
  • Bearpark et al. (1994) M. J. Bearpark, M. A. Robb,  and H. B. Schlegel, Chem. Phys. Lett. 223, 269 (1994).
  • Manaa and Yarkony (1993) M. R. Manaa and D. R. Yarkony, J. Chem. Phys. 99, 5251 (1993).
  • Yarkony (1996) D. R. Yarkony, Rev. Mod. Phys. 68, 985 (1996).
  • Levine et al. (2006) B. G. Levine, C. Ko, J. Quenneville,  and T. J. Martínez, Mol. Phys. 104, 1039 (2006).
  • Huix-Rotllant et al. (2013) M. Huix-Rotllant, M. Filatov, S. Gozem, I. Schapiro, M. Olivucci,  and N. Ferré, J. Chem. Theory Comput. 9, 3917 (2013).
  • Li et al. (2014) S. L. Li, A. V. Marenich, X. Xu,  and D. G. Truhlar, J. Phys. Chem. Lett. 5, 322 (2014).
  • Tuna et al. (2015) D. Tuna, D. Lefrancois, Ł. Wolański, S. Gozem, I. Schapiro, T. Andruniów, A. Dreuw,  and M. Olivucci, J. Chem. Theory Comput. 11, 5758 (2015).
  • Epifanovsky and Krylov (2007) E. Epifanovsky and A. I. Krylov, Mol. Phys. 105, 2515 (2007).
  • Shepler et al. (2008) B. C. Shepler, E. Epifanovsky, P. Zhang, J. M. Bowman, A. I. Krylov,  and K. Morokuma, J. Phys. Chem. A 112, 13267 (2008).
  • Minezawa and Gordon (2009) N. Minezawa and M. S. Gordon, J. Phys. Chem. A 113, 12749 (2009).
  • Andersson et al. (1990) K. Andersson, P.-Å. Malmqvist, B. O. Roos, A. J. Sadlej,  and K. Wolinski, J. Phys. Chem. 94, 5483 (1990).
  • Andersson et al. (1992) K. Andersson, P.-Å. Malmqvist,  and B. O. Roos, J. Chem. Phys. 96, 1218 (1992).
  • Finley et al. (1998) J. Finley, P.-Å. Malmqvist, B. O. Roos,  and L. Serrano-Andrés, Chem. Phys. Lett. 288, 299 (1998).
  • Shiozaki et al. (2011) T. Shiozaki, W. Győrffy, P. Celani,  and H.-J. Werner, J. Chem. Phys. 135, 081106 (2011).
  • Serrano-Andrés et al. (2005) L. Serrano-Andrés, M. Merchán,  and R. Lindh, J. Chem. Phys. 122, 104107 (2005).
  • Granovsky (2011) A. A. Granovsky, J. Chem. Phys. 134, 214113 (2011).
  • MacLeod and Shiozaki (2015) M. K. MacLeod and T. Shiozaki, J. Chem. Phys. 142, 051103 (2015).
  • Vlaisavljevich and Shiozaki (2016) B. Vlaisavljevich and T. Shiozaki, J. Chem. Theory Comput. 12, 3781 (2016).
  • Celani and Werner (2003) P. Celani and H.-J. Werner, J. Chem. Phys. 119, 5044 (2003).
  • Werner et al. (2011) H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby,  and M. Schütz, WIREs Comput. Mol. Sci. 2, 242 (2011).
  • Mori and Kato (2009) T. Mori and S. Kato, Chem. Phys. Lett. 476, 97 (2009).
  • Roos and Andersson (1995) B. O. Roos and K. Andersson, Chem. Phys. Lett. 245, 215 (1995).
  • Fabiano et al. (2008) E. Fabiano, T. W. Keal,  and W. Thiel, Chem. Phys. 349, 334 (2008).
  • Koch et al. (1990) H. Koch, H. J. A. Jensen, P. Jørgensen, T. Helgaker, G. E. Scuseria,  and H. F. Schaefer III, J. Chem. Phys. 92, 4924 (1990).
  • Handy and Schaefer (1984) N. C. Handy and H. F. Schaefer, J. Chem. Phys. 81, 5031 (1984).
  • Győrffy et al. (2013) W. Győrffy, T. Shiozaki, G. Knizia,  and H.-J. Werner, J. Chem. Phys. 138, 104104 (2013).
  • Fatehi and Subotnik (2012) S. Fatehi and J. E. Subotnik, J. Phys. Chem. Lett. 3, 2039 (2012).
  • Birkholz and Schlegel (2016) A. B. Birkholz and H. B. Schlegel, Theor. Chem. Acc. 135, 84 (2016).
  • (55) bagel, Brilliantly Advanced General Electronic-structure Library. under the GNU General Public License.
  • Barbatti et al. (2004) M. Barbatti, J. Paier,  and H. Lischka, J. Chem. Phys. 121, 11614 (2004).
  • Azizi et al. (2006) Z. Azizi, B. O. Roos,  and V. Veryazov, Phys. Chem. Chem. Phys. 8, 2727 (2006).
  • Gozem et al. (2012) S. Gozem, M. Huntress, I. Schapiro, R. Lindh, A. A. Granovsky, C. Angeli,  and M. Olivucci, J. Chem. Theory Comput. 8, 4069 (2012).
  • Gozem et al. (2014) S. Gozem, F. Melaccio, A. Valentini, M. Filatov, M. Huix-Rotllant, N. Ferré, L. M. Frutos, C. Angeli, A. I. Krylov, A. A. Granovsky, R. Lindh, ,  and M. Olivucci, J. Chem. Theory Comput. 10, 3074 (2014).
  • Quenneville and Martínez (2003) J. Quenneville and T. J. Martínez, J. Phys. Chem. A 107, 829 (2003).
  • Minezawa and Gordon (2011) N. Minezawa and M. S. Gordon, J. Phys. Chem. A. 115, 7901 (2011).
  • Ioffe and Granovsky (2013) I. N. Ioffe and A. A. Granovsky, J. Chem. Theory Comput. 9, 4973 (2013).
  • Lei et al. (2014) Y. Lei, L. Yu, B. Zhou, C. Zhu, Z. Wen,  and S. H. Lin, J. Phys. Chem. A 118, 9021 (2014).
  • Harabuchi et al. (2014) Y. Harabuchi, K. Keipert, F. Zahariev, T. Taketsugu,  and M. S. Gordon, J. Phys. Chem. A 118, 11987 (2014).
  • Waldeck (1991) D. H. Waldeck, Chem. Rev. 91, 415 (1991).
  • Tsien (1998) R. Y. Tsien, Annu. Rev. Biochem. 67, 509 (1998).
  • Weber et al. (1999) W. Weber, V. Helms, A. J. McCammon,  and P. W. Langhoff, Proc. Natl. Acad. Sci. U.S.A. 96, 6177 (1999).
  • Martin et al. (2004) M. E. Martin, F. Negri,  and M. Olivucci, J. Am. Chem. Soc. 126, 5452 (2004).
  • Olsen and Smith (2008) S. Olsen and S. C. Smith, J. Am. Chem. Soc. 130, 8677 (2008).
  • Park and Rhee (2016) J. W. Park and Y. M. Rhee, J. Am. Chem. Soc. 138, 13619 (2016).
  • van Thor (2009) J. J. van Thor, Chem. Soc. Rev. 38, 2935 (2009).
  • Polyakov et al. (2010) I. V. Polyakov, B. L. Grigorenko, E. M. Epifanovsky, A. I. Krylov,  and A. V. Nemukhin, J. Chem. Theory Comput. 6, 2377 (2010).
  • Acharya et al. (2017) A. Acharya, A. M. Bogdanov, B. L. Grigorenko, K. B. Bravaya, A. V. Nemukhin, K. A. Lukyanov,  and A. I. Krylov, Chem. Rev. 117, 758 (2017).
  • Morozov and Groenhof (2015) D. Morozov and G. Groenhof, Angew. Chem. Int. Ed. 55, 576 (2015).
  • Minezawa and Gordon (2012) N. Minezawa and M. S. Gordon, J. Chem. Phys. 137, 034116 (2012).
  • Zhang et al. (2014) Q. Zhang, X. Chen, G. Cui, W.-H. Fang,  and W. Thiel, Angew. Chem. Int. Ed. 53, 8649 (2014).
  • Mandal et al. (2004) D. Mandal, T. Tahara,  and S. R. Meech, J. Phys. Chem. B 108, 1102 (2004).
  • Olsen (2010) S. Olsen, J. Chem. Theory Comput. 6, 1089 (2010).
  • Olsen and McKenzie (2009) S. Olsen and R. H. McKenzie, J. Chem. Phys. 130, 184302 (2009).
  • Galloy and Lorquet (1977) C. Galloy and J. C. Lorquet, J. Chem. Phys. 67, 4672 (1977).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description