High-Multiplicity Natural Orbitals in Multireference Configuration

May 28, 2013 - We have previously proposed the use of natural orbitals obtained from a single reference correlated calculation describing a high-multi...
4 downloads 0 Views 2MB Size
Article pubs.acs.org/JPCA

High-Multiplicity Natural Orbitals in Multireference Configuration Interaction for Excited State Potential Energy Surfaces Zhen Lu and Spiridoula Matsika* Department of Chemistry, Temple University, Philadelphia, Pennsylvania 19122, United States ABSTRACT: We have previously proposed the use of natural orbitals obtained from a single reference correlated calculation describing a high-multiplicity state (HMNOs), as a substitute for complete active space self-consistent field (CASSCF) molecular orbitals, in the calculation of singlet excited states using the multi-reference configuration interaction (MRCI) method. In our previous work, we found that the MRCI/HMNO method reliably reproduces vertical excitation energies that can be obtained from the standard MRCI/CASSCF procedure for several representative molecules with an average absolute difference of ca. 0.1 eV. In this paper, we test how the MRCI/HMNO method performs when used to generate excited state potential energy surfaces of three sample systems: ethylene, cyclohexadiene, and uracil. We find that the MRCI/HMNO method gives potential energy surfaces for valence excited states that are in good agreement with those generated with the analogous MRCI/CASSCF method, although there are some cases where special care is needed to avoid discrepancies between the two approaches.



INTRODUCTION In order to understand the photophysics and photochemistry of many molecular systems, the ability to probe the excited state potential energy surface and any conical intersections that lie therein is crucial.1,2 Single reference methods can provide very accurate results on vertical excitations and energies close to the Franck−Condon region; however, they have several limitations. For example, most cannot describe two-electron excitations or bond breaking, and fail at describing conical intersection regions between ground and excited states. One way to remedy these shortcomings is to attempt to augment various single reference methods for specific applications. For example, there are attempts to improve the ability of time-dependent density functional theory (TDDFT) to describe doubly excited states,3 or to describe conical intersections with spin-flip techniques.4−6 Single reference perturbation theory approaches have also been extended to account for near degeneracies and conical intersections.7,8 While these augmentations have been implemented with success, the most reliable approaches rely on multireference methods.9 Multireference methods, however, have a steep scaling with the size of the system, restricting their use to small systems. A common way to reduce the size of the problem is to systematically remove configurations from the full interaction expansion by truncating the orbital manifold. Core electrons can generally be excluded, as their contribution to the total correlation is negligible. One can also achieve large computational savings by truncating the virtual space; however, truncating the molecular orbitals with the highest energy will lead to unpredictable changes in the correlation energy. In order to remedy this problem, one can use natural orbitals instead of canonical molecular orbitals. Natural orbitals are the eigenvectors of the one electron reduced spin density matrix. The most useful property of © 2013 American Chemical Society

natural orbitals that are derived from a correlated calculation is that the occupation number associated with a given orbital is directly related to that orbital’s relative contribution to the correlation energy. This property has been widely exploited ever since natural orbitals were introduced by Löwdin in 1955.10,11 The popular frozen natural orbital (FNO) method was introduced by Barr and Davidson in 1970 in their study of the Ne atom.12 In 1976, Shavitt demonstrated quantitatively that natural orbitals lead to the fastest configuration interaction convergence out of several different orbital transformations.13 Since then, natural orbitals have been used extensively in configuration interaction expansions.14−19 There have also been other applications of natural orbitals of note. In several cases natural orbital methods have been proposed as less expensive alternatives to complete active space self-consistent field (CASSCF).20−22 Natural orbitals have also been demonstrated to help guide the selection of the active space for CASSCF calculations.23,24 Although the topic of natural orbitals is mature, there is a noticeable dearth of work applying natural orbitals to excited states, the few examples being Krylov’s use of FNOs to describe ionized excited states with equation of motion-coupled cluster and the work of Aquilante et al., who used FNO-CASPT2 to describe a few excited states of a glycine dipeptide.25 There has also been work done in the past demonstrating that smooth potential energy surfaces can be produced using natural orbitals,21,22 but these studies limit themselves to ground state potential energy surfaces. Besides Special Issue: Joel M. Bowman Festschrift Received: February 8, 2013 Revised: May 10, 2013 Published: May 28, 2013 7421

dx.doi.org/10.1021/jp401444c | J. Phys. Chem. A 2013, 117, 7421−7430

The Journal of Physical Chemistry A

Article

the natural orbitals, other approaches have also been applied in order to circumvent the CASSCF step.26−30 In a recent work, we have described using natural orbitals obtained from a single reference correlated calculation describing a high-multiplicity state, high multiplicity natural orbitals (HMNOs), with the multireference configuration interaction (MRCI/HMNO) method in order to reduce the scaling and complexity of the standard MRCI/CASSCF method.31 The motivation for using a set of higher multiplicity orbitals comes from the fact that usually the first triplet state in a molecule is described by similar orbitals as the first singlet HOMO−LUMO excited state transition, but simply in a different multiplicity. If more orbitals are required for the description of higher lying valence states, a quintet can also be used with the same basic premise. Although the triplet states tend to be more diradical in nature as compared to the more zwitterionic singlet states,32 the effect may be suppressed by the addition of dynamical correlation at the MRCI step. The fundamental goal behind the MRCI/HMNO method is to generate cheap natural orbitals that still afford respectable energies when used for excited state MRCI. This allows us to replace an expensive multistate averaged CASSCF (SACASSCF) step with a simple, single state calculation that does not require solving the multiconfigurational self-consistent field equations, and simultaneously allows us to truncate a sizable portion of the virtual orbitals at the MRCI step. We have previously found that the MRCI/HMNO method can yield energies that mirror MRCI/CASSCF vertical excitation energies of valence states for a variety of organic molecules with a fraction of the computational effort. However, in order for the HMNO approach to be of greater value, its application must extend to making smooth and accurate slices of excited state potential energy surfaces, as they are required for studying photochemistry. At the full CI limit, the choice of using natural orbitals and molecular orbitals should be inconsequential, as one set is simply a rotated version of the other. However, when we are restricted to small active spaces and are furthermore truncating many high lying virtual orbitals, this is no longer the case. Therefore we must test the method, especially when using orbitals that are highly distorted compared to those of the equilibrium geometry.

