On the accuracy of retinal protonated Schiff base models

On the accuracy of retinal protonated Schiff base models

Jae Woo Park jwpk1201@northwestern.edu    Toru Shiozaki Department of Chemistry, Northwestern University, 2145 Sheridan Rd., Evanston, IL 60208, USA.
July 19, 2019

We investigate the molecular geometries of the ground state and the minimal energy conical intersections (MECIs) between the ground and first excited states of the models for the retinal protonated Schiff base in the gas phase using the extended multistate complete active space second-order perturbation theory (XMS-CASPT2). The biggest model in this work is the rhodopsin chromophore truncated between the and carbon atoms, which consists of 54 atoms and 12-orbital conjugation. The results are compared with those obtained by the state-averaged complete active space self-consistent field (SA-CASSCF). The XMS-CASPT2 results suggest that the minimum energy conical intersection associated with the so-called 13–14 isomerization is thermally inaccessible, which is in contrast to the SA-CASSCF results. The differences between the geometries of the conical intersections computed by SA-CASSCF and XMS-CASPT2 are ascribed to the fact that the charge transfer states are more stabilized by dynamical electron correlation than the diradicaloid states. The impact of the various choices of active spaces, basis sets, and state averaging schemes is also examined.

I Introduction

Photoisomerization is among the non-adiabatic photochemical processes, which undergo electronic transitions near the conical intersections (CIs).Baer (2006) In particular, the photoisomerization of the retinal protonated Schiff base (RPSB) is one of the fastest processes in biology, which occurs within several picoseconds.Schoenlein et al. (1991); Kandori et al. (1995); Ernst et al. (2014); Bassolino et al. (2015); Gozem et al. (2017) The RPSB and related chromophores have been exhaustively studied using theoretical and computational methods. For instance, geometry optimizations, potential energy scans,Cembran et al. (2004); Andruniów et al. (2004); Muñoz-Losa et al. (2011); Ben-Nun et al. (2002); Gozem et al. (2012a); Mori et al. (2010); Valsson and Filippi (2010); Gozem et al. (2012b); Fantacci et al. (2004); Zhou et al. (2014); Garavelli et al. (1998); Gozem et al. (2014); Herbert et al. (2016); Dokukina et al. (2017); Sekharan et al. (2006); Valsson and Filippi (2012); Szymczak et al. (2008); Manathunga et al. (2017); Luk et al. (2015); Schapiro et al. (2011); Manathunga et al. (2016); Gozem et al. (2017) and nonadiabatic molecular dynamics (MD) simulationsFrutos et al. (2007); Liu et al. (2016); Szymczak et al. (2008); Ruckenbauer et al. (2010); Li et al. (2011); Schapiro et al. (2011); Manathunga et al. (2016, 2017); Valentini et al. (2017); Virshup et al. (2009); Luk et al. (2015); Gozem et al. (2017) of truncated model chromophores have been performed in such studies. Because the photoisomerization of RPSB involves transitions between the and electronic states, multireference wave function methods that provide a balanced description of the ground and excited states have often been used to describe the electronic structure of the model chromophores. They include, for instance, state-averaged complete active space self-consistent field (SA-CASSCF),Garavelli et al. (1998); Ben-Nun et al. (2002); Fantacci et al. (2004); Cembran et al. (2004); Andruniów et al. (2004); Frutos et al. (2007); Szymczak et al. (2008); Mori et al. (2010); Valsson and Filippi (2010); Ruckenbauer et al. (2010); Li et al. (2011); Muñoz-Losa et al. (2011); Schapiro et al. (2011); Gozem et al. (2012b, a); Valsson and Filippi (2012); Liu et al. (2016); Manathunga et al. (2016, 2017); Valentini et al. (2017); Luk et al. (2015); Gozem et al. (2017) complete active space second-order perturbation theory (CASPT2),Fantacci et al. (2004); Cembran et al. (2004); Andruniów et al. (2004); Frutos et al. (2007); Szymczak et al. (2008); Mori et al. (2010); Valsson and Filippi (2010); Schapiro et al. (2011); Muñoz-Losa et al. (2011); Valsson and Filippi (2012); Gozem et al. (2012b, a); Zhou et al. (2014); Liu et al. (2016); Manathunga et al. (2016, 2017); Dokukina et al. (2017); Luk et al. (2015); Sekharan et al. (2006); Gozem et al. (2017) and multireference configuration interactions (MRCI).Szymczak et al. (2008); Ruckenbauer et al. (2010); Gozem et al. (2012a, 2017)

