State-specific solvation for restricted active space spin–flip (RAS-SF) wave functions based on the polarizable continuum formalism

The restricted active space spin–flip (RAS-SF) formalism is a particular form of single-reference configuration interaction that can describe some forms of strong correlation at a relatively low cost and which has recently been formulated for the description of charge-transfer excited states. Here, we introduce both equilibrium and nonequilibrium versions of a state-specific solvation correction for vertical transition energies computed using RAS-SF wave functions, based on the framework of a polarizable continuum model (PCM). Ground-state polarization is described using the solvent’s static dielectric constant and in the nonequilibrium solvation approach that polarization is modified upon vertical excitation using the solvent’s optical dielectric constant. Benchmark calculations are reported for well-studied models of photo-induced charge transfer, including naphthalene dimer, C 2 H 4 ⋅ ⋅ ⋅ C 2 F 4 , pentacene dimer, and perylene diimide (PDI) dimer, several of which are important in organic photovoltaic applications. For the PDI dimer, we demonstrate that the charge-transfer character of the excited states is enhanced in the presence of a low-dielectric medium (static dielectric constant ε 0 = 3) as compared to a gas-phase calculation ( ε 0 = 1 ) . This stabilizes mechanistic traps for singlet fission and helps to explain experimental singlet fission rates. We also examine the effects of nonequilibrium solvation on charge-separated states in an intramolecular singlet fission chromophore, where we demonstrate that the energetic ordering of the states changes as a function of solvent polarity. The RAS-SF + PCM methodology that is reported here provides a framework to study charge-separated states in solution and in photovoltaic materials.


I. INTRODUCTION
Multireference wave function methods [1][2][3] are important tools in the quantum chemist's arsenal, for the description of static correlation in systems with stretched bonds or other types of (near-) degeneracies, as well as open-shell systems including radicals and excited states. 2 However, the exponential scaling of computational cost with respect to active space size, for methods such as the complete active-space self-consistent field (CASSCF) theory, limits the size and scope of problems that can be addressed in this way.5][6][7][8][9] Among the SF family of approaches, one of the most important variants is the restricted active space (RAS-)SF method.Configuration interaction (CI) wave functions within the RAS family include all possible configurations within an active space and (at a minimum) single excitations into and out of the active space. 10As such, RAS-SF provides a low-cost, balanced treatment of ground and excited states with multi-radical character. 5,6][7] A recently developed variant of RAS-SF allows analysis of charge-transfer (CT) states, 11 which were previously inaccessible to the theory.The new method can compute not just CT states but also electronic couplings, using a diabatic framework that may facilitate a more in-depth understanding of photophysical processes such as artificial or biological light-harvesting.The position of CT states in the excitation manifold is often sensitive to polarization effects, [12][13][14][15] even in low-dielectric environments such as organic photovoltaic The Journal of Chemical Physics ARTICLE scitation.org/journal/jcpmaterials. 15Therefore, some treatment of the surrounding medium is likely necessary in order to make contact with experiments.
7][18] Specifically, we introduce both equilibrium (staterelaxed) and nonequilibrium versions of the RAS-SF + PCM formalism, in order to account for differential polarization upon electronic excitation of the solute.The nonequilibrium continuum approach 16,[19][20][21][22][23][24] separates the medium's polarization response into fast (electronic) and slow (vibrational and orientational) components, then uses the solvent's optical dielectric constant (ε∞) to describe the fast processes while the static dielectric constant (ε 0 ) describes the total polarization response, including both nuclear and electronic contributions.In this way, nonequilibrium PCMs describe the change in electronic polarization following a sudden change in the charge distribution of the solute, corresponding to vertical electronic excitation [20][21][22] or vertical ionization. 21,23,24In contrast, the state-relaxed algorithm is targeted at equilibrium solvation where all polarization mechanisms are fully relaxed in the excited state, as appropriate for the description of emission.

II. THEORY
The method introduced in this work builds upon the RAS-SF approach for CT states, 11 introducing a PCM framework for solvation effects.The RAS-SF formalism is briefly reviewed in Sec.II A, following which we introduce the state-specific equilibrium and nonequilibrium solvation theories in Secs.II B and II C, respectively.1][22][23][24] Section II D describes how the PCM theory is integrated with the RAS-SF approach to computing the wave function.

A. RAS-SF wave function
A RAS-SF calculation starts from high-spin, restricted openshell Hartree-Fock (ROHF) reference orbitals, from which target states of interest are computed via SF excitations.The RAS scheme divides the orbitals into three subspaces labeled RAS1, RAS2, and RAS3 and the wave function Ansatz for RAS-SF is The RAS2 space is the active space and should include the most important orbitals needed to describe the electronic states of interest; m ∈ RAS2 in Eq. ( 1) includes all configurations within the active space, meaning that RAS-SF resembles a CAS-CI method within the RAS2 space.The RAS1 space contains all doubly occupied orbitals and RAS3 contains virtual orbitals, and we allow single excitations between these two subspaces and RAS2, corresponding to holes h ∈ RAS1 and particles p ∈ RAS3.This particular Ansatz has been called RAS(h,p)-SF. 10,11o describe CT states, 11 the high-spin reference orbitals are first localized onto molecular fragments, where they are then categorized by subspace (RAS1, RAS2, or RAS3) and fragment.The RAS Hamiltonian is then partitioned into CT and non-CT blocks, where the latter will also be called the "locally excited" (LE) block.Diagonalization of these blocks provides either the LE states or the CT states, and the latter are further grouped into left-right vs right-left character, or equivalently A → B vs B → A (Fig. 1).