found in our previous work.31 Note that the virtual space here is defined by the set of unoccupied orbitals in the orbitalgenerating CISD calculation. In our previous work, the HMNO approach was only applied to obtain vertical excitation energies in isolated geometries. In order to be able to produce smooth PESs, some extra steps are needed. When generating natural orbitals, a Hückel guess is made only for the starting geometry of each system. Subsequent calculations along the pathway use the orbitals from the previous generation step instead of a new Hückel guess. This is similar to what is often done in CASSCF calculations and encourages convergence to the same high multiplicity state along the PES. Since the choice of the high spin reference is very important in determining the outcome of the calculation, we make sure that the high spin reference state stays as consistent as possible in terms of orbital character throughout the calculation, and for this reason a previous guess is used. Complications that may arise are discussed in the section Sensitivity to the Reference State. The molecular orbitals that are obtained from SA-CASSCF will generally cause the MRCI expansion to converge faster as these orbitals are a compromise between all the states that are computed. The natural orbitals are generally best suited for the calculation of only the HOMO−LUMO state, and therefore it usually takes more iterations to converge all states of interest. This is shown in Table 1, where the ratio of iterations using

THEORETICAL DETAILS The potential energy surfaces (PESs) of three systems, ethylene, cyclohexadiene, and uracil, are used to check the HMNOs in this work. The HMNO approach consists of a generation and a production step. During the generation step, natural orbitals are generated from a high spin reference (triplet or quintet). In most cases in this study the natural orbitals are obtained using restricted open-shell Hartree−Fock (ROHF) followed by configuration interaction with single and double excitations (CISD). Natural orbitals from unrestricted Møller− Plesset perturbation theory (UMP2) have also been used in some instances for comparisons. Unless noted explicitly, in this work HMNO implies ROHF/CISD natural orbitals. In the production step, the natural orbitals are truncated and used in a singlet MRCI calculation. The occupation percentage preserved criterion is used for the truncation. This criterion involves virtual space truncation by selecting a percentage of occupation that is to be preserved and then truncating all of the orbitals above this threshold. It was first used by Landau et al.,33 and it is the same as the criterion used in our previous publication. More detailed information on the truncation scheme can be

HMNO vs CASSCF orbitals ranges from 2.0 for ethylene to 1.2 for uracil. However, since the HMNOs are intended to be used with a truncated MRCI, the actual wall time for the MRCI to converge will be less for MRCI/HMNO when truncation is included. For example, when using a 99.0% virtual space for ethylene, the MRCI step only takes 18 s as compared to the 123 s it takes for the standard MRCI to converge. This makes the MRCI calculation using natural orbitals more efficient even if more iterations are needed. We have found that the best choice for the multiplicity used to compute the natural orbitals is the lowest multiplicity that includes the orbitals of interest. Therefore, for a molecule like ethylene, we use a triplet that encompasses the π and π* orbitals, while for a molecule like uracil, we must add the lone pair orbital in addition to the π orbitals, as both nπ* and ππ* states are of interest, and thus we use a quintet. In order to check the effect of using a higher multiplicity state, we have attempted to use natural orbitals from a septet state for uracil, and the results are very similar to the ones obtained from the quintet state. The energies obtained using the quintet vs septet NOs for the S1 state are 5.17 vs 5.13 eV, respectively, and for

Table 1. Number of Iterations Required for MRCI Convergence Using Both CASSCF MOs and HMNOs for the Ground State Geometries of the Three Systems Studied in This Worka MRCI iterations ethylene cyclohexadiene uracil

HMNO

CASSCF

61 66 73

31 43 61

a

The MRCI expansion used for ethylene is CAS(2,2)+doubles, the expansion used for CHD is CAS(6,6)+singles, and the expansion used for uracil is CAS(14,10)+singles. The cc-pVTZ basis set was used for ethylene, and cc-pVDZ was used for the other two molecules.



7422

dx.doi.org/10.1021/jp401444c | J. Phys. Chem. A 2013, 117, 7421−7430

The Journal of Physical Chemistry A

Article

excitations from both doubly occupied and CAS orbitals. The reaction path that was considered is a path between the S0 minimum on cyclohexadiene, the S1/S0 conical intersection, and the S0 minimum on hexatriene. A (6,6) CAS was used for both CASSCF and MRCI calculations that included the six π orbitals in hexatriene and four π orbitals in cyclohexadiene as well as the two σ orbitals that become π during the ringopening. HMNOs were generated using a [π,π,π*,π*] quintet reference, including the four π orbitals in CHD. Two MRCI expansions were compared using both MRCI/HMNO and MRCI/CASSCF. The first MRCI expansion is a second order CI allowing for single and double excitations from the CAS only (denoted CAS(6,6)+doubles), while the second allows for single excitations from both the CAS and the doubly occupied orbitals while freezing six core orbitals (denoted CAS(6,6)+allsingles). The CAS(6,6)+allsingles was also used in combination with truncating virtual orbitals. Finally, the nonradiative decay mechanism for the uracil molecule is studied. A joined LLM path connecting the S0 minimum, S2/S1 conical intersection, and S1/S0 conical intersection was created. These geometries were taken from our previous work.34 The quintet state used for uracil has four unpaired electrons in a lone pair orbital and three π orbitals, [n,π,π*,π*]. The active space in the calculations was either a (12,9) or a (14,10). The orbitals included in the (12,9) active space are the 8 π orbitals as well as a lone pair on an oxygen. This active space is in some cases augmented using an additional σ orbital to generate a (14,10) CAS. This addition was necessary in the HMNO calculations using CISD NOs as one of the π orbitals was heavily mixed with the σ orbital. For consistency, some MRCI/CASSCF calculations were also done in a similar (14,10) active space. Two MRCI expansions were used: the first allows for single excitations from the CAS orbitals only (denoted CAS(12,9)+singles or CAS(14,10)+singles) and is equivalent to a first order CI in the conventional notation. The second allows for single excitations from both the CAS and doubly occupied orbitals with eight core orbitals frozen (denoted CAS(12,9)+allsingles or CAS(14,10)+allsingles). A potential energy surface was also constructed using UMP2 HMNOs. For these calculations, the same quintet reference was used in a UMP2 calculation to generate natural orbitals. The MRCI expansion used for the UMP2 PES was a CAS(12,9)+singles. Unless otherwise specified, all calculations were done using the cc-pVDZ basis set. The effect of basis set on the MRCI/ HMNO method has been examined in our previous publication.31 The COLUMBUS computational package was used for all MRCI calculations.35−38 All MP2 optimizations were performed using GAMESS.39,40 Molecular orbitals were viewed using MOLDEN41 and MacMolPlt42 for COLUMBUS and GAMESS calculations, respectively.