The optimization of the molecular geometries of the small RPSB models using CASPT2 was first performed a few decades ago, in which the equilibrium geometries for two protonated Schiff bases ( and ) and the conical intersection of were considered.Page and Olivucci (2003) The ground-state equilibrium structures of the protonated Schiff bases up to (PSB5) have recently been studied using MS-CASPT2.Valsson and Filippi (2010); Dokukina et al. (2017) The nuclear gradients were calculated using finite difference formulas in these studies, and the derivative coupling elements were computed using SA-CASSCF.Page and Olivucci (2003); Valsson and Filippi (2010); Dokukina et al. (2017) The use of numerical gradients, however, does limit the size of the systems to be investigated, because the number of the nuclear degrees of freedom () is multiplied to the computational expense of the underlying quantum chemical methods.

Though several studiesMori et al. (2010); Liu et al. (2016) (using a partially contracted variant, or RS2Celani and Werner (2003); Mori and Kato (2009); Shiozaki et al. (2011)) have demonstrated for small molecules the importance of including dynamical correlation in photodynamics simulations, such simulations for large RPSB models had not been possible until recently, due to the lack of analytical nuclear gradient code for CASPT2 with full internal contraction. Therefore, many of the previous studies have used a hybrid method in which geometries are obtained using CASSCF and energies are computed using CASPT2 (often denoted as CASPT2//CASSCF).Fantacci et al. (2004); Cembran et al. (2004); Andruniów et al. (2004); Frutos et al. (2007); Schapiro et al. (2011); Muñoz-Losa et al. (2011); Gozem et al. (2012b, a); Zhou et al. (2014); Manathunga et al. (2016) For the protonated Schiff bases, the CASPT2 energies at the CASPT2 and CASSCF optimized geometries of the excited states relative to the minimum have shown to differ by 0.11–0.28 eV,Page and Olivucci (2003) and the CASPT2 energy gaps at the conical intersection geometries computed by CASSCF have been around 0.2 eV.Gozem et al. (2012a) Very recently, our group has addressed this problem and reported an efficient parallel implementation of the analytical nuclear gradientsMacLeod and Shiozaki (2015); Vlaisavljevich and Shiozaki (2016); Park and Shiozaki (2017a) and derivative couplingsPark and Shiozaki (2017b) for CASPT2 and its multistate variants. This new development has motivated us to revisit the truncated RPSB models for retinal using CASPT2.

In this work, we have optimized the molecular structures of the RPSB model chromophores using the extended multistate CASPT2 (XMS-CASPT2). The biggest model consists of 54 atoms and conjugation with 12 orbitals. We have shown how the energies change as a function of the conjugation lengths in these model chromophores and have examined the impact of various basis sets, active spaces, and state averaging schemes. The physical origin of such changes is discussed.

Ii Model Chromophores

Figure 1: The structure of rhodopsin and model chromophores used in this study. The main photoisomerization coordinates are highlighted with red (11–12) and blue (13–14) arrows.

Four truncated models of the G protein coupled receptor rhodopsin are studied (see Fig. 1 for their structures and atom labels). The largest model, RPSB6, is a rhodopsin chromophore with a truncated single bond between the and carbons. This RPSB6 model includes all the double bonds in the chromophore, and is used as a reference in this study. The other three models have an additional truncated single bond on the side of cyclohexene ring. These models are designed to have different numbers of double bonds in the carbon chain, while the –N double bond is included in all models. For brevity, we will denote these model chromophores as RPSB in the following discussions, where is the number of double bonds in these chromophores.

These chromophores can be regarded as modifications of the models that have been investigated in previous studies. Compared to the minimal model for the RPSB photoisomerization [PSB3 ()],Page and Olivucci (2003); Szymczak et al. (2008); Mori et al. (2010); Ruckenbauer et al. (2010); Gozem et al. (2012a, 2014); Liu et al. (2016); Zhou et al. (2014); Szymczak et al. (2008); Valsson and Filippi (2010); Fantacci et al. (2004); Gozem et al. (2012b) RPSB3 has two additional methyl groups at the ends of the chromophore and a branched methyl group attached to the atom. RPSB3 corresponds to the model in the recent studies on the QM/MM photodynamicsManathunga et al. (2016) and optomechanical control of photoisomerization.Valentini et al. (2017) RPSB4 without the branch and edge methyl groups corresponds to the model chromophore in Ref. Cembran et al., 2004. RPSB5 is a similar model to that studied in Refs. Ben-Nun et al., 2002 and Herbert et al., 2016, but with a methyl group on the C-terminal. RPSB5 without the branch or edge methyl groups was also considered in some studies.Muñoz-Losa et al. (2011); Garavelli et al. (1998) RPSB6 corresponds to the model chromophore in the previous CASSCF or DFT investigations, although some of them used the model without the methyl group attached to the nitrogen atom.Virshup et al. (2009); Schapiro et al. (2011); Zhou et al. (2014); Gozem et al. (2012b); Valsson and Filippi (2012); Walczak and Andruniów (2015); Manathunga et al. (2017); Ping et al. (2017) RPSB6 is also the model chromophore in several QM/MM studies.Andruniów et al. (2004); Frutos et al. (2007); Luk et al. (2015) RPSB6 model truncated at the edge of the double bond on the cyclohexene ring was also employed in a QM/MM study.Li et al. (2011) The minimum energy structures and pathways in RPSB3, RPSB4, RPSB5 without the edge and branch methyl groups were also optimized using CASSCF and quantum Monte Carlo.Valsson and Filippi (2010)

Figure 2: CASPT2 potential energy landscape of RPSB6 associated with the 11–12 and 13–14 photoisomerization. The energies are reported in eV with respect to the energy at the all-trans geometry. The (12,12) active space consisting of all of the orbitals was used. Three states were included in the calculations, one of which is not shown.

Iii Computational Details

The ground-state equilibrium geometries at the all-trans, 11-cis, and 13-cis geometries, and the MECIs that correspond to the 11–12 and 13–14 photoisomerizations were optimized using the SA-CASSCF and XMS-CASPT2 methods. These MECIs will be denoted as 11–12CI and 13–14CI in the following.

In the SA-CASSCF and XMS-CASPT2 calculations, the states to be included should be specified. At the Franck–Condon (FC) point, the three lowest states of the RPSB are labeled as , , and states, following the notation for polyenes in the symmetry.Gozem et al. (2012a, b); Manathunga et al. (2017); Gozem et al. (2017) The state corresponds to the covalent state at the FC point and the diradicaloid state at the CIs. The state is the charge transfer state, and the state is another diradicaloid state with two or more diradical sites.Gozem et al. (2012a, b); Manathunga et al. (2017); Gozem et al. (2017) The main results were obtained by including the three lowest singlet states in the SA-CASSCF and XMS-CASPT2 calculations, while we also present in Sec. IV.4 the results with a scheme that includes the two lowest singlet states.

Active spaces were the full valence -spaces [(6e,6o), (8e,8o), (10e,10o) and (12e,12o) for RPSB3, RPSB4, RPSB5 and RPSB6, respectively]. Smaller active spaces [(8e,8o) and (10e,10o)] were also tested for RPSB6 in Sec. IV.2. In the XMS-CASPT2 calculations, a vertical shift of was used.Roos and Andersson (1995) The so-called SS-SR contraction scheme was employed.Vlaisavljevich and Shiozaki (2016) The calculations were performed using the cc-pVDZ basis set and its corresponding JKFIT basis set for density fitting unless otherwise specified. The MECIs were optimized using the gradient projection method.Bearpark et al. (1994) All quantum chemistry calculations were performed using the bagel program package.bag (); Shiozaki (2018)

Iv Results and Discussions

Figure 2 shows the XMS-CASPT2 potential energy landscape of RPSB6 associated with the isomerization dynamics. In the following, we will discuss the accuracy of the truncated models and the computational approaches (active spaces, basis sets, and state averaging schemes) using the labels shown in Fig. 2.

iv.1 Truncated Models

Figure 3: Trends of the energies (in eV) with respect to the chromophore length. The energies at the MECIs with respect to the length of the chromophore compared to (a) and (b) energies at the all-trans FC point.

Figure 3 shows the trends of the energies with respect to the length of the models. The Cartesian coordinates of the optimized geometries are compiled in the supplemental online material. As the conjugation of the chromophore becomes longer, the energies at the 11–12CI and 13–14CI become lower and higher, respectively, when computed using XMS-CASPT2 [Fig. 3(a)]. For RPSB5 and RPSB6, the energy at the 13–14CI is higher than the energy at the all-trans geometry [Fig. 3(b)]. This suggests that the 13–14CI is not thermally reachable from the all-trans geometry. This is consistent with the recent gas phase photoisomerization experiment, in which the 11-cis and 9-cis photoisomers were the major products and the 13-cis photoisomer was the minor product.Coghlan et al. (2015) The trend of the XMS-CASPT2 energies at the 13–14CI is not qualitatively reproduced by SA-CASSCF.

Figure 4: Energies (in eV) of the and states computed by XMS-CASPT2 (full) and SA-CASSCF (dotted) along the interpolated coordinates between the XMS-CASPT2 and SA-CASSCF MECIs in RPSB5.

We scanned the XMS-CASPT2 and SA-CASSCF energies for RPSB5 along the interpolated coordinates connecting the MECIs calculated using XMS-CASPT2 and SA-CASSCF (Fig. 4). The diabatic states that cross at the MECI can be assigned to the charge transfer () and diradicaloid () states. It is found that the state is more stabilized by dynamical correlation than the state at both MECIs; the PT2 contributions to the energies for the state is 0.61 eV (11–12CI) and 1.27 eV (13–14CI) larger than those for the state at the SA-CASSCF MECI geometries of RPSB5. This differential dynamical correlation effect is found to result in large structural changes at the MECIs.

The key structural differences between the MECIs obtained by SA-CASSCF and XMS-CASPT2 are the bond lengths. With the PT2 corrections, the and bonds at the 11–12CI geometry elongate by 0.05 Å, whereas the bond contracts by 0.05 Å. These bond length changes correspond to the so-called bond-length alteration (BLA) coordinate, which is a coordinate defined by the elongation of the double bonds and the contraction of the single bonds in the state (see also Fig. 4).Gozem et al. (2012b, 2017) Displacing the molecular geometry toward the FC point along the BLA coordinate stabilizes the diradicaloid diabatic state, as it alters the single and double bonds in the state.Gozem et al. (2012b, 2017)

The bond of the 13–14CI geometry elongates by 0.10 Å with the PT2 corrections, which corresponds to the BLA coordinate for the 13–14CI.Gozem et al. (2012a) This is in agreement with the previous results for PSB3, where the bond length at the 13–14CI geometry was found to be 0.12 Å longer when computed with RS2 than the geometry with CASSCF.Liu et al. (2016) The torsional angle computed by XMS-CASPT2 is about larger than the angle calculated by SA-CASSCF. This change also leads to additional stabilization of the state relative to the state. The shift of the 13–14CI toward the FC point along the BLA coordinate with dynamical correlation was also demonstrated for PSB3 using a variety of the quantum chemical methods [SI-SA-REKS, EOM-SF-CCSD(dT), MRCISD, CASPT2, XMSQDPT2, QD-NEVPT2].Gozem et al. (2012a, b, 2014, 2017) Note that we also observed the same differential dynamical correlation effects in smaller models (RPSB3 and RPSB4); see the supplemental online material for details.

FC 11–12CI111 The difference between the two numbers corresponds to the energy gap in the CASPT2//CASSCF approach. 13–14CI111 The difference between the two numbers corresponds to the energy gap in the CASPT2//CASSCF approach.
Table 1: Comparison between XMS-CASPT2 and XMS-CASPT2//SA-CASSCF. The vertical excitation energies at the FC point and the MECI energies relative to the respective minimum are compared (in eV).

To quantitatively assess the accuracy of the XMS-CASPT2//SA-CASSCF protocol, we calculated the XMS-CASPT2 energies at the SA-CASSCF optimized geometries for all the models (Table 1). The to vertical excitation energy differences between XMS-CASPT2 and XMS-CASPT2//SA-CASSCF at the FC point are relatively small and are in the range of 0.11–0.28 eV. The CI energy differences between the XMS-CASPT2 and XMS-CASPT2//SA-CASSCF energies are in the range of –0.52 eV, and are significantly larger than the differences at the FC point. For RPSB6, and energies at both MECIs calculated using XMS-CASPT2//SA-CASSCF are below the energies computed using XMS-CASPT2. The XMS-CASPT2 gap energies evaluated at the CASSCF optimized MECI geometries are in the range of 0.2–0.6 eV (11–12CI) and 1.0–1.4 eV (13–14CI).

iv.2 Active Spaces

, 11–12CI , 13–14CI
(8e,8o) 3.71 2.18 () 2.25 ()
(10e,10o) 2.83 2.01 () 1.72 ()
(12e,12o) 3.36 2.58 () 2.20 ()
(8e,8o) 2.40 2.17 () 2.76 ()
(10e,10o) 2.39 1.87 () 2.71 ()
(12e,12o) 2.46 1.89 () 2.86 ()
Table 2: Dependence on the size of active spaces for RPSB6. The energies at the FC point and MECIs are shown in eV relative to the energy at the FC point. The MECI energies relative to the energy at the FC point are in parentheses.

The computational cost of multireference calculations drastically increases with the increasing size of the active space. Therefore, for computational efficiency, it is of practical interest to employ as small active spaces as possible when simulating photodynamics. However, the use of such smaller active spaces has to be carefully benchmarked, since it could severely deteriorate the accuracy. In this section, we have tested two smaller active spaces for RPSB6, (8e,8o) and (10e,10o), and compared with the results obtained from the full valence active space [i.e., (12e,12o)].

Table 2 compiles the energies [, ] relative to obtained by SA-CASSCF and XMS-CASPT2 using these active spaces. The main features in (12e,12o) CASPT2 calculations, i.e., at the 13–14CI being higher than , are reproduced by the calculations with smaller active spaces. Using XMS-CASPT2, the deviations in and from the full -valence results are smaller than 0.10 eV and 0.30 eV, respectively, using both active spaces. To the contrary, the SA-CASSCF results are very sensitive to the selection of the active spaces, with 0.82 eV error in the worst case. This encouraging behavior is not surprising because XMS-CASPT2 partially accounts for electronic correlation outside the active space.

For RPSB6, a calculation of XMS-CASPT2 nuclear gradients and derivative couplings with the (8e,8o) active space is about 3 times faster than that with full valence space. Our results show that the trends in the MECI energies are qualitatively reproduced, implying that a smaller active space may be used in future dynamics simulations using XMS-CASPT2. To use smaller active spaces in dynamics simulations, however, further validations are warranted to ensure the energy conservation throughout the dynamics.

iv.3 Basis Sets

, 11–12CI , 13–14CI
cc-pVDZ 2.93 () 2.56 ()
aug-cc-pVDZ 2.86 () 2.59 ()
cc-pVTZ 2.87 () 2.69 ()
cc-pVQZ 2.84 () 2.74 ()
cc-pVDZ 2.18 () 2.77 ()
aug-cc-pVDZ 2.14 () 2.81 ()
cc-pVTZ 2.18 () 2.93 ()
cc-pVDZ 2.02 () 2.80 ()
aug-cc-pVDZ 2.05 () 2.81 ()
cc-pVTZ 2.12 () 2.89 ()
Table 3: Dependence on the size of basis functions. The energies at the FC point and the MECIs are shown in eV relative to the energy at the FC point. The MECI energies relative to the energy at the FC point are in parentheses.

To investigate how the impact of dynamical correlation changes with the size of basis functions, we have performed molecular geometry optimizations of RPSB3, RPSB4, and RPSB5 with larger basis sets. We used the aug-cc-pVDZ, cc-pVTZ, and cc-pVQZ basis sets for RPSB3 and the aug-cc-pVDZ and cc-pVTZ basis sets for RPSB4 and RPSB5. The results are compiled in Table 3. As the basis-set size becomes larger, the vertical excitation energies computed by XMS-CASPT2 at the all-trans point, , relative to becomes smaller. They decrease by 0.01–0.02 eV (cc-pVTZ), 0.04 eV (cc-pVQZ) and 0.05–0.08 eV (aug-cc-pVDZ) compared to that computed using cc-pVDZ.

The MECI energies at the 13–14CI [relative to ] become higher when larger basis sets are used, though the changes are small (about eV). The trend is ascribed to the fact that the differential dynamical correlation effects are more pronounced with larger basis sets, which can describe more dynamical electron correlation contributions. Based on these systematic results, we expect that the same holds for RPSB6, for which we did not perform the calculations with larger basis sets. It is important to note that our results justify the use of small (double-) basis sets in the previous sections and photodynamics simulations in the future.

iv.4 State Averaging Schemes

11–12CI 13–14CI 11–12CI 13–14CI
RPSB3 3.46 () 2.51 () 2.93 () 2.31 ()
RPSB4 2.66 () 2.41 () 2.26 () 2.44 ()
RPSB5 2.20 () 2.37 () 1.79 () 2.55 ()
RPSB6 2.58 ()111Slightly looser threshold was used due to slow convergence. 2.39 () 1.98 () 2.63 ()
RPSB3 3.98 () 2.75 () 2.93 () 2.56 ()
RPSB4 3.08 () 2.52 () 2.18 () 2.77 ()
RPSB5 2.50 () 2.27 () 2.02 () 2.80 ()
RPSB6 2.58 () 2.20 () 1.89 () 2.86 ()
Table 4: Dependence on the state averaging schemes. The energies at the MECIs are shown in eV relative to the energy at the FC point. The MECI energies relative to the energy at the FC point are in parentheses.

There are two commonly used schemes for state averaging in SA-CASSCF and XMS-CASPT2 when studying the RPSB photodynamics: that including the three lowest singlet statesAndruniów et al. (2004); Gozem et al. (2012a); Garavelli et al. (1998); Muñoz-Losa et al. (2011); Cembran et al. (2004); Ruckenbauer et al. (2010); Manathunga et al. (2017); Liu et al. (2016); Frutos et al. (2007); Dokukina et al. (2017); Luk et al. (2015) and that including the two lowest states.Manathunga et al. (2016); Ben-Nun et al. (2002); Gozem et al. (2012a, b); Valsson and Filippi (2010, 2012); Page and Olivucci (2003); Fantacci et al. (2004); Zhou et al. (2014); Szymczak et al. (2008); Mori et al. (2010); Schapiro et al. (2011) We have performed geometry optimizations of all the models using both schemes. For the sake of brevity, we will abbreviate these schemes as SA3 and SA2 schemes. In addition, we will denote the SA-CASSCF and XMS-CASPT2 calculations with these schemes as SA-CASSCF and SA-XMS-CASPT2, where is 3 or 2.

Unlike the results obtained by the SA3 scheme, the correlation energy difference between the states at the 11–12 MECI of RPSB5 using SA2-CASSCF is below 0.01 eV. Consequently, the 11–12CI geometry from SA2-CASSCF closely resembles the SA2-XMS-CASPT2 geometry [root-mean-square deviations between the SA-CASSCF and XMS-CASPT2 geometries are 0.04 Å (SA2) and 0.23 Å (SA3) for RPSB5]. At the 13–14CI obtained by SA2-CASSCF, the difference between the correlation energy contributions for these two states is about 0.54 eV, which is roughly the half of the difference computed using the SA3 scheme. Though smaller than that obtained by the SA3 scheme, this differential dynamical correlation is sufficient to shift the 13–14CI geometry toward the FC point along the BLA coordinate [root-mean-square deviations are 0.63 Å (SA2) and 0.72 Å (SA3) for RPSB5].

The reason why the SA2 correlation energy difference is smaller than that obtained by the SA3 scheme can be explained as follows. The () and () states have similar electronic characters.Gozem et al. (2012a, 2017); Manathunga et al. (2017) In the XMS-CASPT2 calculations, the state-averaged Fock operator is used to define the zeroth-order Hamiltonian, using which the perturbative correction to the energy is computed. Using the SA2 scheme is equivalent to excluding one state from the state-averaged Fock operator in the SA3 scheme, and therefore, the PT2 correction to and will be larger and smaller, respectively, than the correction obtained based on the SA3 scheme. With SA2, therefore, the differences between the PT2 energy for the and states become smaller than with SA3.

Despite such differences between the SA3 and SA2 schemes, the trends in the MECI energies with respect to the chromophore lengths are similar in XMS-CASPT2, and the MECI energies with both schemes agree within 0.4 eV (Table 4). The MECI geometries with both schemes are also similar when XMS-CASPT2 is used: the root-mean-square deviations between them are 0.11 Å and 0.12 Å for the 11–12CI and 13–14CI of RPSB5, respectively. From this observation, it is unclear to us which scheme is superior to the other in the XMS-CASPT2 calculations. Nevertheless, there will be certain situations that one scheme is preferred over another to obtain more accurate results. For example, a recent molecular dynamics studyManathunga et al. (2017) has shown that the mixing between and becomes more important in the polar environment, especially near the FC point, because the states are more stabilized due to the environment than the state. This implies that it is essential to include in the multistate calculations (i.e., use the SA3 scheme) when solvent or protein environments are considered in the simulation.

V Conclusion

In this study, we have optimized the equilibrium and MECI geometries of the RPSB using many structural and computational models. The largest model, RPSB6, consisted of 54 atoms, for which we used the (12e,12o) active space. We found that, as one increases the model size, the XMS-CASPT2 energies at the MECIs (relative to the all-trans minimum) become lower at the 11–12CI and higher at the 13–14CI, respectively. These results are in sharp contrast with those obtained by SA-CASSCF. The differences between the MECI structures obtained by XMS-CASPT2 and SA-CASSCF are ascribed to the fact that the diradicaloid () and charge transfer () states are stabilized by dynamical correlation to different degrees. As a result, the MECI structures are shifted toward the FC point along the BLA coordinate. Various active spaces, basis sets, and state averaging schemes are tested. It is found that the XMS-CASPT2 results are less sensitive to the selection of active spaces than the SA-CASSCF results. Our results suggest that the basis-set dependence is small (though larger basis sets do increase the differential dynamical correlation effect), justifying the use of double- basis sets when simulating the photodynamics of these chromophores. The differential dynamical correlation effects are also shown to be more pronounced when the S () state is included in the SA-CASSCF and XMS-CASPT2 calculations. These results should serve as the basis for choosing the computational models used in the future studies on RPSB photodynamics.

Vi Acknowledgments

This article is part of the special issue in honor of the 80th birthday of Professor Michael Baer. We are grateful to the Air Force Office of Scientific Research Young Investigator Program (Grant No. FA9550-15-1-0031) for generous financial support. The development of the program infrastructure has been in part supported by National Science Foundation [ACI-1550481 (JWP) and CHE-1351598 (TS)].


  • Baer (2006) M. Baer, Beyond Born-Oppenheimer: Electronic Nonadiabatic Coupling Terms and Conical Intersections (Wiley, New Jersey, 2006).
  • Schoenlein et al. (1991) R. W. Schoenlein, L. A. Peteanu, R. A. Mathies,  and C. V. Shank, Science 254, 412 (1991).
  • Kandori et al. (1995) H. Kandori, Y. Katsuta, M. Ito,  and H. Sasabe, J. Am. Chem. Soc. 117, 2669 (1995).
  • Ernst et al. (2014) O. P. Ernst, D. T. Lodowski, M. Elstner, P. Hegemann, L. S. Brown,  and H. Kandori, Chem. Rev. 114, 126 (2014).
  • Bassolino et al. (2015) G. Bassolino, T. Sovdat, A. S. Duarte, J. M. Lim, C. Schnedermann, M. Liebel, B. Odell, T. D. W. Claridge, S. P. Fletcher,  and P. Kukura, J. Am. Chem. Soc. 137, 12434 (2015).
  • Gozem et al. (2017) S. Gozem, H. L. Luk, I. Schapiro,  and M. Olivucci, Chem. Rev. 117, 13502 (2017).
  • Cembran et al. (2004) A. Cembran, F. Bernardi, M. Olivucci,  and M. Garavelli, J. Am. Chem. Soc. 126, 16018 (2004).
  • Andruniów et al. (2004) T. Andruniów, N. Ferré,  and M. Olivucci, Proc. Natl. Acad. Sci. U.S.A. 101, 17908 (2004).
  • Muñoz-Losa et al. (2011) A. Muñoz-Losa, M. E. Martín, I. F. Galván, M. L. Sánchez,  and M. A. Aguilar, J. Chem. Theory Comput. 7, 4050 (2011).
  • Ben-Nun et al. (2002) M. Ben-Nun, F. Molnar, K. Schulten,  and T. J. Martínez, Proc. Natl. Acad. Sci. U.S.A. 99, 1769 (2002).
  • Gozem et al. (2012a) S. Gozem, M. Huntress, I. Schapiro, R. Lindh, A. A. Granovsky, C. Angeli,  and M. Olivucci, J. Chem. Theory Comput. 8, 4069 (2012a).
  • Mori et al. (2010) T. Mori, K. Nakano,  and S. Kato, J. Chem. Phys 133, 064107 (2010).
  • Valsson and Filippi (2010) O. Valsson and C. Filippi, J. Chem. Theory Comput. 6, 1275 (2010).
  • Gozem et al. (2012b) S. Gozem, I. Schapiro, N. Ferré,  and M. Olivucci, Science 337, 1225 (2012b).
  • Fantacci et al. (2004) S. Fantacci, A. Migani,  and M. Olivucci, J. Phys. Chem. A 108, 1208 (2004).
  • Zhou et al. (2014) P. Zhou, J. Liu, K. Han,  and G. He, J. Comput. Chem. 35, 109 (2014).
  • Garavelli et al. (1998) M. Garavelli, T. Vreven, P. Celani, F. Bernandi, M. A. Robb,  and M. Olivucci, J. Am. Chem. Soc. 120, 1285 (1998).
  • 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).
  • Herbert et al. (2016) J. M. Herbert, X. Zhang, A. F. Morrison,  and J. Liu, Acc. Chem. Res. 49, 931 (2016).
  • Dokukina et al. (2017) I. Dokukina, C. M. Marian,  and O. Weingart, Photochem. Photobiol. 93, 1345 (2017).
  • Sekharan et al. (2006) S. Sekharan, O. Weingart,  and V. Buss, Biophys. J. 91, L07 (2006).
  • Valsson and Filippi (2012) O. Valsson and C. Filippi, J. Phys. Chem. Lett. 3, 908 (2012).
  • Szymczak et al. (2008) J. J. Szymczak, M. Barbatti,  and H. Lischka, J. Chem. Theory Comput. 4, 1189 (2008).
  • Manathunga et al. (2017) M. Manathunga, X. Yang, Y. Orozco-Gonzalez,  and M. Olivucci, J. Phys. Chem. Lett. 8, 5222 (2017).
  • Luk et al. (2015) H. L. Luk, F. Melaccio, S. Rinaldi, S. Gozem,  and M. Olivucci, Proc. Natl. Acad. Sci. U.S.A. 112, 15297 (2015).
  • Schapiro et al. (2011) I. Schapiro, M. N. Ryazantsev, L. M. Frutos, N. Ferré, R. Lindh,  and M. Olivucci, J. Am. Chem. Soc. 133, 3354 (2011).
  • Manathunga et al. (2016) M. Manathunga, X. Yang, H. L. Luk, S. Gozem, L. M. Frutos, A. Valentini, N. Ferré,  and M. Olivucci, J. Chem. Theory Comput. 12, 839 (2016).
  • Frutos et al. (2007) L. M. Frutos, T. Andruniów, F. Santoro, N. Ferré,  and M. Olivucci, Proc. Natl. Acad. Sci. U.S.A. 104, 7764 (2007).
  • Liu et al. (2016) L. Liu, J. Liu,  and T. J. Martínez, J. Phys. Chem. B 120, 1940 (2016).
  • Ruckenbauer et al. (2010) M. Ruckenbauer, M. Barbatti, T. Müller,  and H. Lischka, J. Phys. Chem. A 114, 6757 (2010).
  • Li et al. (2011) X. Li, L. W. Chung,  and K. Morokuma, J. Chem. Theory Comput. 7, 2694 (2011).
  • Valentini et al. (2017) A. Valentini, D. Rivero, F. Zapata, C. García-Iriepa, M. Marazzi, R. Palmeiro, I. F. Galván, D. Sampedro, M. Olivucci,  and L. M. Frutos, Angew. Chem. Int. Ed. 56, 3842 (2017).
  • Virshup et al. (2009) A. M. Virshup, C. Punwong, T. V. Pogorelov, B. A. Lindquist, C. Ko,  and T. J. Martínez, J. Phys. Chem. B 113, 3280 (2009).
  • Page and Olivucci (2003) C. S. Page and M. Olivucci, J. Comput. Chem. 24, 298 (2003).
  • Celani and Werner (2003) P. Celani and H.-J. Werner, J. Chem. Phys. 119, 5044 (2003).
  • Mori and Kato (2009) T. Mori and S. Kato, Chem. Phys. Lett. 476, 97 (2009).
  • Shiozaki et al. (2011) T. Shiozaki, W. Győrffy, P. Celani,  and H.-J. Werner, J. Chem. Phys. 135, 081106 (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).
  • Park and Shiozaki (2017a) J. W. Park and T. Shiozaki, J. Chem. Theory Comput. 13, 3676 (2017a).
  • Park and Shiozaki (2017b) J. W. Park and T. Shiozaki, J. Chem. Theory Comput. 13, 2561 (2017b).
  • Walczak and Andruniów (2015) E. Walczak and T. Andruniów, Phys. Chem. Chem. Phys. 17, 17169 (2015).
  • Ping et al. (2017) Y. Ping, T. Xu, R. Momen, A. Azizi, S. R. Kirk, M. Filatov,  and S. Jenkins, Chem. Phys. Lett 685, 222 (2017).
  • Roos and Andersson (1995) B. O. Roos and K. Andersson, Chem. Phys. Lett. 245, 215 (1995).
  • Bearpark et al. (1994) M. J. Bearpark, M. A. Robb,  and H. B. Schlegel, Chem. Phys. Lett. 223, 269 (1994).
  • (46) bagel, Brilliantly Advanced General Electronic-structure Library. http://www.nubakery.org under the GNU General Public License.
  • Shiozaki (2018) T. Shiozaki, WIREs Comput. Mol. Sci. 8, e1331 (2018).
  • Coghlan et al. (2015) N. J. A. Coghlan, B. D. Adamson, L. Gamon, K. Katani,  and E. J. Bieske, Phys. Chem. Chem. Phys. 17, 22623 (2015).
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