B. Equilibrium state-specific solvation
Electrostatic interactions between the electronic ground state of the solute and the continuum solvent are described by a reactionfield operator, R0 .The total electronic energy in the ground state is in which Ĥvac is the vacuum (gas-phase) Hamiltonian.The subscript "0" in R0 indicates that the ground-state wave function is used to polarize the medium.In a ground-state self-consistent reaction-field (SCRF) calculation, the energy that is variationally minimized is not actually E 0 but rather the free energy G 0 = E 0 − W 0 , where is the work required to polarize the continuum.A ground-state SCRF calculation therefore corresponds to minimization of a functional G 0 [Ψ] that is given by 16,25 We will now extend this idea to state-specific solvation of an excited-state wave function |Ψ k ⟩, starting with an equilibrium approach in which the polarization of the medium is fully relaxed with respect to the excited-state charge distribution, using the solvent's static dielectric constant (ε 0 ) to represent all possible polarization mechanisms for the medium.This approach should be valid in the long-time limit after excitation, whereas the case of a sudden (Our notation for excited-state PCM calculations follows that used in a recent review. 16) The quantity Rk is the reaction-field operator for state k, meaning that the charge distribution corresponding to wave function |Ψ k ⟩ is used to polarize the medium.Since Rk depends on |Ψ k ⟩, Eq. ( 5) must be solved iteratively and this constitutes the excited-state SCRF problem.
In analogy to Eq. ( 2), the electronic energy for state k in this fully relaxed approach is using a superscript "eq" to indicate equilibrium solvation.The corresponding free energy is where represents the work required to polarize the continuum using Rk .
The fully relaxed (equilibrium) excitation energy is where is the state-specific eigenvalue difference.The quantity W k − W 0 in Eq. ( 9) represents the differential polarization work between the excited state and the ground state.The iterative procedure required to solve the state-specific eigenvalue problem in Eq. ( 5) is relatively straightforward if there are no quasi-degeneracies, although this equation must be solved separately for each electronic state of interest.However, the presence of near-degeneracies may lead to convergence problems associated with root-flipping, and properties other than the energy (including oscillator strengths) are not entirely well-defined because the excited states are not eigenfunctions of a common Hamiltonian and are therefore not orthogonal to one another. 26To circumvent these difficulties, we next describe a perturbative approach to state-specific solvation, [20][21][22] which furthermore allows for all of the solvent-corrected excitation energies to be obtained from a single calculation, including nonequilibrium corrections.

C. Perturbative state-specific solvation
We next describe a nonequilibrium approach to state-specific solvation that is appropriate for modeling vertical excitation energies. 16The nonequilibrium formalism recognizes that nuclear (vibrational and orientational) degrees of freedom within the (implicit) solvent do not respond quickly enough to remain in equilibrium with a vertical electronic transition, and indeed this is the sense in which the excitation is "vertical."There is also an electronic component to the solvent's polarization response, however, and it ought to remain in equilibrium with the solute.The total polarization in the ground state is thus partitioned into "slow" (nuclear) and "fast" (electronic) components, Operationally this means that the ground-state polarization charge σ 0 (generated by R0 ) is partitioned as σ 0 = σ s 0 + σ f 0 , where 16 σ s 0 = ( The solvent's optical dielectric constant (ε∞), which has been called the "dielectric constant for induced polarization," 27 is used to describe electronic polarization.
Considering the equilibrium expression for E k in Eq. ( 6), where Rk = Rs k + Rf k , a corresponding nonequilibrium expression is obtained by instead using Rs 0 + Rf k as the reaction-field operator for excited state k, corresponding to a state-specific polarization charge density σ k = σ s 0 + σ f k .This substitution affords where is a state-specific Hamiltonian for nonequilibrium solvation.
As discussed in the context of the equilibrium case, the statespecific nature of the Hamiltonian may lead to convergence problems and other formal complexities arising from the nonorthogonality of the excited-state wave functions. 16,26As a result of rootflipping problems, it is probably only realistic to imagine a fully self-consistent solution to the state-specific Schrödinger equation in cases where the state k of interest is spectrally isolated, unless specialized convergence algorithms are employed. 26(One such algorithm is described in Ref. 26.)Even if convergence is not problematic, a separate calculation is required for each excited state.For reasons of simplicity, we desire an approach that can generate the entire spectrum in a single calculation.A perturbative approach to incorporating the solvation correction solves both of these problems.
To obtain such an approach, we start from the unperturbed Hamiltonian that is simply a rewriting of the ground-state Hamiltonian in Eq. ( 2), but in a form that suggests a partition of the state-specific Hamiltonian in Eq. ( 14), namely, Introducing a perturbation parameter λ, this partition can be used to develop a perturbation theory expression for the energy of state k, including nonequilibrium corrections: The Journal of Chemical Physics