the S2 state are 6.28 vs 6.25 eV, respectively. Generating the NOs for the septet state can take considerably longer because the septet state requires many more CSFs to be described. For example, running a septet CISD calculation for uracil with the cc-pVDZ basis set requires 31,790,104 CSFs while running the quintet requires 20,183,226. Furthermore, we have found that the higher the multiplicity of the state that one wishes to use, the harder it is to converge to the desired state for generating the NOs. For these reasons, it does not seem beneficial to go beyond the quintet state. The choice of orbitals in the active space that will be used in the MRCI calculations, beyond the two (four) orbitals singly occupied in the triplet (quintet) wave function, is also very important. This choice is made using chemical intuition similarly to our choice in CASSCF calculations. So, in some cases reordering of the HMNOs is needed to ensure that the active space includes the desired orbitals. In the PES studies, the orbitals are manually reordered once during the initial setup, and then automatic scripts ensure continuity in the active space along the PES. Again, this resembles what is done in CASSCF calculations of PESs, where checking that the active space has the correct orbitals is essential in photochemical studies. The choice of orbitals in the active space in each system studied is given below. Caution should be given to the usage of the term “active space” in this work. Even when we are not carrying out CASSCF calculations, we will have a complete active space in the reference space in the MRCI calculations. Ethylene is the simplest molecule studied. Potential energy surfaces were generated by a fixed rotation around the carbon carbon double bond in two degree increments. The starting ground state minimum geometry was optimized at the MP2/ccpVDZ level. Several MRCI expansions were compared using both the MRCI/HMNO method and the standard MRCI/ CASSCF. HMNOs were generated using a [π,π*] triplet reference. Here the two (or four) orbitals in the brackets define the orbitals with one electron occupation that form the triplet (or quintet) state. These π orbitals were used to construct the (2,2) active space for the MRCI. The notation (m,n) is used here to denote a complete active space with m electrons in n orbitals. The first MRCI expansion allows for single excitations from both the doubly occupied and CAS orbitals, freezing the two core orbitals. This is denoted CAS(2,2)+allsingles. The second expansion allows for single and double excitations from the CAS orbitals, but keeps the other seven occupied orbitals frozen (denoted CAS(2,2)+doubles). This expansion is equivalent to a second order CI in the conventional notation. The final expansion allows for single and double excitations from both the doubly occupied and CAS orbitals, again freezing only the two core orbitals (denoted CAS(2,2)+alldoubles). Vertical excitation energies for the ground state minimum geometry were also obtained using UMP2 and ROMP2 HMNOs. In the vertical excitation energies studies, the CAS(2,2)+doubles expansion was used. The isomerization reaction of cyclohexadiene to produce hexatriene (CHD−HT) is also studied using the MRCI/ HMNO method. The reaction coordinates were created by generating two linear least motion (LLM) paths. This involves linearly interpolating geometry coordinates between important starting and ending geometries using internal coordinates. The S0 minimum geometries for cyclohexadiene and hexatriene were optimized at the MP2/6-31G(d) level while the asymmetric conical intersection geometry was optimized using MRCI/cc-pVDZ with (6,6) CAS and allowing single



RESULTS AND DISCUSSION Ethylene. Ethylene is the smallest prototype system for studying photoinduced cis−trans isomerization about a carbon carbon double bond.43,44 It has a low lying π1π*1 state which, upon bond twisting, leads to a conical intersection. A doubly excited π0π*2 state and several Rydberg states exist above S1. In this work we focus on the two valence states only, as description of Rydberg states is considerably more difficult to describe, and is out of the scope of this work. In our previous paper, we calculated the vertical excitations from the equilibrium geometry as well as the conical intersection energy. 7423

dx.doi.org/10.1021/jp401444c | J. Phys. Chem. A 2013, 117, 7421−7430

The Journal of Physical Chemistry A

Article

CAS(2,2)+alldoubles as described in Theoretical Details. The first excited state is the π1π*1 singly excited state, and the second excited state is either the σ1π*1 or the π0π*2 doubly excited state. The potential energy surfaces generated using the triplet HMNOs are smooth and compare very well with those generated from MRCI/CASSCF. The gradient of the S2 state is discontinuous around 20° because of curve crossing at that geometry. For small angles the S2 state is a σ1π1 state, while for angles between 20° and 160° it is a π0π*2 state. The differences, ΔE = E(HMNO) − E(CASSCF), between the two methods for the same three expansions are shown at the bottom panels of Figure 1. Many interesting observations can be made about these difference plots. First, the differences have a smooth behavior along the twist angle, whose shape and magnitude depend on the state, geometry, and MRCI expansion. The difference is, in most cases, smallest at the S0 minimum and largest at the most distorted 90° region, with a difference of 0−0.4 eV depending on the state and the MRCI expansion. The discontinuity on the S2 surface due to the change of character in the state is also evident in the differences, indicating even more clearly that the character of the state determines how well the HMNO orbitals can reproduce the CASSCF results. ΔE for the S0 and S2 states has anticorrelated behavior, while ΔE for the S1 state shows a smaller variation as a function of geometry. This is presumably because the HOMO and LUMO both change along the path in opposite directions, and this change is reflected more in S0 and S2, which each have two electrons in one of the orbitals and zero in the other. So any change in the HOMO alone or in the LUMO alone will be reinforced. On the other hand, the differences seem to be canceled for S1 as it has one electron in each orbital. The error for CAS+allsingles is smaller than for the other two expansions. Single excitations from all orbitals allow for orbital rotations, and this could compensate for some differences between MOs and NOs. On the other hand, when all orbitals below the active space are frozen in the MRCI expansion, no contribution from them is allowed and any differences between the orbitals in the CAS will be enhanced. Ideally, ΔE should be constant along the PES. In that case the MRCI/CASSCF and MRCI/HMNO PES are parallel to each other, the shape remains the same, and the predicted dynamics will be the same. The S1 surface is the best behaved in this regard, while the fact that the difference is not constant along the PESs of S0 and S2 may create differences in studies of the excited state dynamics when the MRCI/HMNO approach is used compared to MRCI/CASSCF. The reason the S1 state is the best behaved state may be related to our motivation for developing this approach. The shape of the triplet state PES closely resembles the one for the S1 state, as both of them have a π1π*1 configuration. Thus the HMNOs will be best suited for that state. This result seems to indicate that in some cases the HMNOs may be best used for a particular state rather than all the available states. Abrams and Sherrill20 have studied the ground state PES of ethylene along the same torsional coordinate using natural orbitals for the singlet ground state obtained from a variety of methods and used in a CASCI expansion. The main difference with our work is that we use NOs from the triplet state, and we allow for excitations to virtual orbitals. That work highlighted some of the difficulties of obtaining the correct behavior at 90°, where restricted single reference methods will completely fail and produce a cusp at that geometry. For example, using the CISD NOs from the singlet state in the CASCI expansion

