Single-State Single-Reference and Multistate Multireference Zeroth

Article Cite This: J. Chem. Theory Comput. 2019, 15, 3960−3973

pubs.acs.org/JCTC

Single-State Single-Reference and Multistate Multireference ZerothOrder Hamiltonians in MS-CASPT2 and Conical Intersections Jae Woo Park* Department of Chemistry, Chungbuk National University (CBNU), Cheongju 28644, Korea

Downloaded via NOTTINGHAM TRENT UNIV on August 12, 2019 at 16:59:24 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: Multistate complete active space second-order perturbation theory (MSCASPT2) is one of the most successful quantum chemical methods for both static and dynamical correlations in photochemistry. In the literature, there are two definitions of the zeroth-order Hamiltonian depending on the form of the one-electron operator: multistate multireference MS-CASPT2 (MS-MR-MS-CASPT2) and single-state single-reference MSCASPT2 (SS-SR-MS-CASPT2). Here, we implement an analytical gradient and derivative coupling for SS-SR-MS-CASPT2 and test SS-SR-MS-CASPT2 against MS-MR-(X)MSCASPT2 in optimizing the molecular geometry critical points [the most stable geometries and minimum energy conical intersections (MECIs)] of a rhodopsin protein chromophore model (PSB3) and a green fluorescent protein chromophore model. In both cases, the MECIs in MS-MR-XMS-CASPT2 tend to have bridge hydrogen bonds that are more out-of-plane than those in SS-SR-MS-CASPT2, and for PSB3, the topology of potential energy surfaces (PESs) near the conical intersections is different. This result implies that caution is needed in analyzing the simulation results with different zeroth-order Hamiltonians and state-averaging schemes. The MS-MR-XMSCASPT2 theory yields smooth PESs near the MECIs in all the tested cases, while the SS-SR-MS-CASPT2 theory does not. Therefore, we recommend using MS-MR-XMS-CASPT2 in conical intersection simulations, which require smooth PESs.

1. INTRODUCTION Multireference theories, which use a set of electronic configurations as references, are widely used to study a variety of excited states.1 Among many quantum chemical methods [for example, configuration interactions (CI)2 and coupled cluster (CC)3], perturbation theory (PT) has been a powerful tool for including both static and dynamical correlation due to its accuracy and efficiency in providing balanced descriptions of both static and dynamical correlations. In the last three decades, approaches including the multireference second-order Møller−Plesset perturbation theory (MR-MP2) by Hirao,4−6 the complete active space second-order perturbation theory (CASPT2) by Roos and coworkers,7,8 the multireference openshell perturbation theory (MROPT) by Davidson and collaborators,9 the second-order n-electron valence state perturbation theory (NEVPT2) by Angeli, Cimiraglia, and coworkers,10−12 the generalized Van Vleck perturbation theory by Hoffmann and coworkers,13 and the static-then-dynamicthen-static perturbation theory (SDSPT) by Liu, Hoffmann, and coworkers14,15 have been developed. The multistate (diagonalize-perturb-diagonalize) formulations for these theories were developed to include the couplings between the states. For example, the multistate versions of MR-MP2, CASPT2, and NEVPT2 are secondorder multiconfiguration quasidegenerate perturbation theory (MCQDPT2),16,17 multistate CASPT2 (MS-CASPT2),18 and quasidegenerate NEVPT2 (QD-NEVPT2),10,19,20 respectively. Each state-averaged complete active space self-consistent field (SA-CASSCF) reference state, which diagonalizes the electronic Hamiltonian in the complete active space config© 2019 American Chemical Society

uration interaction (CASCI) space, is approximately corrected to the first order with the perturbation theory, and the approximate effective Hamiltonian is then constructed and diagonalized to arrive at the final states. The CASPT2 method is one of the most popular multireference perturbation theories, mainly because of its efficiency owing to internal contraction.7,8,21,22 In the literature, there are two schemes for MS-CASPT2, depending on the form of the zeroth-order Hamiltonian. The first implementation of MS-CASPT2 is called single-state singlereference MS-CASPT2 (SS-SR-MS-CASPT2).18 In SS-SR-MSCASPT2, the partitioning of the Hamiltonian is state specific, M M M M (0) ĤM = Pf̂ ̂ P ̂ + Q̂ f ̂ Q̂

(1)

based on the multipartitioning of the Hamiltonian by Zaitsevskii and Malrieu.23 Here, fM̂ is the Fock operator for state M, whose wave function is |M⟩, and P̂ and Q̂ M are the projection operators to the reference space and the interacting space for state M, respectively. The interacting space is statespecific. This partitioning ensures that the diagonal elements of the effective Hamiltonian are the same as the energies from the single-state CASPT2 (SS-CASPT2) calculations. The multistate multireference MS-CASPT2 (MS-MR-MSCASPT2) theory uses the same zeroth-order Hamiltonian for all states, as follows: SA SA (0) Ĥ = Pf̂ ̂ P ̂ + Q̂ f ̂ Q̂

(2)

Received: January 24, 2019 Published: May 16, 2019 3960

DOI: 10.1021/acs.jctc.9b00067 J. Chem. Theory Comput. 2019, 15, 3960−3973

Article

Journal of Chemical Theory and Computation where fŜ A is the state-averaged Fock operator. The interacting space can be state-specific (with the single-state singlereference (SS-SR) contraction scheme or without the “extension24”) or the same for all the reference states (in the MS-MR-XMS-CASPT2 theory with the multistate multireference (MS-MR) contraction scheme). The use of the state-averaged Fock operator here not only yields mathematical simplicity but also allows a definition of the zeroth-order Hamiltonian that is invariant with the rotations among the states. The definition of the “invariant” zeroth-order Hamiltonian was first introduced in Granovsky’s pioneering work in the context of the MCQDPT2 method. The resulting extended MCQDPT2 (XMCQDPT2) theory24 was implemented in the FIREFLY program suite.25 Then, Shiozaki et al. implemented the extended MS-CASPT2 (XMS-CASPT2) theory, which uses the invariant form of the zeroth-order Hamiltonian in MS-CASPT2, as was suggested by Granovsky.24,26 The XMCQDPT2 and XMS-CASPT2 methods have been successfully applied for near-degeneracy, particularly in removing the “bumps” on the potential energy surfaces (PESs) near conical intersections in MCQDPT2 and MSCASPT2.24,26−28 Current implementations of the partially or fully internally contracted MS-CASPT2 analytical gradient and derivative couplings in the program packages MOLPRO29 and BAGEL30 are based on MS-MR-(X)MS-CASPT2, 26,27,31−35 while MOLCAS,36 one of the most robust, efficient, and popular multireference quantum chemistry programs, uses SS-SR-MSCASPT2 (with finite difference nuclear gradients; note that MOLPRO also implements the SS-SR-MS-CASPT2 energy). The recent implementation of XMS-CASPT2 with a density matrix renormalization group (DMRG) reference function is based on the MS-MR-XMS-CASPT2 theory.37 Researchers employ both schemes for geometry optimizations and molecular dynamics simulations.1,35,38−44 The selection of the zeroth-order Hamiltonian affects the computational results, particularly in conical intersection dynamics, in which energy and wave function differences between the states play important roles. In a recent report by the author and Shiozaki, the selection of the state-averaging and multistate schemes for obtaining the minimum energy conical intersection (MECI) changed the relative energies of MECIs when comparing the CASSCF and CASPT2 conical intersections, as the perturbative energy differences between different states strongly depend on the form of the perturbation operator.45 In this study, we investigate the effect of the selection of the zeroth-order Hamiltonian in MSCASPT2. More specifically, we implement analytical nuclear gradient and derivative coupling within the framework of SSSR-MS-CASPT2, optimize the conical intersections of a rhodopsin protein (RP) chromophore model (PSB3) and a green fluorescent protein (GFP) chromophore model [parahydroxybenzylidene-imidazolinone (pHBI)] using SS-SR-MSCASPT2 and MS-MR-(X)MS-CASPT2, and discuss their differences.

denote virtual orbitals; and x, y, z, and w denote general orbitals. 2.1. MS-CASPT2. In the CASPT2 theory, the Hylleraas functional, (2) (1) ̂ (0) (0) (1) ̂ EM = 2⟨Φ(1) M |H |M ⟩ + ⟨Φ M |HM − EM + Eshift |Φ M ⟩

(3)

is minimized (i.e., the stationary point is found), where Ĥ (0) M is the zeroth-order Hamiltonian, E(0) is the zeroth-order energy, M Φ(1) M is the perturbation-theory based first-order wave function, and Eshift is the level shift to avoid the intruder state problem.46 The first-order wave function is parametrized with the perturbative operator T̂ M as ̂ |Φ(1) M ⟩ = TM |M ⟩ =

∑ TM ,Ω|ΩM ⟩ Ω

(4)

where |ΩM⟩ is a basis function generated by applying an excitation operator Ê Ω to the reference state |M⟩. Here, we employed the so-called SS-SR contraction scheme18 to parametrize the first-order wave function, which is the contraction scheme in the original MS-CASPT2 theory based on the multipartitioning technique.23 Hereafter, the terms MS-MR and SS-SR describe the form of the zeroth-order Hamiltonian, not the contraction scheme of the first-order wave function.27,34,37 Then, the amplitude equation, σΩM = ⟨ΩM |Ĥ |M ⟩ +

(0)

∑ ⟨ΩM|ĤM

(0) − EM + Eshift|Ω′M ⟩TM , Ω′ = 0

Ω′M

(5)