ARTICLE scitation.org/journal/jcp
1][22] Unperturbed states |Ψ (0) k ⟩ are eigenfunctions of Ĥ0 , and the zeroth-order energy for state k is within the nonequilibrium formalism.The corresponding free energy expression is 16 G where W s 0 and W f 0 are defined analogously to W 0 in Eq. ( 3) but with Rs 0 or Rf 0 replacing R0 .Equation ( 18) for E k is a natural generalization of Eq. ( 2) for E 0 and corresponds to solving the Schrödinger equation in the fixed reaction field of the equilibrated ground state.A nonequilibrium correction arises in first-order perturbation theory and is given by where the notation Rf k(0) indicates that this reaction-field operator is constructed using the charge density corresponding to wave function where in analogy to Eq. ( 8).The quantity W 0,k(0) in Eq. ( 21) is a chargeseparation penalty arising from the Coulomb interaction of the initial-and final-state surface charges, 16,21 Here, φ σ s 0 (r) is the electrostatic potential generated by the groundstate polarization charge, σ s 0 (r).As discussed in Refs.16 and 21, the W 0,k(0) term has sometimes been omitted from nonequilibrium polarization treatments, but is necessary when the "Marcus partition" into fast and slow components is used, corresponding to Eq. (12).
Taking the free energy for state k to be the zeroth-order result plus the first-order correction, the ptSS approximation to the excitation energy is where The formula in Eq. ( 25) has a straightforward interpretation. 16The leading term, which is given by Eq. ( 26), is the difference between eigenvalues E k and E 0 of the state-specific Schrödinger equation, where E k is correct through first order in perturbation theory.This eigenvalue difference becomes a free energy upon subtracting W f k(0) − W f 0 , which is the difference between the work required to polarize the fast charge on the excited state relative to that required on the ground state.Finally, W 0,k(0) is the aforementioned charge-separation penalty.
For future reference, we define which is the zeroth-order approximation to E k − E 0 .Equivalently, this is the excitation energy obtained by solving the Schrödinger equation in the fixed reaction field of the ground state, and therefore this approximation does not contain nonequilibrium effects.Equation ( 25) can then be rewritten as The ptSS approximation for the free energy of excitation (ΔG ) is thus obtained by adding the zeroth-order eigenvalue difference ), with all nonequilibrium effects contained in the latter.

D. Implementation
A flow chart for the RAS-SF + PCM algorithm is illustrated in Fig. 2, including both the nonequilibrium ptSS and the staterelaxed equilibrium procedure.The algorithm begins by computing the RAS-SF ground state wave function in the gas phase.Next, one-electron integrals are modified to incorporate R0 and the ground-state total energy is iterated to convergence, as in any SCRF calculation except that the lowest eigenstate of a CI Hamiltonian must be computed at each SCRF iteration.A damping scheme was implemented for the PCM surface charges, 24 such that the charges at each iteration are a convex linear combination of the charges obtained from the current reaction field operator and the previous one.However, the use of this algorithm proves to be unnecessary.Convergence of the ground-state SCRF problem is typically accomplished in 3-4 iterations if the convergence criterion is set to 10 −8 hartree on the change in G 0 .This completes step 1 in Fig. 2.
Once the ground-state eigenfunction has been equilibrated with respect to the medium, excited states (including LE and/or CT states, as desired) are obtained by computing additional eigenstates in the presence of a fixed reaction-field operator R0 .These are the zeroth-order states From there, the corrected free energies for each state can be computing using Eq. ( 19) or excitation free energies using Eq.(25).Note that ptSS corrections can be computed for all excited states at once, without any additional SCRF iterations in step 2 of Fig. 2.
The equilibrium state-relaxed procedure starts from a gasphase RAS-SF calculation to obtain an initial wave function |Ψ k ⟩ for the excited state of interest, from which Rk can be constructed.This procedure bypasses the ground-state SCRF procedure (although the ground-state SCRF problem needs to be solved in advance, to obtain E 0 ), instead proceeding directly to the state-specific Step 1 of the ptSS approach involves an SCRF calculation using the ground-state eigenvector of the CI Hamiltonian (illustrated with a gray box, as in Fig. 1).In step 2, excited-state eigenvectors are computed from the CI Hamiltonian using a frozen reaction field for the ground state ( R0 ), then corrected using the perturbation Rf k − Rf 0 .In this approach, all of the CI eigenvectors can be corrected at once and only a single CI calculation is required in step 2. Alternatively, the fully relaxed equilibrium procedure corresponds to a SCRF calculation in the presence of reaction-field Rk , solving the state-specific Schrödinger equation [Eq.( 5)] one state at a time to obtain fully relaxed energies E k .
eigenvalue problem for E k [Eq.(5)].The solution of this equation requires SCRF iterations in the excited state because each time the RAS-SF Hamiltonian is diagonalized to obtain |Ψ k ⟩, the operator Rk must be modified to reflect how the new excited-state charge density polarizes the medium.This process is iterated to convergence, which typically requires 7-10 iterations with a convergence criterion of 10 −6 hartree.Upon convergence, E eq k is obtained from Eq. ( 6).The proper state energy, however, is G eq k [Eq.( 7)], and the excitation energy for state k is given in Eq. ( 9).The SCRF procedure must be repeated for every state of interest.
The equilibrium state-relaxed approach has some similarities to the state-specific implementation of time-dependent density functional theory (TD-DFT) with PCM, 28,29 especially if TD-DFT is viewed as a CI method with single excitations (CIS), where excited states are obtained iteratively as roots of an eigenvalue problem.Our procedure also bears some similarities to the most complete version of the SCRF approach for correlated wave functions, which is usually called "perturbation to energy and density," 16,22,[30][31][32][33] meaning that all correlation effects are included in the charge density that is used to equilibrate the reaction field.Some differences exist with respect to the present procedure, as RAS-SF is not a post-Hartree-Fock procedure where correlation is added to a single-determinant reference, but rather it is a method in which the ground-state wave function emerges as the lowest eigenfunction of a CI Hamiltonian.
The ptSS-PCM procedure described in Sec.II C has previously been implemented for CIS and TD-DFT wave functions, 20,21 which is technically simpler because in that case the reference state is the physical ground state.The TD-DFT + ptSS-PCM method 20,21 is very similar to the "corrected linear response theory" of Caricato et al. 34