Here, we generate a one-dimensional cut of the potential energy surface by starting with the S0 minimum geometry and twisting about the carbon carbon double bond. When the bond has been twisted so that the participating p type orbitals are orthogonal to one another, we expect to see the S0 and S1 energies approach each other as the π bond is broken, creating two sets of isoenergetic p orbitals. The PESs of the S1 and S2 states reach a minimum at 90°, while the ground state PES has a maximum. As we have only selected one nuclear coordinate, the torsion about the double bond, we only see an avoided crossing; in order to fully describe the conical intersection, the pyramidalization of the CH2 group must also be achieved. Initially, we examine the vertical excitation energies using a variety of orbitals in the MRCI CAS(2,2)+doubles expansion. Table 2 shows the vertical excitation energies for the π1π*1 and Table 2. Vertical Excitation Energies for Ethylene Using a Variety of Natural Orbitalsa

a

state

MRCI/UMP2

MRCI/ROMP2

MRCI/CISD

MRCI/CASSCF

π1π*1 π0π*2

9.66 14.43

9.92 15.08

9.66 14.44

9.57 14.34

The MRCI expansion used is CAS(2,2)+doubles.

π0π*2 states in ethylene using a variety of natural orbitals. The last column shows the energies obtained from the MRCI using the CASSCF orbitals, and the best NOs should match those energies. Natural orbitals obtained from CISD and UMP2 perform the best and give almost identical results. The CISD HMNOs will be used for the subsequent study of the PESs. The top panels of Figure 1 show the potential energy surfaces for the three lowest lying singlet states in ethylene generated from both the MRCI/HMNO and MRCI/CASSCF methods. Three different MRCI expansions are shown from left to right: the CAS(2,2)+allsingles, CAS(2,2)+doubles, and

Figure 1. Top panels: Overlaid comparison between the energies produced by MRCI/CASSCF and MRCI/HMNO for the three lowest singlet states in ethylene using an MRCI expansion described as (a) CAS(2,2)+allsingles, (b) CAS(2,2)+doubles, and (c) CAS(2,2)+alldoubles. Bottom panels: Difference between the energies produced by MRCI/CASSCF and MRCI/HMNO (ΔE = E(HMNO) − E(CASSCF)) for the lowest three singlet states of ethylene for the corresponding cases shown in the top panels. 7424

dx.doi.org/10.1021/jp401444c | J. Phys. Chem. A 2013, 117, 7421−7430

The Journal of Physical Chemistry A

Article

MRCI/HMNO, the absolute energies in hartrees are very comparable in addition to the relative energies in eV. A closer look requires looking at the difference between MRCI/CASSCF and MRCI/HMNO, which is plotted in Figure 2c,d. The difference between the HMNO and CASSCF based MRCI results is always less than 0.2 eV, smaller than the differences seen in ethylene. Furthermore, the differences, especially at the second part of the LLM, show small variations as a function of geometry, which means that excited state dynamics will not be greatly affected if the CASSCF orbitals are substituted with the HMNO ones. The errors do show some discontinuous behavior around the conical intersection. This is partly due to the conical intersection, and partly due to the behavior of the underlying quintet state, which will be addressed further below. Uracil. The largest molecule we study in this work is uracil, which has been well studied and benchmarked. Nucleic acid bases are biological chromophores that participate in important photochemical processes. As a set of suitably complex systems, it would be excellent if the MRCI/HMNO method can accurately describe the potential energy surfaces for these types of molecules. Unlike the previous molecules that we have studied in this work, uracil includes more than just ππ* excited states, so it is a good case to test the ability of the HMNO approach to treat multiple characters of states at once. The additional state in uracil is an nπ* state, which, although still a valence state, has quite different electron distribution compared to the ππ* state. Another difficulty in uracil, not present in our previous test cases but present in many organic molecules, is the fact that a proper active space requires many orbitals. For example, in uracil, the active space that should be used includes either 9 or 10 orbitals. These are all the π,π* orbitals, and the two lone pair orbitals on the oxygens. Previous CASSCF studies showed that there is not much difference between the (12,9) or (14,10) active space,34 while smaller active spaces will give poor results for the excited states. This can be problematic, as there are only four orbitals defining the quintet state. We often have to rearrange the orbitals below or above these four orbitals in order to obtain the desired CAS. In uracil it has been found that radiationless relaxation after excitation to the bright ππ* S2 state occurs through multiple conical intersections. Initially, an S2/S1 conical intersection facilitates decay to the dark nπ* S1 state, and finally the S1/S0 conical intersection leads back to the ground state.34 A simple test of the MRCI/HMNO method is to construct LLM paths between these different key geometries to see if the method gives reasonable potential energy surfaces for this complicated system. Figure 3 shows the MRCI results obtained for the two LLM paths which are combined in one. Point 11 is the S2/S1 conical intersection geometry (denoted ci21), and point 22 is the S1/S0 conical intersection geometry (denoted ci10). These were obtained previously at the same MRCI/CASSCF level as the one shown in Figure 3a (CAS(12,9)+singles).34 In the same figure, the HMNO results are shown as well, where the NOs are taken either from CISD or from UMP2. When using the CISD HMNOs the (14,10) active space had to be used in order to have all the π orbitals in the CAS, as one of them was mixed with a σ, as already described in Theoretical Details. The PESs generated using CISD HMNOs are smooth and continuous, as in the previous molecules. However, in this case, the potential energy surfaces generated for uracil with the MRCI/HMNO method do not compare to the MRCI/

produces a cusp at 90°. Our results show that using the triplet NOs from CISD eliminates this problem. Cyclohexadiene to Hexatriene Reaction. The cyclohexadiene to hexatriene photoisomerization (CHD−HT) reaction has been well studied45−48 as a good model system for photochromism, a nondestructive process involving light induced reversible isomerization. The major pathway for the ring-opening reaction occurs through a conical intersection seam. In the simplest case, the path between the reactants and products retains C2 symmetry, but it has been found that lower energy conical intersections exist, which involve a breaking of the symmetry to achieve degeneracy between S0 and S1. The LLM presented here consists of two parts: the first part is between the CHD S0 minimum and the S1/S0 asymmetric conical intersection, while the second part starts at the conical intersection and ends at the HT S0 minimum. The excited states of CHD are similar to those shown here for ethylene, in that the first excited state is a π1π*1 singly excited state and the second excited state is a π0π*2 doubly excited state. The doubly excited state plays a dominant role in the photochemistry, as an avoided crossing occurs between π1π*1 and π0π*2. Figure 2 compares the PESs along the LLM paths obtained using the MRCI/HMNO and MRCI/CASSCF approaches.