is solved to obtain the amplitudes. The summation over the configurations |ΩM ′ ⟩ is needed because the zeroth-order Hamiltonian is nondiagonal. When the amplitudes for all the reference states are obtained (i.e., the amplitude equations for all the states are solved), the effective Hamiltonian, whose elements are 1 eff ̂ † Ĥ |N ⟩) ̂ N̂ |N ⟩ + ⟨M |TM = ⟨M |Ĥ |N ⟩ + (⟨M |HT HMN 2 ̂ † TM̂ |M ⟩ − δMN Eshift⟨M |TM (6) is constructed. The effective Hamiltonian is then diagonalized to obtain the MS-CASPT2 states, eff RNP = ∑ RMPEP ∑ HMN N

(7)

P

2.2. Definition of the Zeroth-Order Hamiltonian. The ambiguity in the MS-CASPT2 formalism remains in the definition of the zeroth-order Hamiltonian. The MS-MR-MSCASPT2 theory uses the same zeroth-order Hamiltonian for all states,26 while SS-SR-MS-CASPT2 uses zeroth-order Hamiltonians that are specific to the reference states.18 We describe the working expressions for both theories in the following sections. 2.2.1. MS-MR-(X)MS-CASPT2. In MS-MR-MS-CASPT2, the zeroth-order Hamiltonian is defined as SA SA (0) Ĥ = Pf̂ ̂ P ̂ + Q̂ f ̂ Q̂

2. THEORY In this section, we review the MS-CASPT2 theory with both the MS-MR and SS-SR zeroth-order Hamiltonians. In the following, the indices M and N denote reference (SACASSCF) states; P and Q denote (X)MS-CASPT2 states; I, J, and K denote Slater determinants; i, j, k, and l denote closed shell orbitals; r, s, t, and u denote active orbitals; a, b, c, and d

(8)

where the projection operators P̂ and Q̂ are P̂ =

∑ |M⟩⟨M| M

Q̂ = 1 − P ̂ 3961

(9) (10) DOI: 10.1021/acs.jctc.9b00067 J. Chem. Theory Comput. 2019, 15, 3960−3973

Article

Journal of Chemical Theory and Computation and the interacting space is spanned by doubly excited configurations. The state-averaged Fock operator is f̂

SA

=

∑

̂ f xySA Exy

=

∑ [h + g ( d

xy SA dxy

SA

(0) ĤM =

f̂

(0) Ĥ =

SA

SA |M ⟩⟨M | + Q̂ f ̂ Q̂

|M ⟩|N ⟩ (16)

N

Ĥ

=

∑ |M̅ ⟩⟨M̅ |f ̂ M

SA

|M̅ ⟩⟨M̅ | + Q̂ f ̂

SA

∑ |M⟩UMM̅ M

∑ RNP⟨ΩM |Ĥ |N ⟩ N

− 2EshiftRMP[⟨ΩM |TM̂ |M ⟩]

(26)

is solved to obtain such a set. Then, the total Lagrangian is

Q̂

1 1 tr[Z(A − A†)] − tr[X(C†SC − 1)] 2 2 ÉÑ ÄÅ ÑÑ ÅÅ 1 Ñ ÅÅ ref + ∑ WN ÅÅ∑ zI , N ⟨I |Ĥ − EN |N ⟩ − xN (⟨N |N ⟩ − 1)ÑÑÑ ÑÑ ÅÅ 2 ÑÖ ÅÇ I N

LP = L PT2,P +

(18)

The amplitude equations are then solved for each “rotated” reference state, |M̅ ⟩ =

(25)

is stationary with a suitable set of λ, where σΩM is the residual in eq 5. The so-called λ-equation,

where U is a unitary transformation matrix. Then, the zerothorder Hamiltonian can be rewritten as (0)

M

(17)

MN

(24)

∑ λM ,ΩσΩ

(0) (0) ̂ 0 = ⟨ΩM |(ĤM − EM )λM |M ⟩ +

(0) |N ⟩UNN̅ = δ MN ̅ ̅ E M̅

Q̂

M ,Ω

while f is not diagonal on the basis of M. To include the offdiagonal elements while maintaining the reference function as an eigenfunction, one diagonalizes the state-averaged Fock operator as follows: SA

O

∑ Q̂ f ̂ O

L PT2,P = EP +

Ŝ A

∑ UMM̅ ⟨M|f ̂

O

|N ⟩⟨N |f ̂ |M ⟩⟨M | +

which is physically the same as eq 15 when we sum over O and the reference weights are the same. One cannot define the state-specific zeroth-order Hamiltonian without losing the simplicity of the MS-CASPT2 theory (i.e., solving the CASPT2 equation for each state and mixing these states by constructing the effective Hamiltonian), as the Fock operator for each state should be separately diagonalized to include all the off-diagonal terms. We tried to define an approximately invariant zerothorder Hamiltonian to incorporate the XMS procedure into SSSR-MS-CASPT2, which we present in the Supporting Information. 2.3. Nuclear Gradient Theory. Now, we address the nuclear gradient theory for SS-SR-MS-CASPT2. (The nuclear gradient theory for MS-MR-(X)MS-CASPT2 is presented elsewhere.26,32,34) The MS-CASPT2 energy is not stationary with respect to the amplitudes. The CASPT2 Lagrangian,26,34

(15)

∑ ⟨N |f ̂

∑ MNO

This expression is invariant with respect to the rotations among the states, assuming that all the weights for the references in the CASSCF calculations are the same. When we use this operator as the zeroth-order Hamiltonian, the reference function will not be an eigenfunction of the zerothorder Hamiltonian, i.e., (0) Ĥ |M ⟩ =

(23)

The Fock operator used in the zeroth-order Hamiltonian now varies with respect to the rotations among the states. The invariant form is

and this definition ensures that |M⟩ is an eigenfunction of Ĥ (0). To make the zeroth-order Hamiltonian invariant with respect to the rotations of the reference functions, Granovsky24 included the off-diagonal contributions in the zeroth-order Hamiltonian as follows:

MN

(22)

M

(14)

SA

(21)

xy

(0) EM = ⟨M |f ̂ |M ⟩

SA |M ⟩⟨M | + Q̂ f ̂ Q̂

M

∑ |N ⟩⟨N |f ̂

∑ fxyM Exŷ = ∑ [h + g(dM)]xy Exŷ

This definition ensures that the zeroth-order wave function |M⟩ is an eigenfunction of Ĥ (0) M with an eigenvalue

and WM is the weight in the SA-CASSCF calculation. Conventional MS-CASPT2 includes only the diagonal elements of the first term in eq 8, as follows:

(0) Ĥ =

=

̂ |M ⟩ dxyM = ⟨M |Exy

(13)

SA

M

xy

where h is the one-electron integral, dSA is the state-averaged one-electron density matrix, ÄÅ ÉÑ Å 1 SA SA SA Ñ [g(dSA)]xy = ∑ ÅÅÅÅ(xy|zw)dzw − (xw|zy)(dzw + dwz )ÑÑÑÑ ÅÇ ÑÑÖ 4 zw Å

∑ |M⟩⟨M|f ̂

(20)

(12)

M

(0) Ĥ =

M M M |N ⟩⟨N | + Q̂ f ̂ Q̂

where the Fock operator is state-specific, (11)

̂ |M ⟩ = ∑ WM⟨M |Exy

M

N

̂ )]xy Exy

xy

∑ |N ⟩⟨N |f ̂

closed frozen

+

(19)

∑ ∑ i

resulting in the XMS-CASPT2 theory. We note that the state weights should be the same for each state in XMS-CASPT2, and using the SS-SR internal contraction scheme instead of MS-MR breaks the invariance.37 2.2.2. SS-SR-MS-CASPT2. In SS-SR-MS-CASPT2, the zeroth-order Hamiltonian for state M is

j

c

zijccf ijSA c

(27)

Here, Z, z, X, and x are the Lagrangian multipliers for the CASSCF orbital convergence, full CI convergence, orbital orthonormality, and CI vector orthonormality, respectively. The condition for inactive frozen core orbitals is specified in the last term, whose multiplier is zc. 3962


Journal of Chemical Theory and Computation The Lagrangian multipliers are calculated by solving the socalled Z-vector equation.31 The source terms for the Z-vector equation are the derivatives of the CASPT2 Lagrangian with respect to the SA-CASSCF orbital rotation coefficient κxy and the SA-CASSCF CI coefficient cI,M, Yrs =

(2), M dxy

NM = ⟨M |λM̂ TM̂ |M ⟩

(28)

∂cI , M

(29)

These terms can be evaluated by recasting the CASPT2 Lagrangian into a compact form with the molecular integrals and reduced density matrices (RDMs) as follows: L PT2, P = tr(hd) +

∑ tr[g(d(0), M)d(2), M] + ∑ tr(KklDlk ) M

+

kl

jst

where we introduce a matrix form of the two-electron integrals and second-order density matrices, (31) (32)

yI , M =

There are additional terms associated with the level shifts46 in the compact form of the Lagrangian, eq 30. These terms do not contribute to the orbital derivative and thus are omitted. The total one-electron and two-electron density matrices d and D are (1)

+d

(2)

=

ÉÑ ÑÑ kl lk Ñ + ∑ K D ÑÑÑ ÑÑ kl ÑÖ

(44)

(45)

(0), M ̂ |M ⟩ dxy = ⟨M |Exy

†

†

†

†

†

(0)

(0) + ⟨M |λM̂ (ĤM − EM + Eshift)TM̂ |I ⟩ † (0) (0) + ⟨I |λM̂ (ĤM − EM + Eshift)TM̂ |M ⟩

+ 2 ∑ ∑ ⟨I |Erŝ |M ⟩[g(d(2), M) − NM f M]rs M

rs

(46)

(36)

ln

†

∑ RMPRNP(2⟨I |Ĥ |N ⟩ + ⟨I |TM̂ Ĥ |N ⟩ + ⟨N |TM̂ Ĥ |I ⟩)

+ ⟨M |λM̂ Ĥ |I ⟩ + ⟨I |λM̂ Ĥ |M ⟩

(35)

̂ |N ⟩ = ∑ RLPRNP⟨L|Exy

∂cI , M

†

(34)

+D

∂L PT2, P

2 ̂ Ĥ |I ⟩ + ⟨I |TM ̂ Ĥ |M ⟩) (⟨M |TM + RMP

where the zeroth- and first-order contributions are

(0) dxy

M

M

N

(33)

+d

(1)

D=D

∑ g(d(0), M)d(2), M

∑ g(d(2), M)d(0), M

The CI derivative with respect to the CI coefficient cI,M is

Diljk = Dij , kl

(0)

∑ Dstij (rs|tj)

ÄÅ ÅÅ 1 Å Yra = ÅÅÅhd + ÅÅ 2 ÅÇ

(30)

Kxyzw = (xz|yw)

(43)

The orbital derivatives are then ÄÅ ÅÅ 1 Å Yri = ÅÅÅhd + ∑ g(d(0), M)d(2), M + ÅÅ 2 ÅÇ M ÉÑ ÑÑ kl lk Ñ + ∑ K D ÑÑÑ ÑÑÑÖ kl ri

∂L PT2, P

d=d

(42)

†

∂κxy

(0)

l (2), M o ̂ |M ⟩ x , y ∈ r , s o − NM⟨M |Exy ̅ o o dxy =m o o o d ̅ (2), M o otherwise o n xy

and the correlated norm is

∂L PT2, P

yI , M =

Article

c

(0) Dxyzw =

̂ |N ⟩ ∑ RLPRNP⟨L|Exyzw

(37)

ln (1) dxy

=

The constraint z is calculated before solving the Z-vector equation and added to the source terms as follows: zijcc = −

† ̂ |N ⟩ RLPRNP⟨L|TL̂ Exy

∑ ln

+

L (1) Dxyzw =

(38)

̂ |N ⟩ ∑ RLPRNP⟨L|TL̂ †Exyzw ln

+

†

† 2 ̂ |L⟩ + ⟨L|λL̂ Exyzw ̂ |L ⟩ ⟨L|TL̂ Exyzw ∑ RLP L

∑ d(2),M (40)

M

where (2), M

dxy ̅

†

̂ Exy ̂ λM̂ |M ⟩ = ⟨M |TM

(47)

3. RESULTS AND DISCUSSION In this section, we present the optimization of the MECIs of a retinal model chromophore, penta-2,4-dieniminium cation (PSB3), and 4-para-hydroxybenzylidene-imidazolin-5-one (pHBI) using SA-CASSCF, MS-MR-(X)MS-CASPT2, and SS-SR-MS-CASPT2. We used the so-called SS-SR contraction scheme34 with a level shift46 of 0.2 Eh. We searched for the MECIs using a gradient projection method.47 We carried out the calculations using a modified version of the program

(39)

and the second-order contribution is d(2) =

fiiSA − f jSA c c j

The Z-vector equation is then solved with these source terms using the method described in ref 31. The extension to derivative coupling is straightforward.27 Although our program is not yet fully optimized, the computational cost of the SS-SRMS-CASPT2 gradient is almost the same as that of the usual MS-MR-(X)MS-CASPT2 gradient.

†

† 2 ̂ |L⟩ + ⟨L|λL̂ Exy ̂ |L ⟩ ⟨L|TL̂ Exy ∑ RLP

Yijc − Yjc i

(41) 3963


Article

Journal of Chemical Theory and Computation package BAGEL.30 We used the cc-pVDZ basis set and its corresponding density-fitting basis set. Note that we set the S0 energy at the Franck−Condon (FC) point to zero. For the sake of brevity, we describe a PES as “more stabilized” when the PT2 correction stabilizes the PES more than the S0 state at the FC point. Below, the discussion of MS-MR-(X)MSCASPT2 applies to both MS-MR-MS-CASPT2 and MS-MRXMS-CASPT2. 3.1. Conical Intersection of PSB3. PSB3 (Figure 1) is one of the most widely used computational models of the

The representative geometrical parameters for each conical intersection are shown in Table 4. The main difference between the SS-SR and MS-MR conical intersections is the twisting angle of the C2C3 bond, or the “reaction coordinate” of the photoisomerization,48,53−55 while the CASSCF and CASPT2 MECIs differ mainly in the lengths of N−C1 and C1−C2 bonds, which compose the “bond length alternation (BLA)” coordinate. The BLA coordinate of the SSSR-MS-CASPT2 MECI is between those of the MS-MR-XMSCASPT2 and SA-CASSCF MECIs. The MS-MR MECI occurs at a geometry that is more twisted than that of the SS-SR MECI and is more displaced from the SA-CASSCF MECI. The SS-SR-MS-CASPT2 twisting angle compares favorably with the previous MRCISD (−92.6°)64 and XMCQDPT2 (−92.2°) results.53 This displacement from the SA-CASSCF MECI is due to the differential correlation effect. To inspect the differential correlation effect for PSB3, we computed the potential energies at the interpolated coordinates between the MECIs by SA-CASSCF, SS-SR-MS-CASPT2, and MS-MR-XMSCASPT2. The resulting potential energy scan is shown in Figure 3. In terms of the diabatic states,45,48,59 the charge transfer (1Bu) state is more stabilized than the diradicaloid (1Ag) state by PT2 correction in both cases. However, the amounts of this differential correlation effect are approximately 0.7 eV (SS-SR-XMS-CASPT2) and 1.0 eV (MS-MR-XMSCASPT2). This difference makes the distance between the SSSR-MS-CASPT2 and SA-CASSCF MECIs smaller than the distance between the MS-MR-XMS-CASPT2 and SA-CASSCF MECIs. Indeed, the SS-SR-MS-CASPT2 surfaces cross between the SA-CASSCF and MS-MR-XMS-CASPT2 MECIs (Figure 3b). The MS-MR-MS-CASPT2 theory yields PESs with nonphysical “bumps” near the conical intersections for PSB3. 2 7 , 2 8 , 5 3 , 5 5 Using the extended theories, e.g., XMCQDPT253,55 and XMS-CASPT2,27,28 successfully removes these bumps near the conical intersections. We inspected the behavior of MS-MR-XMS-CASPT2 and SS-SRMS-CASPT2 by scanning the PESs near the MECIs on the two-dimensional coordinate that consists of the gradient difference vector g and the interstate coupling vector h in the range of [−0.4 Å, +0.4 Å]. The resulting PES contour plots are displayed in Figure 4, while we show these vectors in Figure 5. As expected, the MS-MR-XMS-CASPT2 PES does not exhibit unphysical behavior near the MECI. The SS-SRMS-CASPT2 PES shows sawtooth-like bumps near the MECI, as expected. One of the important properties of conical intersections is their topology.65 The peaked and sloped topologies are possible, depending on the relative orientation of the crossing surfaces. The SS-SR-MS-CASPT2 conical intersections give a peaked topology, which can yield two or more relaxation pathways to

Figure 1. Chemical structure of PSB3.

rhodopsin chromophore.27,28,40,48−58 A recent benchmark study has shown that XMS-CASPT2 can yield the correct topology of the conical intersections of PSB3, comparable to that of MRCISD+Q.28 The active space of (6e,6o), which includes all the π- and π*-orbitals, was used. We included the three lowest singlet states in the SA-CASSCF and MSCASPT2 calculations. The mixing between the S2 state and the rest is not dominant in PSB3, as the S1−S2 gap energy is larger than 1 eV at the FC point. However, this mixing becomes more severe for a larger RP chromophore model or in the presence of solvents59 or electron-withdrawing functional groups.60 Although including the third root is necessary for studies of the larger RP chromophore models, we selected the three-state model in this section for consistency. We additionally optimized the FC point using a two-state averaging scheme and found that the vertical excitation energies do not change with respect to the selection of the state-averaging scheme (Table 2). Twisting around the C2C3 bond is the lowest-energy channel for photoisomerization in PSB3.40,49,53,55 We optimized the structure of the ground-state equilibrium geometry and the MECI that corresponds to C2C3 bond twisting using SA-CASSCF, MS-MR-(X)MS-CASPT2, and SS-SR-MSCASPT2. The energies of the ground and excited states at the FC point and MECI are shown in Table 1. All the energies at the MECI are lower than the S1 energy at the FC point. This finding is consistent with previous results (see Tables 2 and 3), including a series of benchmark studies by Gozem, Olivucci, and coworkers28,53,55,61−63 and the CASPT2 dynamics trajectories that yielded rapid decay to the ground statistics (∼100 fs) with an initial velocity at room temperature.44 The MECIs optimized using MS-MR-XMS-CASPT2 and SS-SR-MS-CASPT2 are shown in Figure 2. The root-meansquare (rms) distances to the SA-CASSCF MECI geometries are 0.086 Å (SS-SR) and 0.132 Å (MS-MR) after alignment.

Table 1. Energies of the S0, S1, and S2 States of PSB3 at the FC Point and the MECIs with Respect to the S0 Energy at the FC Point (eV) Calculated by SA-CASSCF, SS-SR-MS-CASPT2, MS-MR-MS-CASPT2, and MS-MR-XMS-CASPT2a SA-CASSCF EFC ECI ECI−EFC

SS-SR-MS-CASPT2

MS-MR-MS-CASPT2

MS-MR-XMS-CASPT2

S0

S1

S2

S0

S1

S2

S0

S1

S2

S0

S1

S2

0.00 3.20

5.05 3.20 1.86

5.73 6.60

0.00 2.41

3.99 2.41 1.58

5.15 5.68

0.00 2.53

4.12 2.53 1.59

5.32 5.81

0.00 2.52

4.16 2.52 1.64

5.31 5.79

a

The S1 energy differences between the FC and MECI points are also shown. 3964


Article

Journal of Chemical Theory and Computation Table 2. Comparison of the Vertical Absorption Energy of PSB3 Calculated with Multireference Methods method

basis set

MS-MR-MS-CASPT2 MS-MR-MS-CASPT2 CASSCFc MRCISD+Q//CASSCF SS-SR-MS-CASPT2//CASSCF SS-SR-MS-CASPT2 (IPEA = 0.25)d//CASSCF MRCISD+Q//CASSCFc MRCISD SA-CASSCF MS-MR-MS-CASPT2 SA-CASSCF MS-MR-XMS-CASPT2 MS-MR-XMS-CASPT2 SS-SR-MS-CASPT2 MS-MR-XMS-CASPT2 SS-SR-MS-CASPT2

6-31G* 6-31G* 6-31G* 6-31G* 6-31G* 6-31G* 6-31G* 6-31+G** 6-31G 6-31G* 6-31G* 6-31G* cc-pVDZ cc-pVDZ cc-pVDZ cc-pVDZ

active space 6e, 6e, 6e, 6e, 6e, 6e, 6e, 4e, 6e, 6e, 6e, 6e, 6e, 6e, 6e, 6e,

6o 6o 6o 6o 6o 6o 6o 4o 6o 6o 6o 6o 6o 6o 6o 6o

Nstatea

ΔE (FC)b

ref

2 3 2 2 2 2 2 2 3 3 2 2 3 3 2 2

3.96 4.28 4.78 4.53 3.99 4.38 4.40 4.47 4.91 4.26 4.95 4.19 4.16 3.99 4.11 4.02

58 49 53, 62 28, 53, 63 53 53 62 64 40 40 28 28 this work this work this work this work

Number of states included in the SA-CASSCF and MS-CASPT2 calculations. bS0−S1 energy gap at the Franck−Condon point. ccis conformer. The IPEA shift parameter92 for the systematic correction of CASPT2. This parameter is zero for the CASPT2 results if no value is mentioned.

a

d

Table 3. MECI Energies of the PSB3 Calculated with Multireference Methods method

basis set

SS-SR-MS-CASPT2 MS-MR-MS-CASPT2 MS-MR-MS-CASPT2 SA-CASSCF SS-SR-MS-CASPT2, QD-NEVPT2, XMCQDPT2//CASSCF MRCISD SA-CASSCF MS-MR-MS-CASPT2 SA-CASSCF MS-MR-XMS-CASPT2 SS-SR-MS-CASPT2

6-31G* 6-31G* 6-31G* 6-31G* 6-31G* 6-31+G** 6-31G 6-31G* ANO-RCC cc-pVDZ cc-pVDZ

active space 6e, 6e, 6e, 6e, 6e, 4e, 6e, 6e, 6e, 6e, 6e,

6o 6o 6o 6o 6o 4o 6o 6o 6o 6o 6o

Nstatea

Eb

ΔE (CI)c

ref

2 2 3 2 2 2 3 3 2 3 3

* 2.04 2.6 2.59 *e 2.38 2.87 2.61 *f 2.52 2.41

0.04

57 58 49 53 53 64 40 40 89 this work this work

d

*e

a

Number of states included in the SA-CASSCF and MS-CASPT2 calculations. bS1 energy of CI relative to the S0 energy at the Franck−Condon point. cS0−S1 gap energy at the optimized MECI due to the differential correlation effect. dNot presented. eSee Figure 2 in ref 53. fThe S0 energy was not reported. The absolute energy is −248.341553 Eh.

the ground state. On the other hand, the MS-MR-XMSCASPT2 and SA-CASSCF conical intersections give a sloped topology with a single relaxation pathway. We can see this difference in conical intersection topology more clearly in Figure 3b. The SA-CASSCF, SS-SR-MS-CASPT2, and “MSMR-XMS-CASPT2 MECIs are located sequentially along the BLA coordinate. The SA-CASSCF conical intersection is sloped; the SS-SR-MS-CASPT2 conical intersection is at the point at which two PESs have a peaked topology, and the MSMR-XMS-CASPT2 conical intersection is again sloped. The branching plane vectors (Figure 5) show a behavior consistent with this result: the g and h vectors at the SA-CASSCF MECI correspond to the twisting along the C2C3 bond (reaction coordinate) and the stretching (or contraction) of the C1−N bond (BLA coordinate), respectively, as in work by Gozem and coworkers.53 The SS-SR-MS-CASPT2 branching plane vectors are less relevant to the BLA coordinate, and both vectors are related to the reaction coordinate. Then, again, the MS-MR-MS-CASPT2 branching plane vectors are similar to those at the SA-CASSCF MECI, although the sign of the h vector is reversed. We also note that with the MS-MR-XMSCASPT2 method, these vectors are interchanged (i.e., the g

Figure 2. MECI geometries of PSB3 optimized by (a) SS-SR-MSCASPT2 and (b) MS-MR-XMS-CASPT2. For comparison, the MECIs optimized by SA-CASSCF are also shown in gray.

Table 4. Representative Geometrical Parameters in the Optimized MECIsa

C1−C2−C3−C4 torsion (deg) N−C1 bond (Å) C1−C2 bond (Å)

SACASSCF

SS-SR-MSCASPT2

MS-MR-MSCASPT2

MS-MRXMSCASPT2

−94.96

−92.05

−96.85

−97.88

1.37 1.33

1.34 1.39

1.31 1.45

1.31 1.44

a

The atomic labels are shown in Figure 1.

3965


Article

Journal of Chemical Theory and Computation

Figure 4. PES scans near the MECIs of PSB3. The vectors g and h denote the gradient difference and interstate coupling vectors, respectively. For MS-MR-(X)MS-CASPT2, g and h are interchanged for visual clarity. The vectors g and h are shown in Figure 5.

Figure 3. PESs on the interpolated coordinates between the (a) SACASSCF and SS-SR-MS-CASPT2 MECIs, (b) SA-CASSCF and MSMR-XMS-CASPT2 MECIs, and (c) SS-SR-MS-CASPT2 and MSMR-XMS-CASPT2 MECIs. The SS-SR-MS-CASPT2 (red), MS-MRXMS-CASPT2 (blue), and SA-CASSCF (gray) potential energies are shown.

and h vectors correspond to the BLA and reaction coordinates, respectively) due to the rotation between the reference states. Finally, we note that previous works with a two-state averaging scheme showed peaked MECIs,28,44,53,54 differing from our and others’ three-state-averaged results.40,49,59 These observations imply that the nature of the decays (e.g., the approach to the conical intersection and the decays in the dynamics simulations) might differ by the selection of the state-averaging scheme and the zeroth-order Hamiltonian, and thus, these schemes should be chosen carefully when setting up dynamics simulations. 3.2. Conical Intersections of pHBI. GFP is a widely used imaging probe for biological systems.66−76 Notably, the isolated GFP chromophore undergoes twisting followed by rapid nonadiabatic decays when photoexcited, resulting in a low fluorescence quantum yield (lower than 10−4),68−74 and is a useful theoretical benchmark system for studying photochemical phenomena.77−80 There are two twisting channels around the imidazolinone (I)- and phenoxy (P)-bridge bonds. In this study, we optimized the molecular geometries of the FC point and MECIs of pHBI (Figure 6). There are many benchmark studies on the vertical excitation energies (Table 5) and MECIs (Table 6) of GFP chromophore models.27,64,69,75,80−89 Various active spaces and state-averaging schemes have been employed. We decided

Figure 5. Branching plane vectors of PSB3. The red and blue arrows indicate the g (gradient difference) and h (interstate coupling) vectors, respectively.

Figure 6. Chemical structure of pHBI.

to use a three-state-averaging scheme with (4e, 3o) active space.77 The resulting energies of the FC point and MECIs are shown in Table 7. An advantage of employing this scheme is the ability to analyze the results based on the diabatic states. These diabatic states “preserve” their electronic structure, and therefore, their chemical properties with respect to nuclear motion and the response to the perturbative correction remains similar on the PESs when the zeroth-order Hamiltonian remains the same. This assumption is only approximately correct in the SS-SR-MS-CASPT2 case, as the zeroth-order Hamiltonian differs in the adiabatic reference state. We can mix the three lowest adiabatic states of pHBI to form three diabatic states, |P⟩, |B⟩, and |I⟩, and these states 3966


Article


Table 5. Comparison of the Vertical Absorption Energy of the GFP Chromophore Models Calculated with Multireference Methods model

method

basis set

active space

Nstatea

ΔE (FC)b

pHBI pHBI pHBDIc pHBI pHBI pHBI pHBI pHBI pHBDIc pHBDIc pHBI pHBDIc pHBI pHBI pHBI pHBI pHBI pHBDI

SS-SR-MS-CASPT2 //CASSCF SS-CASPT2 //CASSCF MS-MR-MS-CASPT2 //CASSCF MS-MR-MS-CASPT2 SS-SR-MS-CASPT2 (IPEA = 0.25)d //BLYP SS-SR-MS-CASPT2 (IPEA = 0.25)d //BLYP SS-SR-MS-CASPT2 //BLYP/cc-pVTZ MS-MR-MS-CASPT2 //CASSCF aug-MCQDPT2//PBE0 XMCQDPT2 //PBE0 MRCISD MS-MR-XMS-CASPT2 SS-SR-MS-CASPT2 (IPEA = 0.25)d //MP2/cc-pVTZ SS-SR-MS-CASPT2 SS-SR-MS-CASPT2 MS-MR-XMS-CASPT2 MS-MR-XMS-CASPT2 exptl

6-31G* 6-31G DZP 6-31G** cc-pVTZ cc-pVTZ 6-31G* 6-31G (aug)-cc-pVDZe cc-pVDZ 6-31G** cc-pVDZ ANO-S cc-pVDZ cc-pVDZ cc-pVDZ cc-pVDZ

12e, 11o 2e, 2o 4e, 3o 2e, 2o 4e, 4o 14e, 14o 14e, 14o 6e,5o 16e, 14o 14e, 12o 2e, 2o 4e, 3o 16e, 14o 4e, 3o 4e, 3o 4e, 3o 4e, 3o

2 2 3 2 2 2 2 2 2 2 2 3 5 3 2 3 2

2.67 2.66 2.69 3.16 2.97 2.88 2.63 2.61 2.54 2.51 3.88 2.50 3.50 2.28 2.47 2.56 2.60 2.59

ref 75 81 80 82 83 83 83 84 85 69 64 27 86 this this this this 76

work work work work

Number of states included in the SA-CASSCF and/or MS-CASPT2 calculations. bS0−S1 energy gap at the Franck−Condon point. cparaHydroxybenzylidene-dimethyl-imidazolinone (pHBDI), which has two methyl groups attached to the imidazolinone ring. dThe IPEA shift parameter92 for the systematic correction of CASPT2. This parameter is zero for the CASPT2 results if no value is mentioned. eDiffusion functions added only to the oxygen atom. a

Table 6. MECI Energies of the GFP Chromophore Models Calculated with Multireference Methods model pHBI pHBI pHBDIf pHBDIf pHBDIf pHBDIf pHBDIf pHBI pHBI pHBI pHBI pHBI pHBI pHBI pHBDIf pHBDIf pHBI pHBI pHBI pHBI

method SS-SR-MS-CASPT2//CASSCF MS-MR-MS-CASPT2 CASSCF CASSCF MS-MR-MS-CASPT2//CASSCF MS-MR-MS-CASPT2//CASSCF CASSCF MS-MR-XMS-CASPT2//CASSCF CASSCF CASSCF CASSCF MRCISD MRCISD CASSCF MS-MR-XMS-CASPT2 MS-MR-XMS-CASPT2 MS-MR-XMS-CASPT2 MS-MR-XMS-CASPT2 SS-SR-MS-CASPT2 SS-SR-MS-CASPT2

basis set

active space

6-31G* 6-31G** DZP DZP DZP DZP cc-pVDZ 6-31G 6-31G** 6-31G* 6-31G* 6-31G** 6-31G** ANO-RCCh cc-pVDZ cc-pVDZ cc-pVDZ cc-pVDZ cc-pVDZ cc-pVDZ

12e, 11o 2e, 2o 4e, 3o 4e, 3o 4e, 3o 4e, 3o 12e, 11o 6e, 5o 2e, 2o 2e, 2o 2e, 2o 2e, 2o 2e, 2o 2e, 2o 4e, 3o 4e, 3o 4e, 3o 4e, 3o 4e, 3o 4e, 3o

Nstatea 2 2 3 3 3 3 2 2 2 3 3 2 2 2 3 3 3 3 3 3

conf.b e

HT P I P I P I HTe P I P I P I I P I P I P

Ec

ΔE (CI)d

2.74 2.56 2.82 3.38 2.52 3.14 2.99 2.34 *g 2.78 3.41 2.87 3.20 *i 2.45 2.94 2.40 2.88 2.00 2.44

0.20

0.65 0.85 0.38

ref 75 82 80 80 80 80 69 84 87 88 88 64 64 89 27 27 this this this this

work work work work

a

Number of states included in the SA-CASSCF and/or MS-CASPT2 calculations. bConformation of MECI. cS1 energy of CI relative to the S0 energy at the Franck−Condon point. dS0−S1 gap energy at the optimized MECI due to the differential correlation effect. eHula-twisted conical intersection. fpHBDI, which has two methyl groups attached to the imidazolinone ring. gThe S0 energy was not reported. The absolute energy is −641.303574 Eh. hWith double-ζ-plus-polarization contraction. iThe S0 energy was not reported. The absolute energy is −641.829002 Eh.

indicate that the charges are localized in phenoxy, bridge, and imidazolinone moieties, respectively.77,90 The adiabatic potential energies are transformed into the diabatic potential energy matrix D as follows: ij DPP DPI DPB yz jj zz j z † D = R VR = jjjj DIP DII DIB zzzz jj zz jD z k BP DBI DBB {

where V contains the adiabatic potential energies as diagonal elements, and the transformation matrix R is R = C−1(CC†)1/2

(49)

Here, C contains the CI vectors for the rotated configuration state functions (CSFs) |ppbi⟩, |iipb⟩, |bbpi⟩, in which the electrons occupy the localized orbitals p, b, and i formed from the active orbitals (see Figure 7). The orbital localization was performed using the Pipek−Mezey scheme.91 For example, at

(48) 3967


Article


Table 7. Energies of the S0, S1, and S2 States of pHBI at the FC Point and the MECIs with Respect to the S0 Energy at the FC Point (eV) Calculated by SA-CASSCF, SS-SR-MS-CASPT2, MS-MR-MS-CASPT2, and MS-MR-XMS-CASPT2a SA-CASSCF EFC ECI,P ECI,I ECI,P−EFC ECI,I−EFC

SS-SR-MS-CASPT2

MS-MR-MS-CASPT2

S0

S1

S2

S0

S1

S2

S0

0.00 3.44 2.95

4.14 3.44 2.95 −0.70 −1.19

6.01 7.02 6.72

0.00 2.44 2.00

2.28 2.44 2.00 0.16 −0.28

3.51 5.17 4.75

0.00 2.89b 2.38

S1 2.55 2.89b 2.38

MS-MR-XMS-CASPT2 S2

S0

S1

S2

3.57 5.08b 4.57

0.00 2.88 2.40

2.56 2.88 2.40 0.32 −0.16

3.57 5.07 4.58

−0.17

a

The S1 energy differences between the FC and MECI points are also shown. bA looser convergence criterion was used.

values with SS-SR-MS-CASPT2, despite recent discussions on the effectiveness of the IPEA shift.93 The energy of the P-twisted MECI is higher than that of the I-twisted MECI. This result agrees well with the previous studies (Table 6), which yield higher P-twisted MECI energy than I-twisted energy, ranging from 0.33 eV (MRCISD)64 to 0.63 eV (CASSCF).88 The major difference between the current results and most of the previous results is the relative energy of MECI compared to the S1 state at the FC point, which might be attributed to the minimal active space [(2e, 2o)] employed in many of these works. Previous studies with a larger active space yield conical intersections higher than the FC point, although they are based on the CASPT2//CASSCF scheme.75,80 We will reserve the relationship between the size of the active space and the relative energies of the MECI and FC point for future investigations. Now, let us discuss the structures of MECIs (Figure 8). The main difference among these MECI geometries is the out-of-

Figure 7. Localized active orbitals of pHBI at the Franck−Condon point optimized by SS-SR-MS-CASPT2. The Pipek−Mezey localization scheme was used.91 The orbital contours were generated with the computer program IboView.94,95

the FC point, the diabatic potential energy matrices DFC in eV are DFC,SS ‐ SR

ij1.26 0.96 0.75 yz jj zz = jjjj 0.96 1.36 − 0.74 zzzz jj zz k 0.75 − 0.74 3.18 {

(50)

ij1.38 1.11 0.74 yz jj zz DFC,MS ‐ MR = jjj 1.11 1.50 − 0.71zzz jj z j 0.74 − 0.71 3.25 zz (51) k { Here, the diagonal elements are reported with respect to the S0 energy at the FC point. The off-diagonal element DIP is 0.15 eV larger by MS-MR-XMS-CASPT2 than by SS-SR-MSCASPT2, and the S0−S1 gap energy increases by 0.28 eV, which is approximately twice that amount. We also optimized the FC geometry using the two-state averaging scheme. In contrast to the PSB3 case, the S0−S1 gap energy increases by 0.19 eV when the two-state averaging scheme is used. Interestingly, the resulting gap energy is not significantly affected by the use of the MS-MR-XMS-CASPT2 theory (increases by 0.04 eV). The vertical absorption energy with the MS-MR-XMS-CASPT2 theory is closer to the experimental value (2.59 eV).76 From the differences between the gap energy values with and without the IPEA shift92 in the previous studies (see Table 5),75,83 there is a possibility that adding the standard IPEA shift (0.25 Eh) may improve these

Figure 8. Optimized MECI geometries of pHBI. (a and b) P-twisted MECIs and (c and d) I-twisted MECIs. The geometries in red (a, c) were optimized by SS-SR-MS-CASPT2, while the geometries in blue (b, d) were optimized by MS-MR-XMS-CASPT2. For comparison, the MECIs optimized by SA-CASSCF are also shown in gray.

plane twisting of the bridge hydrogen atom. The MS-MRXMS-CASPT2 theory locates this hydrogen atom approximately 10° farther out-of-plane than the SS-SR-MS-CASPT2 theory. In terms of this coordinate (ϕ = H−C3−C4−C6 torsion), the SS-SR-MS-CASPT2 geometry (|ϕP‑twisted| = 129.11°, |ϕI‑twisted| = 11.24°) compares more favorably than MS-MR-(X)MS-CASPT2 geometries (|ϕP‑twisted| = 135.87°, |ϕI‑twisted| = 23.45°) with the previous MRCISD geometries (|ϕP‑twisted| = 125.25°, |ϕI‑twisted| = 12.95°), although it should be noted that the active spaces in both this work and the MRCISD study are quite small. The rms distances to the SACASSCF MECI geometries are 0.064 Å (SS-SR) and 0.112 Å (MS-MR) and 0.126 Å (SS-SR) and 0.054 Å (MS-MR) for the P-twisted and I-twisted MECIs, respectively. In other words, 3968


Article

Journal of Chemical Theory and Computation these two conical intersections exhibit opposite behaviors in the selection of the zeroth-order Hamiltonian: at the P-twisted MECI, MS-MR-XMS-CASPT2 dislocates the MECI more, while at the I-twisted MECI, SS-SR-XMS-CASPT2 dislocates the MECI more. To unravel the reason behind these changes to the MECIs, we again computed the potential energies by SS-SR-MSCASPT2 and MS-MR-XMS-CASPT2 on the interpolated coordinates between the MECIs computed by SS-SR-MSCASPT2 and MS-MR-XMS-CASPT2. A significant fraction of this interpolated coordinate is the out-of-plane motion of the bridge hydrogen atom. The resulting PESs are shown in Figure 9. In both cases, the MECIs are relocated by MS-MR-XMS-

Figure 9. PESs on the interpolated coordinates between the SS-SRXMS-CASPT2 and MS-MR-XMS-CASPT2 MECIs at the (a) P- and (b) I-twisted MECIs. The potential energies computed by SS-SR-MSCASPT2 (red solid), MS-MR-MS-CASPT2 (blue dashed), and MSMR-XMS-CASPT2 (blue solid) are shown. The SA-CASSCF PES (gray) is also shown for comparison.

Figure 10. Elements of the diabatic potential energy matrix D on the interpolated coordinates between the SS-SR-MS-CASPT2 and MSMR-XMS-CASPT2 MECIs near the (a) P-twisted and (b) I-twisted MECIs. The potential energies were computed by SS-SR-MSCASPT2 (red), MS-MR-MS-CASPT2 (blue), and SA-CASSCF PES (gray).

CASPT2 to a geometry that has higher energy, and this result is due to the differences in the amount of stabilization of the ground and excited states. The ground state is more stabilized in MS-MR than in SS-SR, and the excited state is more stabilized in SS-SR than in MS-MR near both MECIs. The CASSCF conical intersection occurs between the SS-SR-MSCASPT2 and MS-MR-XMS-CASPT2 MECIs, and this finding means that SS-SR-MS-CASPT2 and MS-MR-XMS-CASPT2 relocate the conical intersection in opposite directions. An additional intriguing picture can be obtained by using the diabatic potential energy matrix. The diabatic potential energy matrix elements for the interpolated coordinate are shown in Figure 10 (i.e., the potential energies shown in Figure 9 are transformed into the elements of D). At the P-twisted MECI, the adiabatic states are |S 0 ⟩ = |P ⟩

|S1⟩ =

1 (|I ⟩ + |B⟩) 2

(53)

|S 2 ⟩ =

1 (|I ⟩ − |B⟩) 2

(54)

and the diagonal elements DPP (the S0 state) are more stabilized with MS-MR-XMS-CASPT2. Additionally, the MSMR-MS-CASPT2 coupling between the B and I states is smaller than the SS-SR-MS-CASPT2 coupling, while the diagonal elements DBB and DII are almost the same. Consequently, the MS-MR-MS-CASPT2 S1 state is less stabilized than that of SS-SR-MS-CASPT2, resulting in the MECI having higher energy, as shown in Figure 9. At the I-twisted MECI, the adiabatic states are |S 0 ⟩ = |I ⟩

(52) 3969

(55) DOI: 10.1021/acs.jctc.9b00067 J. Chem. Theory Comput. 2019, 15, 3960−3973

Article

Journal of Chemical Theory and Computation |S1⟩ =

1 (|P⟩ + |B⟩) 2

(56)

|S 2 ⟩ =

1 (|P⟩ − |B⟩) 2

(57)

The dependence of the differential correlation effect on the zeroth-order Hamiltonian near the I-twisted MECI is similar to that at the P-twisted MECI. The only difference is that MSMR-MS-CASPT2 stabilizes the I state (again, the S0 state) more than the P state in this case. Overall, the selection of the zeroth-order Hamiltonian changes the energies of the conical intersections in the GFP chromophore case. In particular, the conical intersections were relocated to points with higher energy when we used MS-MR-MS-CASPT2, and these intersections became less thermally reachable from the FC point. In other words, MS-MR-MS-CASPT2 could slow the nonadiabatic decays more than SS-SR-MS-CASPT2. We note that MS-MR-XMS-CASPT2 and SS-SR-MS-CASPT2 also yield similar results with an approximate state rotation scheme (see the Supporting Information). Finally, we discuss the topology of the conical intersections. The PES scans near the conical intersections, which are constructed in the same way as in Figure 4, are shown in Figure 11. All the conical intersections optimized with SS-SRMS-CASPT2, MS-MR-XMS-CASPT2, and SA-CASSCF have sloped topology. The unphysical bumps near the CASSCF degeneracy are more severe than in the PSB3 case. Thus, for dynamics simulations, in which smooth PESs are desirable, we recommend using the MS-MR-XMS-CASPT2 theory, which always produces smooth PESs near conical intersections due to its construction.

4. CONCLUSIONS In this work, we compared two common zeroth-order Hamiltonians (SS-SR and MS-MR) in the MS-CASPT2 theory for calculating conical intersections. To this end, we developed the analytical gradient theory for SS-SR-MS-CASPT2. The conical intersections for two popular model biological chromophores, PSB3 and pHBI, were optimized. These intersections involve twisting along the double bonds. The structural and energy differences are due to the changes in the differential correlation effects. The 1Bu state in PSB3 and the P and I states at the P- and I-twisted MECIs in pHBI are additionally stabilized by MS-MR-XMS-CASPT2 because of the state-averaged Fock operator in the zeroth-order Hamiltonian. The SS-SR-MS-CASPT2 yielded nonphysical bumps near the CASSCF degeneracies, as expected in its mathematical formalism and from the previous investigations.24,53,55 Therefore, we highly recommend that the MSMR-XMS-CASPT2 theory be used in dynamics simulations to ensure that the PESs are always smooth near the conical intersections. Nevertheless, the energies and structures of the MECIs and the topology of the PESs near the conical intersections vary significantly with the zeroth-order Hamiltonian: for example, the geometrical parameters of the SS-SRMS-CASPT2 MECIs compared more favorably with the previous MRCISD results64 than did the parameters of the MS-MR-(X)MS-CASPT2 MECIs. Therefore, the results from optimization and dynamics simulations should be analyzed with caution.

Figure 11. PES scans near the MECIs of pHBI. The vectors g and h denote the gradient difference and interstate coupling vectors, respectively. For MS-MR-(X)MS-CASPT2, g and h are interchanged for visual clarity. The vectors g and h are shown in the Supporting Information.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.9b00067. 3970


Article


■

(11) Angeli, C.; Cimiraglia, R.; Evangelisti, S.; Leininger, T.; Malrieu, J.-P. Introduction of n-electron Valence States for Multireference Perturbation Theory. J. Chem. Phys. 2001, 114, 10252− 10264. (12) Angeli, C.; Cimiraglia, R.; Malrieu, J.-P. n-electron Valence State Perturbation Theory: A Spinless Formulation and an Efficient Implementation of the Strongly Contracted and of the Partially Contracted Variants. J. Chem. Phys. 2002, 117, 9138−9153. (13) Khait, Y. G.; Song, J.; Hoffmann, M. R. Explication and Revision of Generalized Van Vleck Perturbation Theory for Molecular Electronic Structure. J. Chem. Phys. 2002, 117, 4133−4145. (14) Liu, W.; Hoffmann, M. R. iCI: Iterative CI toward Full CI. J. Chem. Theory Comput. 2016, 12, 1169−1178. (15) Lei, Y.; Liu, W.; Hoffmann, M. R. Further Development of SDSPT2 for Strongly Correlated Electrons. Mol. Phys. 2017, 115, 2696−2707. (16) Nakano, H. Quasidegenerate perturbation theory with multiconfigurational self-consistent-field reference functions. J. Chem. Phys. 1993, 99, 7983−7992. (17) Nakano, H. MCSCF Reference Quasidegenerate Perturbation Theory with Epstein–Nesbet Partitioning. Chem. Phys. Lett. 1993, 207, 372−378. (18) Finley, J.; Malmqvist, P. Å.; Roos, B. O.; Serrano-Andrés, L. The Multi-State CASPT2Method. Chem. Phys. Lett. 1998, 288, 299− 306. (19) Reynolds, R. D.; Shiozaki, T. Zero-Field Splitting Parameters from Four-Component Relativistic Methods. J. Chem. Theory Comput. 2019, 15, 1560−1571. (20) Suo, B.; Lian, Y.; Zou, W.; Lei, Y. Electronic Structure of OsSi Calculated by MS-NEVPT2 with Inclusion of the Relativistic Effects. J. Phys. Chem. A 2018, 122, 5333−5341. (21) Yanai, T.; Kurashige, Y.; Saitow, M.; Chalupský, J.; Lindh, R.; Malmqvist, P.-Å. Influence of the Choice of Projection Manifolds in the CASPT2 Implementation. Mol. Phys. 2017, 115, 2077−2085. (22) Celani, P.; Werner, H.-J. Multireference perturbation theory for large restricted and selected active space reference wave functions. J. Chem. Phys. 2000, 112, 5546−5557. (23) Zaitsevskii, A.; Malrieu, J.-P. Multi-Partitioning Quasidegenerate Perturbation Theory. A New Approach to Multireference MøllerPlesset Perturbation Theory. Chem. Phys. Lett. 1995, 233, 597−604. (24) Granovsky, A. A. Extended Multi-Configuration QuasiDegenerate Perturbation Theory: The New Approach to MultiState Multi-Reference Perturbation Theory. J. Chem. Phys. 2011, 134, 214113. (25) Granovsky, A. A. Firefly version 8. http://classic.chem.msu.su/ gran/firefly/index.html, 2016. (26) Shiozaki, T.; Gyòrffy, W.; Celani, P.; Werner, H. J. Communication: Extended Multi-State Complete Active Space Second-Order Perturbation Theory: Energy and Nuclear Gradients. J. Chem. Phys. 2011, 135, 081106. (27) Park, J. W.; Shiozaki, T. Analytical Derivative Coupling for Multistate CASPT2 Theory. J. Chem. Theory Comput. 2017, 13, 2561−2570. (28) Sen, S.; Schapiro, I. A Comprehensive Benchmark of the XMSCASPT2Method for the Photochemistry of a Retinal Chromophore Model. Mol. Phys. 2018, 116, 2571−2582. (29) Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M. Molpro: a General-Purpose Quantum Chemistry Program Package. WIREs Comput. Mol. Sci. 2012, 2, 242−253. (30) Shiozaki, T. BAGEL: Brilliantly Advanced General Electronicstructure Library. WIREs Comput. Mol. Sci. 2018, 8, e1331. (31) Celani, P.; Werner, H.-J. Analytical Energy Gradients for Internally Contracted Second-Order Multireference Perturbation Theory. J. Chem. Phys. 2003, 119, 5044−5057. (32) Gyòrffy, W.; Shiozaki, T.; Knizia, G.; Werner, H.-J. Analytical Energy Gradients for Second-Order Multireference Perturbation Theory using Density Fitting. J. Chem. Phys. 2013, 138, 104104. (33) MacLeod, M. K.; Shiozaki, T. Communication: Automatic Code Generation Enables Nuclear Gradient Computations for Fully

Derivation of the suggested SS-SR-XMS-CASPT2 method and its gradient theory, the benchmark results using the SS-SR-XMS-CASPT2 method, the branching plane vectors (g and h) at the optimized MECIs of pHBI, and the optimized molecular geometries and resulting energies (PDF)

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Jae Woo Park: 0000-0002-4701-8801 Funding

This work was supported by a research grant from the Chungbuk National University in 2018 (Grant 2018100984) and the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (Grant 2019R1C1C1003657). Notes

The author declares no competing financial interest.

■

ACKNOWLEDGMENTS We thank Dr. Saemee Song (Northwestern) for fruitful scientific discussions and continuous support. We also thank Prof. Young Kee Kang (CBNU) for providing us with computational resources, and we thank Prof. Toru Shiozaki (Northwestern), Dr. Rachael Al-Saadon (Northwestern), Prof. Young Min Rhee (KAIST), and Dr. Chang Woo Kim (KAIST) for stimulating discussions.

■

REFERENCES

(1) Lischka, H.; Nachtigallová, D.; Aquino, A. J. A.; Szalay, P. G.; Plasser, F.; Machado, F. B. C.; Barbatti, M. Multireference Approaches for Excited States of Molecules. Chem. Rev. 2018, 118, 7293−7361. (2) Szalay, P. G.; Muller, T.; Gidofalvi, G.; Lischka, H.; Shepard, R. Multiconfiguration Self-Consistent Field and Multireference Configuration Interaction Methods and Applications. Chem. Rev. 2012, 112, 108−181. (3) Lyakh, D. I.; Musiał, M.; Lotrich, V. F.; Bartlett, R. J. Multireference Nature of Chemistry: The Coupled-Cluster View. Chem. Rev. 2012, 112, 182−243. (4) Hirao, K. Multireference Møller-Plesset Method. Chem. Phys. Lett. 1992, 190, 374−380. (5) Hirao, K. Multireference Møller-Plesset Perturbation Theory for High-Spin Open-Shell Systems. Chem. Phys. Lett. 1992, 196, 397− 403. (6) Hirao, K. State-Specific Multireference Møller-Plesset Perturbation Treatment for Singlet and Triplet Excited States, Ionized States and Electron Attached States of H2O. Chem. Phys. Lett. 1993, 201, 59−66. (7) Andersson, K.; Malmqvist, P.-Å.; Roos, B. O. Second-Order Perturbation Theory with a Complete Active Space Self-Consistent Field Reference Function. J. Chem. Phys. 1992, 96, 1218−1226. (8) Andersson, K.; Malmqvist, P.-Å.; Roos, B. O.; Sadlej, A. J.; Wolinski, K. Second-Order Perturbation Theory with a CASSCF Reference Function. J. Phys. Chem. 1990, 94, 5483−5488. (9) Kozlowski, P. M.; Davidson, E. R. Considerations in Constructing a Multireference Second-Order Perturbation Theory. J. Chem. Phys. 1994, 100, 3672−3682. (10) Angeli, C.; Borini, S.; Cestari, M.; Cimiraglia, R. A Quasidegenerate Formulation of the Second Order n-Electron Valence State Perturbation Theory Approach. J. Chem. Phys. 2004, 121, 4043−9. 3971


Article

Journal of Chemical Theory and Computation Internally Contracted Multireference Theory. J. Chem. Phys. 2015, 142, 051103. (34) Vlaisavljevich, B.; Shiozaki, T. Nuclear Energy Gradients for Internally Contracted Complete Active Space Second-Order Perturbation Theory: Multistate Extensions. J. Chem. Theory Comput. 2016, 12, 3781−3787. (35) Park, J. W.; Shiozaki, T. On-the-Fly CASPT2 Surface-Hopping Dynamics. J. Chem. Theory Comput. 2017, 13, 3676−3683. (36) Aquilante, F.; Autschbach, J.; Carlson, R. K.; Chibotaru, L. F.; Delcey, M. G.; De Vico, L.; Fdez Galván, I.; Ferré, N.; Frutos, L. M.; Gagliardi, L.; Garavelli, M.; Giussani, A.; Hoyer, C. E.; Li Manni, G.; Lischka, H.; Ma, D.; Malmqvist, P. Å.; Muller, T.; Nenov, A.; Olivucci, M.; Pedersen, T. B.; Peng, D.; Plasser, F.; Pritchard, B.; Reiher, M.; Rivalta, I.; Schapiro, I.; Segarra-Martí, J.; Stenrup, M.; Truhlar, D. G.; Ungur, L.; Valentini, A.; Vancoillie, S.; Veryazov, V.; Vysotskiy, V. P.; Weingart, O.; Zapata, F.; Lindh, R. Molcas 8: New capabilities for multiconfigurational quantum chemical calculations across the periodic table. J. Comput. Chem. 2016, 37, 506−541. (37) Yanai, T.; Saitow, M.; Xiong, X. G.; Chalupsky, J.; Kurashige, Y.; Guo, S.; Sharma, S. Multistate Complete-Active-Space SecondOrder Perturbation Theory Based on Density Matrix Renormalization Group Reference States. J. Chem. Theory Comput. 2017, 13, 4829− 4840. (38) Mai, S.; Marquetand, P.; González, L. Intersystem Crossing Pathways in the Noncanonical Nucleobase 2-Thiouracil: A TimeDependent Picture. J. Phys. Chem. Lett. 2016, 7, 1978−1983. (39) Richter, M.; Mai, S.; Marquetand, P.; González, L. Ultrafast Intersystem Crossing Dynamics in Uracil Unravelled by ab initio Molecular Dynamics. Phys. Chem. Chem. Phys. 2014, 16, 24423− 24436. (40) Liu, L.; Liu, J.; Martínez, T. J. Dynamical Correlation Effects on Photoisomerization: Ab Initio Multiple Spawning Dynamics with MSCASPT2 for a Model trans-Protonated Schiff Base. J. Phys. Chem. B 2016, 120, 1940−1949. (41) Mai, S.; Marquetand, P.; González, L. Nonadiabatic Dynamics: The SHARC Approach. WIREs Comput. Mol. Sci. 2018, 8, e1370. (42) Tao, H.; Levine, B. G.; Martínez, T. J. Ab Initio Multiple Spawning Dynamics Using Multi-State Second-Order Perturbation Theory. J. Phys. Chem. A 2009, 113, 13656−13662. (43) Coe, J. D.; Levine, B. G.; Martínez, T. J. Ab Initio Molecular Dynamics of Excited-State Intramolecular Proton Transfer Using Multireference Perturbation Theory. J. Phys. Chem. A 2007, 111, 11302−11310. (44) Manathunga, M.; Yang, X.; Luk, H. L.; Gozem, S.; Frutos, L. M.; Valentini, A.; Ferre, N.; Olivucci, M. Probing the Photodynamics of Rhodopsins with Reduced Retinal Chromophores. J. Chem. Theory Comput. 2016, 12, 839−850. (45) Park, J. W.; Shiozaki, T. On the Accuracy of Retinal Protonated Schiff Base Models. Mol. Phys. 2018, 116, 2583−2590. (46) Roos, B. O.; Andersson, K. Multiconfigurational Perturbation Theory With Level Shift - the Cr2 Potential Revisited. Chem. Phys. Lett. 1995, 245, 215−223. (47) Bearpark, M. J.; Robb, M. A.; Schlegel, H. B. A Direct Method for the Location of the Lowest Energy Point on a Potential Surface Crossing. Chem. Phys. Lett. 1994, 223, 269−274. (48) Gozem, S.; Luk, H. L.; Schapiro, I.; Olivucci, M. Theory and Simulation of the Ultrafast Double-Bond Isomerization of Biological Chromophores. Chem. Rev. 2017, 117, 13502−13565. (49) Mori, T.; Nakano, K.; Kato, S. Conical Intersections of Free Energy Surfaces in Solution: Effect of Electron Correlation on a Protonated Schiff Base in Methanol Solution. J. Chem. Phys. 2010, 133, 064107. (50) Fantacci, S.; Migani, A.; Olivucci, M. CASPT2//CASSCF and TDDFT//CASSCF Mapping of the Excited State Isomerization Path of a Minimal Model of the Retinal Chromophore. J. Phys. Chem. A 2004, 108, 1208−1213. (51) Ben-Nun, M.; Molnar, F.; Schulten, K.; Martínez, T. J. The Role of Intersection Topography in Bond Selectivity of cis-trans

Photoisomerization. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 1769− 1773. (52) Zhou, P.; Liu, J.; Han, K.; He, G. The Photoisomerization of 11-cis-Retinal Protonated Schiff Base in Gas Phase: Insight from SpinFlip Density Functional Theory. J. Comput. Chem. 2014, 35, 109−120. (53) Gozem, S.; Huntress, M.; Schapiro, I.; Lindh, R.; Granovsky, A. A.; Angeli, C.; Olivucci, M. Dynamic Electron Correlation Effects on the Ground State Potential Energy Surface of a Retinal Chromophore Model. J. Chem. Theory Comput. 2012, 8, 4069−4080. (54) Gozem, S.; Schapiro, I.; Ferré, N.; Olivucci, M. The Molecular Mechanism of Thermal Noise in Rod Photoreceptors. Science 2012, 337, 1225−1228. (55) Gozem, S.; Melaccio, F.; Valentini, A.; Filatov, M.; HuixRotllant, M.; Ferre, N.; Frutos, L. M.; Angeli, C.; Krylov, A. I.; Granovsky, A. A.; Lindh, R.; Olivucci, M. Shape of Multireference, Equation-of-Motion Coupled-Cluster, and Density Functional Theory Potential Energy Surfaces at a Conical Intersection. J. Chem. Theory Comput. 2014, 10, 3074−3084. (56) Ruckenbauer, M.; Barbatti, M.; Muller, T.; Lischka, H. Nonadiabatic Photodynamics of a Retinal Model in Polar and Nonpolar Environment. J. Phys. Chem. A 2013, 117, 2790−2799. (57) Page, C. S.; Olivucci, M. Ground and Excited State CASPT2 Geometry Optimizations of Small Organic Molecules. J. Comput. Chem. 2003, 24, 298−309. (58) Keal, T. W.; Wanko, M.; Thiel, W. Assessment of Semiempirical Methods for the Photoisomerisation of a Protonated Schiff Base. Theor. Chem. Acc. 2009, 123, 145−156. (59) Manathunga, M.; Yang, X.; Orozco-Gonzalez, Y.; Olivucci, M. Impact of Electronic State Mixing on the Photoisomerization Time Scale of the Retinal Chromophore. J. Phys. Chem. Lett. 2017, 8, 5222− 5227. (60) Manathunga, M.; Yang, X.; Olivucci, M. Electronic State Mixing Controls the Photoreactivity of a Rhodopsin with all-trans Chromophore Analogues. J. Phys. Chem. Lett. 2018, 9, 6350−6355. (61) Gozem, S.; Krylov, A. I.; Olivucci, M. Conical Intersection and Potential Energy Surface Features of a Model Retinal Chromophore: Comparison of EOM-CC and Multireference Methods. J. Chem. Theory Comput. 2013, 9, 284−292. (62) Huix-Rotllant, M.; Filatov, M.; Gozem, S.; Schapiro, I.; Olivucci, M.; Ferré, N. Assessment of Density Functional Theory for Describing the Correlation Effects on the Ground and Excited State Potential Energy Surfaces of a Retinal Chromophore Model. J. Chem. Theory Comput. 2013, 9, 3917−3932. (63) Xu, X.; Gozem, S.; Olivucci, M.; Truhlar, D. G. Combined SelfConsistent-Field and Spin-Flip Tamm−Dancoff Density Functional Approach to Potential Energy Surfaces for Photochemistry. J. Phys. Chem. Lett. 2013, 4, 253−258. (64) Nikiforov, A.; Gamez, J. A.; Thiel, W.; Huix-Rotllant, M.; Filatov, M. Assessment of Approximate Computational Methods for Conical Intersections and Branching Plane Vectors in Organic Molecules. J. Chem. Phys. 2014, 141, 124122. (65) Atchity, G. J.; Xantheas, S. S.; Ruedenberg, K. Potential Energy Surfaces near Intersections. J. Chem. Phys. 1991, 95, 1862−1876. (66) Tsien, R. Y. The Green Fluorescent Protein. Annu. Rev. Biochem. 1998, 67, 509−544. (67) van Thor, J. J. Photoreactions and dynamics of the green fluorescent protein. Chem. Soc. Rev. 2009, 38, 2935−2950. (68) Meech, S. R. Excited state reactions in fluorescent proteins. Chem. Soc. Rev. 2009, 38, 2922−2934. (69) Polyakov, I. V.; Grigorenko, B. L.; Epifanovsky, E.; Krylov, A. I.; Nemukhin, A. V. Potential Energy Landscape of the Electronic States of the GFP Chromophore in Different Protonation Forms: Electronic Transition Energies and Conical Intersections. J. Chem. Theory Comput. 2010, 6, 2377−2387. (70) Bravaya, K. B.; Grigorenko, B. L.; Nemukhin, A. V.; Krylov, A. I. Quantum Chemistry Behind Bioimaging: Insights from Ab Initio Studies of Fluorescent Proteins and Their Chromophores. Acc. Chem. Res. 2012, 45, 265−275. 3972


Article

Journal of Chemical Theory and Computation (71) Mandal, D.; Tahara, T.; Meech, S. R. Excited-State Dynamics in the Green Fluorescent Protein Chromophore. J. Phys. Chem. B 2004, 108, 1102−1108. (72) Jonasson, G.; Teuler, J. M.; Vallverdu, G.; Mérola, F.; Ridard, J.; Lévy, B.; Demachy, I. Excited State Dynamics of the Green Fluorescent Protein on the Nanosecond Time Scale. J. Chem. Theory Comput. 2011, 7, 1990−1997. (73) Virshup, A. M.; Punwong, C.; Pogorelov, T. V.; Lindquist, B. A.; Ko, C.; Martínez, T. J. Photodynamics in Complex Environments: ab initio Multiple Spawning Quantum Mechanical/Molecular Mechanical Dynamics. J. Phys. Chem. B 2009, 113, 3280−3291. (74) Martínez, T. J. Insights for Light-Driven Molecular Devices from Ab Initio Multiple Spawning Excited-State Dynamics of Organic and Biological Chromophores. Acc. Chem. Res. 2006, 39, 119−126. (75) Martin, M. E.; Negri, F.; Olivucci, M. Origin, Nature, and Fate of the Fluorescent State of the Green Fluorescent Protein Chromophore at the CASPT2//CASSCF Resolution. J. Am. Chem. Soc. 2004, 126, 5452−5464. (76) Nielsen, S. B.; Lapierre, A.; Anderson, J. U.; Pedersen, U. V.; Tomita, S.; Andersen, L. H. Absorption Spectrum of the Green Fluorescent Protein Chromophore Anion In Vacuo. Phys. Rev. Lett. 2001, 87, 228102. (77) Olsen, S.; McKenzie, R. H. A Diabatic Three-State Representation of Photoisomerization in the Green Fluorescent Protein Chromophore. J. Chem. Phys. 2009, 130, 184302. (78) Olsen, S.; Lamothe, K.; Martínez, T. J. Protonic Gating of Excited-State Twisting and Charge Localization in GFP Chromophores: A Mechanistic Hypothesis for Reversible Photoswitching. J. Am. Chem. Soc. 2010, 132, 1192−1193. (79) Olsen, S. A Modified Resonance-Theoretic Framework for Structure-Property Relationships in a Halochromic Oxonol Dye. J. Chem. Theory Comput. 2010, 6, 1089−1103. (80) Olsen, S.; Smith, S. C. Bond Selection in the Photoisomerization Reaction of Anionic Green Fluorescent Protein and Kindling Fluorescent Protein Chromophore Models. J. Am. Chem. Soc. 2008, 130, 8677−8689. (81) Toniolo, A.; Olsen, S.; Manohar, L.; Martínez, T. J. Conical intersection dynamics in solution: The chromophore of Green Fluorescent Protein. Faraday Discuss. 2004, 127, 149−163. (82) Levine, B. G.; Coe, J. D.; Martínez, T. J. Optimizing Conical Intersections without Derivative Coupling Vectors: Application to Multistate Multireference Second-Order Perturbation Theory (MSCASPT2). J. Phys. Chem. B 2008, 112, 405−413. (83) Filippi, C.; Zaccheddu, M.; Buda, F. Absorption Spectrum of the Green Fluorescent Protein Chromophore: A Difficult Case for ab Initio Methods? J. Chem. Theory Comput. 2009, 5, 2074−2087. (84) Zhao, L.; Zhou, P. W.; Li, B.; Gao, A. H.; Han, K. L. Nonadiabatic Dynamics of Isolated Green Fluorescent Protein Chromophore Anion. J. Chem. Phys. 2014, 141, 235101. (85) Bravaya, K. B.; Bochenkova, A. V.; Granovskii, A. A.; Nemukhin, A. V. Modeling of the Structure and Electronic Spectra of Green Fluorescent Protein Chromophore. Russ. J. Phys. Chem. B 2008, 2, 671−675. (86) Conyard, J.; Heisler, I. A.; Chan, Y.; Bulman Page, P. C.; Meech, S. R.; Blancafort, L. A New Twist in the Photophysics of the GFP Chromophore: A Volume-Conserving Molecular Torsion Couple. Chem. Sci. 2018, 9, 1803−1812. (87) Hohenstein, E. G.; Bouduban, M. E. F.; Song, C.; Luehr, N.; Ufimtsev, I. S.; Martínez, T. J. Analytic First Derivatives of Floating Occupation Molecular Orbital-Complete Active Space Configuration Interaction on Graphical Processing Units. J. Chem. Phys. 2015, 143, 014111. (88) Mori, T.; Martínez, T. J. Exploring the Conical Intersection Seam: The Seam Space Nudged Elastic Band Method. J. Chem. Theory Comput. 2013, 9, 1155−1163. (89) Fdez Galván, I.; Delcey, M. G.; Pedersen, T. B.; Aquilante, F.; Lindh, R. Analytical State-Average Complete-Active-Space SelfConsistent Field Nonadiabatic Coupling Vectors: Implementation

with Density-Fitted Two-Electron Integrals and Application to Conical Intersections. J. Chem. Theory Comput. 2016, 12, 3636−3653. (90) Park, J. W.; Rhee, Y. M. Diabatic Population Matrix Formalism for Performing Molecular Mechanics Style Simulations with Multiple Electronic States. J. Chem. Theory Comput. 2014, 10, 5238−5253. (91) Pipek, J.; Mezey, P. G. A Fast Intrinsic Localization Procedure Applicable for ab initio and Semiempirical Linear Combination of Atomic Orbital Wave Functions. J. Chem. Phys. 1989, 90, 4916−4926. (92) Ghigo, G.; Roos, B. O.; Malmqvist, P.-Å. A Modified Definition of the Zeroth-Order Hamiltonian in Multiconfigurational Perturbation Theory (CASPT2). Chem. Phys. Lett. 2004, 396, 142−149. (93) Zobel, J. P.; Nogueira, J. J.; González, L. The IPEA dilemma in CASPT2. Chem. Sci. 2017, 8, 1482−1499. (94) Knizia, G. Intrinsic Atomic Orbitals: An Unbiased Bridge between Quantum Theory and Chemical Concepts. J. Chem. Theory Comput. 2013, 9, 4834−4843. (95) Knizia, G.; Klein, J. E. M. N. Electron Flow in Reaction MechanismsRevealed from First Principles. Angew. Chem., Int. Ed. 2015, 54, 5518−5522.

3973


Single-State Single-Reference and Multistate Multireference Zeroth

Recommend Documents