III. RESULTS AND DISCUSSION
The RAS-SF + PCM procedure has been implemented in a locally modified version of Q-Chem 5.3. 35

A. Computational details
All RAS calculations use a (4,4) active space and the holeparticle (h,p) algorithm that was outlined in Sec.II A, starting from a quintet ROHF reference state.In order to use the CT version of RAS-SF, the ROHF orbitals were localized using the Pipek-Mezey procedure. 36Core orbitals are not correlated (i.e., frozen core approximation).The basis set is 6-31G * , cc-pVDZ, or cc-pVTZ, as indicated in the discussion that follows.Q-Chem's RAS algorithm uses a resolution-of-identity (RI) approximation for the integrals. 5For the auxiliary (density fitting) basis set, we use either the ones designed for use with cc-pVXZ 37 or else those designed for use with the Ahlrichs SVP basis set. 38In Q-Chem, these two basis sets are called RIMP2-cc-pVDZ and RIMP2-VDZ, respectively.
For the solvent model, we use the "conductor-like" (C-)PCM. 16,39The alternative "integral equation formalism" (IEF)-PCM [40][41][42][43] is a formally more complete treatment of dielectric boundary conditions, 16 and therefore in principle might be a better model for low-dielectric environments since C-PCM introduces errors of O(1/ε) in the dielectric boundary conditions. 18,44In practice, however, this difference is usually numerically unimportant and the two models afford similar solvation energies even in nonpolar solvents. 16,43,45,46The more important consideration is that both models contain an implicit correction for outlying charge, 47 i.e., for penetration of the wave function beyond the cavity.From a technical standpoint, IEF-PCM is somewhat more challenging to work with, and in particular its integral operators are more sensitive to discretization as compared to the surface potential operator in C-PCM. 16,48For these reasons, we have opted for the simpler C-PCM approach.
A van der Waals cavity is used to define the interface with the continuum. 16,49That cavity is constructed using Bondi's atomic radii (as modified in Ref. 50), scaled by a factor of 1.2.Atomic spheres were discretized using 194 Lebedev grid points per sphere.The switching/Gaussian (SwiG) method is used for the cavity surface discretization. 25,51or intermolecular CT states, we report some test calculations using constrained density functional theory (cDFT) 52 to move an electron from one monomer to the other using a groundstate formalism.In these calculations, which were performed at the B3LYP/cc-pVTZ level, the charge constraint is implemented using Becke populations with atomic size corrections. 53,54 State-relaxed equilibrium procedure All calculations in this section use the state-specific equilibrium procedure that is described in Sec.II B. The convergence criterion for the ground-state SCRF iterations (used to determine R0 ) is set The Journal of Chemical Physics ARTICLE scitation.org/journal/jcpto 10 −8 hartree and that for the excited-state SCRF iterations (to determine Rk ) is set to 10 −6 hartree.