Figure 2. Top panels: Overlaid comparison between the energies produced by MRCI/CASSCF and MRCI/HMNO for the three lowest singlet states in the CHD−HT reaction using an MRCI expansion described as (a) CAS(6,6)+allsingles expansion and (b) the CAS(6,6)+doubles expansion. The LLM path is connecting the equilibrium structure of CHD to the ring-opening conical intersection (ci10) with 11 intermediate geometries, and ci10 to the equilibrium structure of hexatriene with 11 more geometries. Bottom panels: Difference between the energies produced by MRCI/CASSCF and MRCI/ HMNO (ΔE = E(HMNO) − E(CASSCF)) for the lowest three singlet states of CHD−HT for the corresponding cases shown in the top panels.

Two MRCI expansions were used here as described in Theoretical Details. Figure 2a shows the CAS(6,6)+allsingles expansion while Figure 2b shows the CAS(6,6)+doubles. The MRCI PESs seen in Figure 2 are continuous and well shaped within each LLM path. As in ethylene, the MRCI/HMNO results seem to be behaving very well. It is furthermore interesting to note that, when compared to the results from 7425

dx.doi.org/10.1021/jp401444c | J. Phys. Chem. A 2013, 117, 7421−7430

The Journal of Physical Chemistry A

Article

occupied and active orbitals. Figure 3b shows the potential energy surface generated using a (14,10) CAS and allowing single excitations from both the σ orbitals and the CAS orbitals. Although there is still a 0.1−0.4 eV energy difference between the two methods between steps 5 and 17, the curve is in much better agreement than in the previous calculation where all the σ orbitals were frozen. There is also a much smaller energy gap in the conical intersection regions, indicating that this approach is able to predict the conical intersections accurately. Allowing for single excitations from the σ orbitals proves to be crucial in the HMNO calculations. We first saw this in ethylene where the difference between MRCI/HMNO and MRCI/CASSCF was reduced in CAS+allsingles. In uracil the improvement is even more dramatic. Without inclusion of single excitations the performance is very poor, while it improves considerably when they are included. Since σ correlation is believed to have a large effect on the excitation energies for uracil to begin with,49 an ideal calculation should include excitations from the σ orbitals. One of the goals of developing the MRCI/HMNO method is to make such large calculations more feasible since it is much cheaper to run a truncated MRCI/HMNO calculation than the corresponding MRCI/CASSCF calculation. Sensitivity to the Reference State. We have so far assumed that the PES of the triplet or quintet state that is being used to generate the HMNOs is smooth along the reaction coordinate. Figure 4 shows the PESs of the triplet state for ethylene and the quintet states for CHD and uracil. The curves for ethylene and uracil are indeed very smooth and well behaved. As was already mentioned when discussing ethylene, the triplet state behaves very similarly to the singlet S1 state, having a minimum at 90°. The curve for CHD, however, is much more complicated, as seen in Figure 4c. The lowest energy quintet state encounters two curve crossings along the reaction coordinate, and, as a result, the curve has a rather strange shape (green curve in Figure 4c). In fact, the results that we get for the quintet change depending on how we choose the initial guess at each point along the path. The approach that we commonly implement in this work is to use the previous solution as a starting guess for the next points, in order to ensure a smooth behavior. This works well in the cases where the high multiplicity state PES is well behaved. When this approach is applied to CHD, however, the results are not very smooth. These are shown in blue in Figure 4c. After the first crossing is encountered at geometry 7, using the previous wave function as a guess leads to a solution that is not the lowest root but rather the state that retains diabatically the character of the previous point. After two geometries, however, the solution converges back to the lowest root, creating a discontinuity in the PES. Even more interestingly, when we generate the path starting from the product HT (red curve in Figure 4c), the PES obtained is very different, and consequently different natural orbitals are produced. The situation becomes clear if one looks at excited quintet states, in addition to the ground quintet state, which are shown in Figure 4d. It is obvious that several states are coupled and interacting, and this is the reason the ground state looks discontinuous. Therefore, proper description requires many states to be included. This is troublesome because we have previously stated that we attempt to keep the quintet state from changing character by using the previous step’s solution as a new guess. This works quite well for most systems, but it is not a great idea for CHD. Interestingly, the discontinuities seen in the quintet state do not impact the excitation energies greatly

Figure 3. Top panels: Overlaid comparison between the energies produced by MRCI/CASSCF, MRCI/HMNO obtained from CISD, and MRCI/HMNO obtained from UMP2, for the three lowest singlet states in uracil using an MRCI expansion described as (a) CAS(12,9)+singles (CAS(14,10)+singles for the MRCI/CISD) and (b) CAS(14,10)+allsingles. The LLM path connects the equilibrium structure of uracil to ci21 with 11 intermediate geometries, and ci21 to ci10 with 11 more geometries. Bottom panels: Difference between the energies produced by MRCI/CASSCF and MRCI/HMNO from CISD (ΔE = E(HMNO) − E(CASSCF)) for the lowest three singlet states of uracil for the corresponding cases shown in the top panels.

CASSCF ones as favorably, and the results depend heavily on the MRCI expansion. Using the CAS(14,10)+singles expansion (shown in Figure 3a) there is a large deviation between the two methods as the molecule is distorted away from the S0 minimum geometry. The first excited nπ* state is the most well behaved, but the ππ* S2 and ground states are both elevated by about 0.5 eV by the sixth geometry, which results in a delayed S2/S1 conical intersection. The most concerning issue is that the error is not systematic; the deviation between the two methods increases as steps are taken away from the Franck−Condon region, causing a qualitative difference in the potential energy surfaces. The deviations from the CASSCF results are shown in panel c. Since the location of the ci21 conical intersection is at different geometries for the MRCI/ HMNO and the MRCI/CASSCF, ΔE will show discontinuities in this area. The HMNO potential energy surface for uracil indicates the existence of a small barrier where the MRCI/ CASSCF results do not show one. This difference can create major differences when the excited state dynamics are studied, as the existence of a barrier will have major effects in lifetimes and efficiency of radiationless decay. The UMP2 orbitals do much worse, especially for the S2 state. The PES is not smooth and the barrier here is very high, and in general the PES is unacceptable compared to the MRCI/CASSCF one. The major reason for these problems is that being able to select the proper active space along the path was not possible. There were σ orbitals mixing with the desired π orbitals, and we could not isolate the desired π orbitals unless we increased the active space considerably. A much better potential energy surface can be generated for uracil if orbital rotations are allowed between the doubly 7426

dx.doi.org/10.1021/jp401444c | J. Phys. Chem. A 2013, 117, 7421−7430

The Journal of Physical Chemistry A

Article

Figure 4. (a) PES of the triplet state of ethylene along the twist coordinate, (b) PES of the quintet state in uracil along the LLM path, (c) PES of the quintet state of CHD along the CHD−HT reaction LLM path, and (d) PES of first 10 quintet states of CHD along the CHD−HT reaction LLM path calculated using CAS(6,6)+doubles with cc-pVDZ basis set. The Reaction Coordinate label along the x axis implies the twist angle in degrees for panel a, and LLM path points in panels b−d.

dot products for ethylene are almost perfectly diagonal, indicating that the MOs and HMNOs are very similar. In CHD, the orbitals involved in the active space (20−25) are also very similar between the MOs and HMNOs. The orbital that shows mixing with the lower σ orbitals is orbital 20, which is a σ orbital. This can be seen in both Figure 5 and Figure 6. Finally, uracil, being the most complicated molecule in our series, shows the most complicated patterns. Figure 6c shows that orbital 23 of the HMNOs is the equivalent of orbital 24 in the MOs, and this is why we had to include both in active space for HMNO, as discussed in the previous section. Furthermore, the dot products are primarily block diagonal, with the σ orbitals mixing with each other and the CAS orbitals mixing with each other, except for orbital 24, which mixes heavily with all the doubly occupied σ orbitals. This is the reason our results of the PES improve substantially when moving from the (14,10)+singles to the (14,10)+allsingles expansion. Allowing for single excitations from all the doubly occupied σ orbitals compensates for this mixing. Although the plots in Figure 6 show dot products for the equilibrium geometries only, we have looked at these dot products for some selected geometries along the PESs as well. In most cases, the mixing is increased when the molecules are distorted. We further briefly discuss the similarity between the wave functions that are generated from the MRCI/HMNO and MRCI/CASSCF calculations. For CHD and ethylene, a dominant CSF describes each state with a coefficient of 0.91−0.94 in both the CASSCF and HMNO based MRCI. In uracil there is more mixing between the CSFs when the HMNOs are used, but, because the orbitals in the active space

for the CHD−HT reaction, but they do cause some minor discontinuities which may pose problems when gradients are desired. The place where the strange quintet behavior is evident is in the discontinuities of the errors seen in Figure 2c,d. Between geometries 13 and 14 the error of S0 and S2 has a discontinuity. This is the same geometry where the quintet state changes character as well. In conclusion, keeping track of the behavior of the reference state is very important. Comparison of the Wave Functions. Although the main focus of this paper is the performance of the MRCI/HMNO approach in producing smooth PESs, a few words about the underlying wave functions are appropriate. Figure 5 shows the

Figure 5. Orbitals generated with quintet CISD and (6,6) CASSCF for the ground state cyclohexadiene geometry.

orbitals included in the active spaces of the MRCI calculations generated by CASSCF and HMNO CISD for the ground state geometry of CHD. Qualitatively, the orbitals look similar, particularly the π orbitals. In order to present a more quantitative comparison, we examine the dot products between the two sets of orbitals. These are shown in a pictorial way in Figure 6 for the equilibrium geometries of the molecules. The 7427

dx.doi.org/10.1021/jp401444c | J. Phys. Chem. A 2013, 117, 7421−7430

The Journal of Physical Chemistry A

Article

Figure 6. Color plots for dot products between all combinations of HMNO and CASSCF MOs (excluding the frozen core and virtual orbitals) for each molecule: (a) ethylene, (b) CHD, and (c) uracil. All NO generating steps use ROHF-CISD with cc-pVDZ basis set.

truncation of the virtual orbitals. We have already addressed how truncation of the virtual space affects the vertical excitation energies of several molecules near the Franck−Condon region. In this section we examine how truncation affects the PESs and their smoothness. We performed truncated calculations on the ethylene molecule using the larger cc-pVTZ basis set. As previously found, using a larger basis set allows for a better truncation scheme. By having a constant occupation preserved as the criterion, the actual number of virtual orbitals that are truncated may vary along the PES. This is not the case for ethylene, but for CHD, for example, that number is between 21 and 22 when preserving 99.0% occupation and between 15 and 16 when preserving 99.9% occupation. The truncated ethylene results can be found in Figure 7. It is apparent that the energy changes very little; the energy curves largely cannot be differentiated from one another in panel a. In panel b it can be seen the average energy difference between results using the 99.0% and 100% virtual space is 0.025 eV and the average difference between the 99.9% and 100% approaches 0. The largest differences occur when the S2 state is a σπ* state. The important part is that the energy difference curve is smooth. Small discontinuities in the difference curve may cause problems when analytic derivatives are implemented. Discontinuities are only seen when there is a state crossing, in which case this is true for both CASSCF and HMNO approaches, and derivative couplings are needed to address this region properly.

rotate with respect to each other, this does not mean that the wave functions are different. A better way to compare the wave functions is to look at the actual density of each state. We compared the densities of the ππ* excited states by using the natural orbitals produced at the MRCI level. For all molecules, there were two dominant NOs with occupation numbers close to one. The dot product between the NOs of the HMNO produced state and the NOs of the CASSCF produced state is 0.99−0.999 in all cases, indicating that the densities of the wave functions are similar. We have also compared other observables derived from the wave function besides the energy. The magnitudes of the dipole moments for the S0 and S1 states in CHD and the S0, S1, and S2 states in uracil obtained from the MRCI/HMNO and the MRCI/CASSCF are shown in Table 3. Their values are very Table 3. Magnitude of Dipole Moments in Debye for Several States at the Ground State Geometry of CHD and Uracila cyclohexadiene S0 S1 ππ* uracil S0 S1 nπ* S2 ππ*

HMNO

CASSCF

0.32318139 0.46164835

0.32078907 0.56408670

4.17508477 1.69554017 5.63310660

4.04010203 1.75089880 5.39089316



a

The MRCI expansion used for CHD is CAS(6,6)+singles and for uracil is CAS(14,10)+singles.

CONCLUSIONS In this work we have explored the efficacy of using our previously developed MRCI/HMNO scheme in obtaining smooth PESs. The results show that the MRCI/HMNO method in most cases predicts smooth PESs with the correct overall behavior. The differences between the MRCI/CASSCF and MRCI/HMNO energies range from 0 to 1 eV, but with proper care the average error is 0.2 eV. The PESs of the ground

similar, with the deviation being between 0.7% and 18%, although most deviations are less than 6%. This is another indication that properties beyond the energy can be calculated accurately with the HMNO approach. Truncation. The main advantage of using the HMNO approach is the ability to save computational effort by 7428

dx.doi.org/10.1021/jp401444c | J. Phys. Chem. A 2013, 117, 7421−7430

The Journal of Physical Chemistry A



Article

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS We gratefully acknowledge support from the National Science Foundation under grant CHE-1213614.



REFERENCES