Distance dependence for CT excitation
][60][61] Both dimers are arranged in a parallel face-to-face orientation.Solvation corrections for singlet CT from C 2 H 4 to C 2 F 4 , and from one naphthalene monomer to the other, are shown in Fig. 3 as a function of reciprocal intermolecular separation, 1/R.The quantity plotted is the difference between excitation energies computed using ε 0 = 78.4vs ε 0 = 1.0, i.e., the solvent shift in water.RAS-SF calculations on C 2 H 4 ⋅ ⋅ ⋅C 2 F 4 use the cc-pVTZ basis set and RIMP2-cc-pVTZ auxiliary basis set while those on (naphthalene) 2 use the 6-31G * basis set and RIMP2-cc-pVDZ basis set.
For a point of comparison, we also computed the CT energies using a ΔSCF approach in conjunction with cDFT, a comparison that was also used in Ref. 11 to test the gas-phase RAS-SF CT procedure.With cDFT, both of the energies required for the ΔSCF calculation can be computed as ground-state DFT + PCM calculations, and comparison to an equilibrium model of the solvent response is therefore appropriate.
For both systems, the behavior of the solvent shift is very similar in both the RAS-SF and the cDFT calculations, despite significant differences in methodology: RAS is a correlated wave function approach that uses excited-state SCRF iterations to compute the state-specific solvation energy, while the cDFT-based ΔSCF procedure uses two ground-state calculations, one of them with a charge constraint that forces the monomers to integrate to ±1 charge.Despite these differences, the solvent shift is found to be a linear function of 1/R, in accordance with the Born model of two wellseparated charges in a dielectric medium. 16For C 2 H 4 ⋅ ⋅ ⋅C 2 F 4 , the shift varies by about 0.6 eV (14 kcal/mol) over a 1.5 Å change in R, at both levels of electronic structure theory, despite the fact that the absolute shifts are offset by about 0.1 eV.The slope (solvent shift vs 1/R) that is obtained from RAS-SF + PCM is within 10% of that obtained using cDFT + PCM for both dimers.These tests suggest that the combination of a PCM with the RAS-SF Ansatz for CT states provides a meaningful description of solvation for a photo-excited ion pair.

Perylene diimide dimer
3][64][65] The overall kinetics of the singlet fission process are solvent-dependent, [65][66][67] and while CT states in PDI materials appear rapidly in polar solvents, high-yield singlet fission occurs equally rapidly in nonpolar solvents. 65,68,69ere, we use RAS-SF + PCM to explore the impact of solvent polarity on the singlet fission process.
Previous work has suggested that CT contributions are mixed into the key multi-exciton (ME) state of interest in singlet fission. 70,71he ME state is nominally a singlet-coupled pair of triplet excitations, 15,72,73 However, it may also contain small contributions from chargeseparated determinants |A + B − ⟩ and |A − B + ⟩.These sometimes appear in the form of charge-resonance (CR) states, 15,59 Despite its charge-separated character, the CR state may have a small or vanishing dipole moment if c 1 ≈ −c 2 , which is often the case for high-symmetry (or even near-symmetry) dimers. 15Gas-phase calculations suggest that the charge-separated states in (PDI) 2 are indeed of CR form. 71ere, we consider two geometries of (PDI) 2 that are taken from Ref. 70 and which afford the fastest and slowest singlet fission rates of the eight (PDI) 2 geometries that were examined in Ref. 71.The geometries of these two dimer models are illustrated on the left side of Fig. 4, and their frontier molecular orbitals are illustrated in the middle part of Fig. 4. The two dimers have somewhat different offset-stacking arrangements for the cofacial PDI monomers and also different face-to-face distances.We will examine the S 1 ( 1 ππ * ) state of these dimers and also two different ME states that will be designated 1 ME (lower state) and 1 ME ′ (upper state), following Ref.71.
Figure 5 reports the percentage of CT character for these three states, as quantified using the RAS-SF CT method, 11 as a function of solvent polarity.We consider the values ε 0 = 1 (corresponding to vacuum boundary conditions), ε 0 = 3 [characteristic of typical thin films, such as poly(methyl methacrylate), PMMA], and ε 0 = 37.5 (representing acetonitrile).For model 1 of the dimer, the CT character varies only slightly as a function of these very significant changes in ε 0 , but for model 2 we observe significant changes in CT character in polar vs nonpolar environments.In particular, the S 1 state of model 2 exceeds 50% CT character in acetonitrile as compared to 30% CT character in the gas phase.
The molecular orbitals that are primarily responsible for the CT character are illustrated in the middle part of Fig. 4 that represent the two dominant configurations in the wave function.The right side of Fig. 4 quantifies the CT contribution in these two configurations, as a function of ε 0 .In the gas phase, the two CT configurations appear in essentially equal percentages, in both FIG. 5. Percentage CT character for the S 1 , 1 ME, and 1 ME ′ states in two different geometries of (PDI) 2 , called models 1 and 2 (or "MO" and "C8") in Fig. 4. The 1 ME ′ state is defined to be a higher-energy state of pair of ME states that emerge from the calculation.
structures of (PDI) 2 , leading in both cases to a CR state with a vanishing dipole moment.For model 2, however, a dielectric constant of ε 0 = 3 is sufficient to break the symmetry, resulting in a net dipole moment.In contrast, symmetry-breaking is not observed in model 1 even in acetonitrile.This observation has potential ramifications for singlet fission.4][75] Model 2, which transfers electrons in only one direction when ε 0 ≳ 3, might act as an energetic trap in a dielectric medium.This trap (an excimer state) would then compete with singlet fission and would lower the yield of the latter process.This suggests that one possible route to improving the efficiency of singlet fission in PDI dimers is to exploit geometries in which singledirectional CT is inhibited.This difference highlights the fact that calculations carried out using vacuum boundary conditions are illpositioned to address whether CT character plays a role in the singlet fission mechanism.