(1) Domcke, W.; Yarkony, D. R.; Köppel, H. Conical Intersections: Theory, Computation and Experiment; World Scientific: Singapore, 2011. (2) Matsika, S.; Krause, P. Nonadiabatic Events and Conical Intersections. Annu. Rev. Phys. Chem. 2011, 62, 621−643. (3) Maitra, N. T.; Zhang, F.; Cave, R. J.; Burke, K. Double Excitations within Time-Dependent Density Functional Theory Linear Response. J. Chem. Phys. 2004, 120, 5932−5937. (4) Krylov, A. I. Size-Consistent Wave Functions for Bond Breaking: The Equation-Of-Motion Spin-Flip Model. Chem. Phys. Lett. 2001, 338, 375−384. (5) Krylov, A. I. Spin-Flip Equation-Of-Motion Coupled-Cluster Electronic Structure Method for a Description of Excited States, Bond Breaking, Diradicals, and Triradicals. Acc. Chem. Res. 2006, 39, 83−91. (6) Minezawa, N.; Gordon, M. S. Optimizing Conical Intersections by Spin-Flip Density Functional Theory: Application to Ethylene. J. Chem. Phys. 2009, 113, 12749−12753. (7) Head-Gordon, M.; Oumi, M.; Maurice, D. Quasidegenerate Second-Order Perturbation Corrections to Single-Excitation Configuration Interaction. Mol. Phys. 1999, 96, 593−602. (8) Laikov, D.; Matsika, S. Inclusion of Second-Order Correlation Effects for the Ground and Singly-Excited States Suitable for the Study of Conical Intersections: The CIS(2) Model. Chem. Phys. Lett. 2007, 448, 132−137. (9) Szalay, P. G.; Müller, T.; Gidofalvi, G.; Lischka, H.; Shepard, R. Multiconfiguration Self-Consistent Field and Multireference Configuration Interaction Methods and Applications. Chem. Rev. 2012, 112, 108−181. (10) Löwdin, P.-O. Quantum Theory of Many-Particle Systems. I. Physical Interpretations by Means of Density Matrices, Natural SpinOrbitals, and Convergence Problems in the Method of Configurational Interaction. Phys. Rev. 1955, 97, 1474−1489. (11) Löwdin, P.; Shull, H. Natural Orbitals in the Quantum Theory of Two-Electron Systems. Phys. Rev. 1956, 101, 1730−1739. (12) Barr, T. L.; Davidson, E. R. Nature of the ConfigurationInteraction Method in Ab Initio Calculations. I. Ne Ground State. Phys. Rev. A 1970, 1, 644−658. (13) Shavitt, I. I.; Rosenberg, B. J.; Palalikit, S. Comparison of Configuration Interaction Expansions Based on Different Orbital Transformations. Int. J. Quantum Chem. 1976, 10, 33−46. (14) Bender, C. F.; Davidson, E. R. Theoretical Calculation of the Potential Curves of the Be2 Molecule. J. Chem. Phys. 1967, 47, 4972− 4979. (15) Desjardins, S. J.; Bawagan, A. D. O.; Liu, Z. F.; Tan, K. H.; Wang, Y.; Davidson, E. R. Correlation States of Ethylene. J. Chem. Phys. 1995, 102, 6385−6400. (16) Buenker, R. J.; Peyerimhoff, S. D. Individualized Configuration Selection in CI Calculations with Subsequent Energy Extrapolation. Theor. Chim. Acta 1974, 35, 33−58. (17) Thunemann, K. H.; Römelt, J.; Peyerimhoff, S. D.; Buenker, R. J. A Study of the Convergence in Iterative Natural Orbital Procedures. Int. J. Quantum Chem. 1977, 11, 743−752. (18) Sherrill, C. D.; Schaefer, H. F. The Configuration Interaction Method: Advances in Highly Correlated Approaches. In Advances in Quantum Chemistry; Academic Press: 1999; Vol. 34, pp 143−269.

Figure 7. (a) PES of the triplet state of ethylene along the twist coordinate using 100%, 99.9%, and 99.0% occupation and cc-pVTZ basis set with CAS(2,2)+doubles MRCI expansion. (b) Difference in energy between the calculations using truncated and nontruncated virtual space, Δ = E(truncated) − E(nontruncated).

state as well as different types of excited states on three molecules have been examined, and it is found that all of them can be described simultaneously, making the use of high multiplicity natural orbitals a suitable substitute for the more expensive state averaged CASSCF molecular orbitals. The deviation from MRCI/CASSCF and how this deviation changes along the PES depends on the nature of the state. Regions of conical intersections can also be described accurately. An important advantage of using natural orbitals is the ability to truncate the virtual space without loss of accuracy. Our tests indicate that truncating the virtual space has very small effects on the energies of the PESs. Although the MRCI/HMNO method does very well when applied to potential energy surfaces, the underlying behavior of the reference triplet or quintet state should be taken into account when choosing this approach. If the reference state itself exhibits curve crossings and nonadiabatic behavior, this can cause problems in the PESs of the singlet states which are the goal of the study. In systems where many orbitals have to be included in the active space in addition to the two or four orbitals defining the high multiplicity state, one has to be extra careful. In the case of the uracil molecule, the method does not correctly predict the correct qualitative shape of the ππ* excited state if the commonly used (12,9) active space is used, due to the highly mixed nature of the occupied natural orbitals. This becomes problematic as more and more orbitals are required to create a good active space. The effect of this can be mitigated by simply increasing the size of the active space, but this also makes the calculation more expensive. In these cases, using a very large active space while truncating a vast portion of the virtual manifold may help mitigate the added cost of increasing the active space to allow for more orbital rotation. Overall, the MRCI/HMNO approach is promising for studying PESs but great care has to be used when applying it. 7429

dx.doi.org/10.1021/jp401444c | J. Phys. Chem. A 2013, 117, 7421−7430

The Journal of Physical Chemistry A

Article

(19) Ivanic, J.; Ruedenberg, K. Deadwood in Configuration Spaces. II. Singles + Doubles and Singles + Doubles + Triples + Quadruples Spaces. Theor. Chim. Acta 2002, 107, 220−228. (20) Abrams, M. L.; Sherrill, C. D. Natural Orbitals as Substitutes for Optimized Orbitals in Complete Active Space Wavefunctions. Chem. Phys. Lett. 2004, 395, 227−232. (21) Bofill, J. M.; Pulay, P. The Unrestricted Natural OrbitalComplete Active Space (UNO-CAS) Method: An Inexpensive Alternative to the Complete Active Space-Self-Consistent-Field (CAS-SCF) Method. J. Chem. Phys. 1989, 90, 3637−3647. (22) Grev, R. S.; Schaefer, H. F. Natural Orbitals From Single and Double Excitation Configuration Interaction Wave Functions: Their use in Second-Order Configuration Interaction and Wave Functions Incorporating Limited Triple and Quadruple Excitations. J. Chem. Phys. 1992, 96, 6850−6856. (23) Jensen, H. J. A.; Jørgensen, P.; Ågren, H.; Olsen, J. SecondOrder Møller-Plesset Perturbation Theory as a Configuration and Orbital Generator in Multiconfiguration Self-Consistent Field Calculations. J. Chem. Phys. 1988, 88, 3834−3839. (24) Pulay, P.; Hamilton, T. P. UHF Natural Orbitals for Defining and Starting MC-SCF Calculations. J. Chem. Phys. 1988, 88, 4926− 4933. (25) Aquilante, F.; Todorova, T. K.; Gagliardi, L.; Pedersen, T. B.; Roos, B. O. Systematic Truncation of the Virtual Space in Multiconfigurational Perturbation Theory. J. Chem. Phys. 2009, 131, 34113−34119. (26) Potts, D. M.; Taylor, C. M.; Chaudhuri, R. K.; Freed, K. F. The Improved Virtual Orbital-Complete Active Space Configuration Interaction Method, A ‘Packageable’ Efficient Ab Initio Many-Body Method for Describing Electronically Excited States. J. Chem. Phys. 2001, 114, 2592. (27) Chaudhuri, R. K.; Chattopadhyay, S.; Mahapatra, U. S.; Freed, K. F. Molecular Applications of Analytical Gradient Approach for the Improved Virtual Orbital-Complete Active Space Configuration Interaction Method. J. Chem. Phys. 2010, 132, 034105. (28) Slavicek, P.; Martinez, T. J. Ab Initio Floating Occupation Moluecular Orbital-Complete Active Space Configuration Interaction: An Efficient Approximation to CASSCF. J. Chem. Phys. 2010, 132, 234102. (29) Iwata, S. Valence Type Vacant Orbitals for Configuration Interaction Calculations. Chem. Phys. Lett. 1981, 83, 134−138. (30) Okada, K.; Iwata, S. Ab initio MO study of the A, D and third 2 states of CO+. J. Electron. Spectrosc. Relat. Phenom. 2000, 108, 225− 234. (31) Lu, Z.; Matsika, S. High-Multiplicity Natural Orbitals in Multireference Configuration Interaction for Excited States. J. Chem. Theory Comput. 2012, 8, 509−517. (32) Salem, L.; Rowland, C. The Electronic Properties of Diradicals. Angew. Chem., Int. Ed. 1972, 11, 92−111. (33) Landau, A.; Khistyaev, K.; Dolgikh, S.; Krylov, A. I. Frozen Natural Orbitals for Ionized States within Equation-Of-Motion Coupled-Cluster Formalism. J. Chem. Phys. 2010, 132, 014109. (34) Matsika, S. Radiationless Decay of Excited States of Uracil through Conical Intersections. J. Phys. Chem. A 2004, 108, 7584−7590. (35) Lischka, H.; Shepard, R.; Brown, F. B.; Shavitt, I. New Implementation of the Graphical Unitary Group Approach for Multireference Direct Configuration Interaction Calculations. Int. J. Quantum Chem. 1981, 20, 91−100. (36) Shepard, R.; Shavitt, I.; Pitzer, R. M.; Comeau, D. C.; Pepper, M.; Lischka, H.; Szalay, P. G.; Ahlrichs, R.; Brown, F. B.; Zhao, J. A Progress Report on the Status of the COLUMBUS MRCI Program System. Int. J. Quantum Chem. 1988, 34, 149−165. (37) Lischka, H.; Shepard, R.; Pitzer, R. M.; Shavitt, I.; Dallos, M.; Müller, T.; Szalay, P. G.; Seth, M.; Kedziora, G. S.; Yabushita, S.; et al. New High-Level Multireference Methods in the Quantum-Chemistry Program System COLUMBUS: Analytic MR-CISD and MR-AQCC Gradients and MR-AQCC-LRT for Excited States, GUGA Spin-Orbit CI, and Parallel CI Density. Phys. Chem. Chem. Phys. 2001, 3, 664− 673.

(38) Lischka, H.; Shepard, R.; Shavitt, I.; Pitzer, R. M.; Dallos, M.; Müller, T.; Szalay, P. G.; Brown, F. B.; Ahlrichs, R.; Böhm, H. J.; et al. COLUMBUS, An Ab Initio Electronic Structure Program, Release 5.9.1; 2006. (39) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; et al. Computation of Conical Intersections by Using Perturbation Techniques. J. Comput. Chem. 1993, 14, 1347−1363. (40) Gordon, M. S.; Schmidt, M. W. Theory and Applications of Computational Chemistry: The First Forty Years; Elsevier: Amsterdam, 2005; pp 1167−1189. (41) Schaftenaar, G.; Noordik, J. H. Molden: A Pre and PostProcessing Program for Molecular and Electronic Structures. J. Comput.-Aided Mol. Des. 2000, 14, 123−134. (42) Bode, B. M.; Gordon, M. S. Macmolplt: A Graphical User Interface for GAMESS. J. Mol. Graphics Modell. 1998, 16, 133−138. (43) Levine, B. G.; Martnez, T. J. Isomerization Through Conical Intersections. Annu. Rev. Phys. Chem. 2007, 58, 613−634. (44) Barbatti, M.; Paier, J.; Lischka, H. Photochemistry of Ethylene: A Multireference Configuration Interaction Investigation of the Excited-State Energy Surfaces. J. Chem. Phys. 2004, 121, 11614− 11624. (45) Tamura, H.; Nanbu, S.; Nakamura, H.; Ishida, T. A Theoretical Study of Cyclohexadiene/Hexatriene Photochemical Interconversion: Multireference Configuration Interaction Potential Energy Sufraces and Transition Probabilities for the Radiationless Decays. Chem. Phys. Lett. 2005, 401, 487−491. (46) Tamura, H.; Nanbu, S.; Ishida, T.; Nakamura, H. Ab Initio Nonadiabatic Quantum Dynamics of Cyclohexadiene/Hexatriene Ultrafast Photoisomerization. J. Chem. Phys. 2006, 124, 84313. (47) Patel, P. D.; Masunov, A. E. Theoretical Study of Photochromic Compounds: Part 3. Prediction of Thermal Stability. J. Phys. Chem. C. 2011, 115, 10292−10297. (48) Deb, S.; Weber, P. M. The Ultrafast Pathway of Photon-Induced Electrocyclic Ring-Opening Reactions: The Case of 1,3-Cyclohexadiene. Annu. Rev. Phys. Chem. 2011, 62, 19−39. (49) Epifanovsky, E.; Kowalski, K.; Fan, P.-D.; Valiev, M.; Matsika, S.; Krylov, A. I. On the Electronically Excited States of Uracil. J. Phys. Chem. A 2008, 112, 9983−9992.

7430

dx.doi.org/10.1021/jp401444c | J. Phys. Chem. A 2013, 117, 7421−7430