C. State-specific perturbation theory
In this section, we consider the ptSS approach for vertical excitation energies.We first examine the model systems C 2 H 4 ⋅ ⋅ ⋅C 2 F 4 and (naphthalene) 2 that were considered above, but we also investigate some pentacene dimers of interest in organic photovoltaic applications.These include a model system extracted from an intramolecular singlet fission chromophore, 75 and we note that The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp7][78] For these calculations, the convergence criterion on the SCRF iterations was set to 10 −10 hartree, with convergence achieved in 12-16 cycles.

Distance and solvent dependence of CT energies
We first consider the distance dependence of CT excitation energies in C 2 H 4 ⋅ ⋅ ⋅C 2 F 4 and in (naphthalene) 2 , similar to the tests reported in Sec.III B 1 using the equilibrium solvation approach.Comparison to cDFT calculations is not appropriate in the present case, however, because we use a nonequilibrium formulation of the solvation correction in order to describe vertical excitation energies, and there is no analogous ΔSCF calculation.RAS-SF excitation energies in vacuo and in water are shown in Fig. 6 as a function of 1/R.These calculations use the 6-31G * basis set and RIMP2-VDZ auxiliary basis set.
As expected, the excitation energies vary linearly with 1/R and the solvent shifts are quite large, e.g., ΔE = −1.33 eV for C 2 H 4 ⋅ ⋅ ⋅C 2 F 4 at R = 5 Å (RAS-SF/cc-pVDZ level).This is comparable to the shift of ΔE = −1.15eV that was reported previously for this system using TD-DFT + ptSS-PCM. 21able I lists the lowest CT excitation energy that is obtained for C 2 H 4 ⋅ ⋅ ⋅C 2 F 4 in six different solvents including cyclohexane (cHex), toluene, tetrahydrofuran (THF), chlorobenzene (ClBz), benzonitrile (BzCN), and water.Because we consider only the electrostatic part of the solvation energy, each solvent is fully characterized by the pair of values ε 0 and ε∞.Although nonelectrostatic solvation models are reasonably well established for ground-state solvation 16 (e.g., using the SMx models), 79 there has been much less work on nonelectrostatic effects for excited states and we do not consider that topic here.(For a brief overview of this topic, see Ref.  [Eq. ( 25)].c Equation (21).
in nonelectrostatic interactions upon electronic excitation, including dispersion and Pauli repulsion, are expected to be much smaller than changes in the electrostatic solvation energy.
Listed in Table I is the first-order solvation correction to the excitation energy, G ptSS( 1) k [Eq.( 21)], along with the zeroth-order approximation to the excitation energy itself, ΔE ptSS(0) k [Eq.( 27)], and the first-order approximation to the same excitation energy, ΔG ptSS( 1) k [Eq.( 25)].These three quantities are related according to Eq. ( 28).The zeroth-order quantity ΔE ptSS(0) k is the excitation energy computed in the frozen reaction field of the ground state and does not contain nonequilibrium corrections, i.e., it depends on ε 0 but not ε∞.
We note that the nonequilibrium correction G ptSS(1) k exceeds 1 eV for the CT state that is considered in Table I, even in nonpolar solvents.Perhaps counterintuitively, this correction is smaller in water than it is in less polar solvents.This is ultimately a consequence of the partition between "slow" and "fast" (or "inertial" and "noninertial") contributions to the solvent response, 16,80 with the consequence that for a solvent like water where ε 0 ≫ ε∞, a much greater proportion of the solvent response is frozen upon vertical excitation.This can be understood in terms of the "Pekar factor," ε −1 ∞ − ε −1 0 , which replaces the familiar factor of 1 − ε −1 0 in the Born (or generalized Born) solvation energy expression, in cases where vertical excitation energies are involved. 80For example, ε −1 ∞ − ε −1 0 appears in the Marcus theory expression for the outer-sphere reorganization energy, which is the earliest version of a nonequilibrium dielectric continuum theory. 16

Pentacene dimer models
In thinking about CT states that might be relevant in organic photovoltaic applications, 15,72 it is interesting to note that even nonpolar solvents (with ε 0 ≲ 5) afford a large nonequilibrium correction to a CT excitation energy.This can be seen for C 2 H 4 ⋅ ⋅ ⋅C 2 F 4 in Table I, for example.The nonequilibrium correction comes from the optical dielectric constant, which is related to the electronic polarizability of the solvent molecules and has values ε∞ = 1.7-2.3 for most common solvents, 16 even while the static dielectric constant varies from ε 0 = 2-4 for nonpolar hydrocarbons up to ε 0 = 78 for water. of Chemical Physics With this in mind, we next consider an example relevant to organic photovoltaics, namely, a pentacene dimer model of a chromophore that undergoes intramolecular singlet fission, 75 which is depicted in Fig. 7(a).Calculations are reported here using the unlinked dimer model that is shown in Fig. 7(b).We have confirmed, via CIS calculations, that the model system and the full chromophore afford excitation energies within about 0.3 eV of one another, which is sufficient for our purposes.The same is true for two other truncated models (not shown) in which the side chains or the spacer group were removed, but not both, indicating that the tri-isobutyl silyl side chains in the full chromophore [Fig.7(a)] function mainly for solubility while the spacer moiety functions mostly to position the pentacene chromophores relative to one another.The singlet fission dynamics can be quite sensitive to the relative position and orientation of the pentacene moieties, 81 so the atomic positions of the pentacene dimer core are the same in the model system as they are in the full chromophore.Calculations on the model are reported below in six different solvents, and in each case we compute six RAS-SF states: the ground state, the first LE state (which is T 1 ), two left CT states (LCT1 and LCT2), and two right CT states (RCT1 and RCT2).These calculations use the 6-31G * basis set and RIMP2-VDZ auxiliary basis set.Selected results are shown in Table II, focusing on the LCT1 state but also including LE1 for comparison.
We observe that stabilization of the ground state increases with solvent polarity, and at the level of the zeroth-order eigenvalues E (0) k , which are computed in the reaction field of the ground state, approximately the same stabilization is obtained for the LCT1 state as for the ground state itself, to within a few millihartree.The first-order corrected excitation energies show a different trend, however.Here, both the static and optical dielectric constants are in play, and the value of the first-order correction G ptSS(1) k is smaller in water than it is in nonpolar solvents.As explained in Sec.III C 1, this behavior results from the fact that water has a much larger fraction of its total polarization response frozen upon vertical excitation since [Eq.(28)].
The Journal of Chemical Physics ARTICLE scitation.org/journal/jcpε 0 ≫ ε∞.This is a feature of the Marcus partition of the polarization into fast and slow components. 16,21,80In nonpolar or weakly polar solvents, where ε 0 is much closer to ε∞, the nonequilibrium correction amounts to a greater fraction of the total solvation energy of the excited state, and the trend in first-order corrections G ptSS(1) k is governed primarily by the value of ε∞.As noted elsewhere, 16,21,82 the partition between inertial and noninertial solvent response need not be done in this way, and the total solvation correction is essentially the same in Pekar's alternative partition scheme but is apportioned differently into zeroth-order and corrective terms. 16,82For the LE1 state, we also observe a reduction in the excitation energy as solvent polarity increases, but the effect is smaller as compared to the LCT1 state because LE1 does not involve significant charge separation.
We next compare RAS-SF + ptSS-PCM results to the corresponding calculation at the CIS level, which will illustrate the ease of reaching target states with the RAS-SF procedure.Table III shows results for the pentacene dimer model in Fig. 7(b), in two different solvents.In cyclohexane, we observe a correspondence between the RAS-SF states LCT1 (triplet) and LCT2 (singlet) and the 37th and 38th excited states obtained from the CIS calculation.Comparing excitation energies, the RAS-SF and the CIS results agree within 0.04 eV and the dipole moments agree within 0.5 D. A similar correspondence is observed between the RAS-SF states RCT1 and RCT2, and the 39th and 40th excited states obtained from the CIS calculation.The principle pair of natural transition orbitals (NTOs) for states 37 and 39 are plotted in Fig. 8, from the CIS calculation, which confirms the CT character of the states in question.(The RAS-SF states have CT character by construction.)While the excitation energies are similar, indicative of the lack of significant dynamical correlation in the RAS-SF wave function, 7 the CIS approach necessitates the calculation of a large number of excited states, which must then be analyzed in terms of orbitals and amplitudes in order to deduce which states possess CT character.With RAS-SF, on the other hand, one may target the CT states directly (essentially by fiat), so that these emerge as the lowest eigenvalues of the CT block of the Hamiltonian. 11teady state absorption measurements performed on the intramolecular singlet fission chromophore in Fig. 7(a) indicate a red shift in the absorption maxima as solvent polarity increases, which was taken as evidence that the excited state is more polar than the ground state. 75We observe a similar pattern for the unlinked dimer  model, as discussed above, with a decrease in the excitation energy as a function of increasing solvent polarity, although such a shift is observed in both the LE and the CT states.Table IV compares RAS-SF results for the unlinked dimer model to multireference perturbation theory calculations on the full chromophore, from Ref. 75.The latter calculations were performed in vacuo so the same is true for the RAS-SF calculations in Table IV.Both sets of calculations reveal four closely spaced excited states of alternating multiplicity (singlet and triplet), with significant CT character as indicated by excited-state dipole moments μ ≈ 65 D. This comparison suggests at least a qualitative correspondence between the two levels of theory, even if the gas-phase excitation energies are rather different.

TABLE V.
Excitation energies (ΔE) and excited-state dipole moments (μ), computed in benzonitrile solution, for the full chromophore in Fig. 7(a) vs the unlinked dimer model in Fig. 7 In Ref. 75, the effects of solvent polarity were investigated by comparing excitation energies in vacuum to those computed in benzonitrile, using a SCRF model in conjunction with a semiempirical Hamiltonian.These results are reproduced in Table V alongside the corresponding RAS-SF + ptSS-PCM results.Semiempirical calculations in Ref. 75  state hardly shifts at all in benzonitrile solvent, consistent with T 1 results from Ref. 75.The present calculations therefore corroborate the results in Ref. 75, indicating that "polar" states (meaning those with CT character) are considerably shifted whereas "nonpolar" (LE) states are not.Note that the nonpolar S 1 bright state should have essentially the same charge density as T 1 and thus will not be significantly shifted by the solvent model.

IV. CONCLUSIONS
We have formulated, implemented, and tested a state-specific approach to solvation for RAS-SF wave functions.Both equilibrium (self-consistent and fully relaxed) and nonequilibrium (vertical excitation) versions are reported, the latter within a perturbative "ptSS" approximation introduced in earlier work. 20,21This new methodology opens the possibility to describe dielectric continuum effects on CT states within a cost-effective and easy-to-use wave function model for strong correlation, which can describe multi-exciton states and other types of excitations that are beyond the reach of (or poorly described by) single-excitation methods such as CIS and TD-DFT.Benchmark tests indicate good agreement with other methods for well-separated CT excitations (ion-pair states), in systems such as C 2 H 4 ⋅ ⋅ ⋅C 2 F 4 and (naphthalene) 2 .The ptSS approach to vertical excitation energies agrees with literature results for stabilization of CT states vs LE states, in a model of a covalently linked pentacene dimer that undergoes intramolecular singlet fission.
Calculations on (PDI) 2 in two different geometries provide a rationale to understand the very different singlet fission rates that have previously been calculated for these systems.These calculations suggest that dielectric stabilization creates excitonic trap states in one geometry but not the other, by breaking symmetry and localizing charge.This observation has potentially important implications for the design of materials with good singlet fission rates, as the excitation spectra of model chromophores can now be easily tested as a function of solvent polarity.It also suggests that gas-phase calculations, even with highly correlated wave function models, may not offer a realistic description of singlet fission in solution, where polarization-induced charge localization is observed in solvents with dielectric constants as small as ε 0 = 3, but not in vacuum (ε 0 = 1).

SUPPLEMENTARY MATERIAL
See the supplementary material for coordinates for all structures considered.

FIG. 1 .
FIG. 1.Schematic representation of the states that emerge in RAS-SF from a highspin quintet reference configuration.Note that each state (including the ground state) emerges from the diagonalization of a CI Hamiltonian, indicated by the gray box.

FIG. 4 .
FIG.4.RAS-SF/cc-pVDZ + PCM results for two geometries of (PDI) 2 that have been called (a) "MO" and (b) "C8."70,71Structures are shown on the left and correspond to slightly different offsets in the lateral plane as well as different stacking distances, 3.46 Å for model 1 and 3.28 Å for model 2.71 Both models give rise to a pair of 1 ME states, and the middle panel illustrates the orbitals involved in the two largest CT configurations for both states.The tables on the right quantify the percentage CT character in each of these two configurations, along with the corresponding dipole moment, as a function of the solvent's dielectric constant.The state labeled 1 ME ′ is the higher-energy state.

FIG. 7 .
FIG. 7. Structural models for an intramolecular singlet fission chromophore based on pentacene dimer: (a) the original chromophore from Ref. 75, which includes a 1,4-diethynylbicyclo[2.2.2]octane spacer and tri-isobutyl silyl side chains, vs (b) a truncated and unlinked dimer model, in which the side chains are replaced by SiH 3 and the spacer has been replaced by ethynyl side chains.Positions of the atoms that constitute the pentacene dimer core are the same in both systems.

FIG. 8 .
FIG. 8.Principle NTO pair for each of CIS states 37 and 39, depicting a CT excitation between pentacene fragments.
afford four closely spaced singlet excited states (in vacuum), in the range ΔE = 2.63-2.97eV, with dipole moments ranging from μ = 55.0-62.1 D and which are stabilized by 1.4-1.7 eV in benzonitrile whereas other nonpolar states are stabilized by ≤0.01 eV.With RAS-SF and nonequilibrium ptSS-PCM solvation, we obtain four closely spaced states ranging from ΔE = 3.75-3.91eV with dipole moments μ = 64.8-65.0D, which are stabilized by 1.33 eV in benzonitrile.By contrast, the LE1 (T 1 )

The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp (
vertical) excitation requires a nonequilibrium approach, as developed in Sec.II C. For equilibrium solvation on excited state k, one must solve the state-specific Schrödinger equation

TABLE II .
RAS-SF/6-31G * + ptSS-PCM energies for the pentacene dimer model in Fig.7(b), computed in various solvents.Values in parenthesis are shifts relative to the corresponding calculation in vacuum.

TABLE III .
RAS-SF/ and CIS/6-31G * excitation energies and excited-state dipole moments for the pentacene dimer model in Fig.7(b), computed in various solvents using the nonequilibrium ptSS-PCM approach.
a CIS excited-state dipole moments are unrelaxed.

TABLE IV .
Excitation energies (ΔE) and excited-state dipole moments (μ) for the full chromophore in Fig.7(a) vs the unlinked dimer model in Fig.7(b).All calculations are performed in vacuum.
a Solvent shift relative to vacuum value.b CISD calculations using the AM1 Hamiltonian and a SCRF solvation model, from Ref. 75.c RAS-SF/6-31 * G with ptSS-PCM solvation.