{'fraction_non_alphanumeric': 0.030380169427867962, 'fraction_numerical': 0.021449735008137544, 'mean_word_length': 4.74814104101703, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Describing ground and excited state potential energy surfaces for molecular photoswitches using coupled cluster models

In this article, we use two extensively studied systems, a retinal model system and azobenzene, to explore the use of coupled cluster models for describing ground and singlet excited state potential energy surfaces of photoswitchable systems. While not being suitable for describing nuclear dynamics of photoisomerization, coupled cluster models have useful attributes, such as the inclusion of dynamical correlation, their black box nature, and the systematic improvement offered by truncation level. Results for the studied systems show that when triple excitations (here through the CC3 model) are included, ground and excited state potential energy surfaces for isomerization paths may reliably be generated, also for states of doubly excited character. For ground state equilibrium cis‐ and trans‐azobenzene, the molecular geometry and basis set is seen to significantly impact the vertical excitation energies for the two lowest excited states. Efficient implementations of coupled cluster models can therefore constitute valuable tools for investigating photoswitchable systems and can be used for preliminary black box studies to gather information before more complicated excited state dynamics approaches are pursued.


| INTRODUCTION
Photoswitches constitute a class of molecular machines that undergo reversible isomerization upon stimulus by light, so-called photoisomerization. Computational approaches for photoisomerization processes are considered to be of great value for the study of photoswitches, especially in conjunction with spectroscopic techniques. [1][2][3] To get a more complete overview of the isomerization mechanism itself and the interplay of several electronic states, accurate ab initio molecular dynamics are needed. However, the computational costs of such calculations are high, and they are complicated by the frequent occurrence of conical intersections, 4 which require special attention to the chosen method. From ground and excited state electronic structure calculations, the mechanism itself cannot be fully elucidated, but electronic structure theory provides valuable information about the potential energy surfaces and may indicate possible reaction paths for the isomerization process.
Many features of the potential energy surface are important for the machinery of photoswitches. Ground state surfaces give information on relative stability of isomers, relaxation processes following decay from excited state and energy barriers for thermal isomerization. The energy barriers will indicate whether a molecule may function as a photoswitch at all (high enough energy barrier in one direction) and if the process is thermally reversible. Excited states naturally play a key role in the mechanisms of photoisomerization, with curvature and shapes giving information on the efficiency of possible decay pathways to the ground state and whether the process is photochemically reversible. The energy of the excited states relative to the ground state energy will give information on wavelengths needed for initiating (and for certain photoswitches, reversing) the isomerization process. Mapping out potential energy surfaces for ground and excited states can also yield information on and characterize conical intersections, depending on the choice of method.
The coupled cluster Jacobian, which is diagonalized to obtain excited states, is nonsymmetric and may become defective at conical intersections. 53 This precludes the standard models from being used for ab initio excited state dynamics of photochemical systems. However, coupled cluster has advantages for the characterization of ground and excited state surfaces, such as the inclusion of dynamical correlation and the possibility for systematic improvement offered by choice of truncation level. In addition, coupled cluster models are black box in nature and the user does not have to define active spaces or parameters. However, due to the computational cost of including triple (or approximate triple) excitations to the hierarchy, there are few such results for potential energy surfaces of photoswitches in the literature. A notable exception is the use of EOM spin-flip (SF) CCSD with perturbative triples corrections (EOM-SF-CCSD(dT) and EOM-SF-CCSD(fT)) for the penta-2,4-dieniminium cation (PSB3) model system of retinal, which shows excellent agreement with MRCISD+Q results. 54 The spin-flip approach exploits the fact that a high-spin triplet reference is well behaved across the bond torsion and generates singlet ground and excited states as spin-flip excitations from the triplet reference. From this study, it is clear that inclusion of triples leads to a significantly better agreement with MRCISD+Q compared to EOM-SF-CCSD; however, it is unclear whether the high-spin reference is necessary to obtain this accuracy. CC3 has also been used for the PSB3 system, 44 but without comparison to high-level multireference methods.
For PSB3, the extensive set of MRCISD+Q calculations by Olivucci and collaborators 27 have proven indispensable for benchmarking other computational methods. 54,55 Comparisons have been made between these MRCISD+Q results and EOM variants of coupled cluster such as EOM-CCSD, EOM-SF-CCSD, EOM-SF-CCSD (dT), and EOM-SF-CCSD(fT). 54 Based on this comparison, EOM-SF-CCSD(dt/ft) was shown to agree well with MRCISD+Q. The importance of dynamical electron correlation has also been shown for PSB3, where CASSCF gave a barrierless first singlet excited state (S 1 ) potential energy surface while CASPT2 yields an S 1 surface with a small barrier. 56 Introduction of dynamical correlation by performing single point calculations along the CASSCF S 1 surface did not correct the shape of the potential energy surface. In addition, Send et al. 44 investigated excited state surfaces generated by rotation around all bonds in PSB3 using TDDFT, CC2, and CC3. From their results, it appears that there are significant differences between the coupled cluster methods and TDDFT, particularly for rotations around single bonds.
Another photoswitch that has garnered much attention is azobenzene and its derivatives. Due to the size of azobenzene, limitations are imposed on choice of computational method and/or basis set. Zhu et al. 31,32 optimized equilibrium states, transition states and conical intersections with five state averaged (5SA) CASSCF (6,6)/6-31G, 57 for the ground state and the first three singlet excited states. The focus of the work in Reference 32 was primarily to carry out molecular dynamics, however, a rotational path was generated by interpolation, which can be used to compare other computational methods with CASSCF. It should be noted that the CASSCF excitation energies of the cis and trans structures in Refs. 31 and 32 deviate significantly from experimental values, which can be explained by the fact that CASSCF excitation energies have been found to have a large deviation from more accurate methods. 58 Recently, a thorough mapping of the ground state and first excited state potential energy surfaces of azobenzene was performed by Aleotti et al. 33 using RASPT2.
In this case, the excitation energies of the cis and trans structures are in much better agreement with experimental values. This shows the importance of dynamical electron correlation when computing excitation energies for comparison with experimental results. Casellas et al. 30 further argue that dynamical electron correlation must also be included for geometry optimization of critical points on the lowest singlet excited states. In addition to the studies mentioned, there are several studies of azobenzene using CASSCF and other multireference techniques containing dynamical electron correlation. 3,18,[21][22][23]25,28,59 With respect to single-reference methods, CC2 and CCSD calculations of excitation energies have been performed for azobenzene by Hättig and coworkers. 38,39 However, high-level single-reference calculations with at least triples corrections are not presented in literature for azobenzene, and hence experimental results 60 provide the basis for evaluating computational results.
The focus of this article is to explore the capabilities of the coupled cluster models for mapping singlet potential energy surfaces of photoswitches. We utilize the newly released and efficient coupled cluster code in the eT program 61 to generate potential energy surfaces for CC2, CCSD, and CC3. As already mentioned, CC2 and CCSD will not give good results for excitations of strong doubly excited character, but for single-excitation dominated processes they are expected to give good results. 62 In this study, the performance of CC2, CCSD, and CC3 are presented for torsional angles that include potentially problematic regions such as isomeri- The article is structured as follows. In Section 2.1, we provide computational details. In Section 2.2, we present potential energy surfaces of PSB3 generated by CC2, CCSD, and CC3 and compare to MRCISD+Q 27 and spin-flip coupled cluster methods. 54 Further, we explore the basis set effect for the potential energy surfaces for PSB3. In Section 2.3, we explore how coupled cluster model, basis set, molecular geometry, and solvent effects affect vertical excitation energies for cis and trans ground state equilibrium azobenzene structures. Finally, CC3 potential energy surfaces for isomerization through CNNC torsion for azobenzene are presented in Section 2.4. Concluding remarks are given in Section 3.

| Computational details
In this section, the software and settings used for computations are presented. CC2, CCSD, and CC3 excitation energies were computed with the eT program. 61 All coupled cluster calculations presented are computed using frozen core. The decomposition of the electron repulsion integrals 65 threshold was set to 10 À10 . The energy and gradient thresholds for Hartree-Fock and the energy and residual thresholds for the coupled cluster ground state were all set to 10 À8 . The energy and residual threshold for the coupled cluster excited states were set to 10 À4 . In coupled cluster calculations where the solvent effects of methanol were accounted for by a polarizable continuum model (PCM), 66 the static permittivity was set to 32.63, the optical permittivity was set to 1.758 and the probe radius was set to 1.855 Å.

| The penta-2,4-dieniminium cation (PSB3)
In this section, we discuss CC2, CCSD, and CC3 results for CASSCF optimized penta-2,4-dieniminium cation (PSB3) structures presented by Olivucci and collaborators. 27 The three relevant isomerization paths are depicted schematically in Figure 1 To assess the shape of surfaces generated by coupled cluster models, the nonparallelity value is computed for the S 0 and S 1 surfaces for all paths. The nonparallelity value is given by, where ΔE max and ΔE min is the maximum and minimum energy difference for a potential energy surface computed in two different man-

| BLA path
In this section, we consider results for the BLA path using 6-31G*.
CC2, CCSD, and CC3 results are plotted in Figure 2 together with the MRCISD+Q benchmark results of Reference 27, and energies are relative to cis-PSB3. Additionally, nonparallelity values of the coupled cluster models relative to MRCISD+Q are presented in Table 1. From Figure 2, we first of all see that CCSD cannot describe the BLA path. The deviation in vertical excitation energies relative to MRCISD+Q is found to vary between 0.13 and 0.28 eV.
The conical intersection cannot be seen for the CCSD surfaces, as it is shifted beyond TS DIR and outside the range of the plot. However, the nonparallelity values for the S 0 and S 1 surfaces for the BLA path computed with CCSD are only of 0.02 eV and 0.01 eV (see Table 1 Hence, we see that accurate results may be obtained without using a well-behaved high-spin triplet reference.

| CT path
In this section, we consider results for the CT path using 6-31G*.
CC2, CCSD, and CC3 results are plotted in Figure 3 together with the MRCISD+Q benchmark results of Reference 27, and energies are relative to cis-PSB3. From Figure 3, we see that CCSD provides too large excitation energies, as was the case for the BLA path. The difference in vertical excitation energies between CCSD and MRCISD+Q is in the range 0.23-0.28 eV, which is similar to that of the BLA path (0.13-0.28 eV). However, CCSD is doing a worse job reproducing the F I G U R E 2 S 0 and S 1 energies for BLA-path relative to cis-PSB3 computed with CC2, CCSD and CC3 with 6-31G*. MRCISD+Q with 6-31G* benchmark results from Reference 27 are included for comparison F I G U R E 1 Schematic representation of the CT and DIR path, which are minimum energy paths in the cis ! trans isomerization, and the BLA path, which is a path that interpolates the transition states of the CT and DIR paths shapes of the S 0 and S 1 surfaces than for the BLA path, as can be seen from the nonparallelity values (see Table 1

| DIR path
The results for the DIR path for CC2, CCSD and CC3 using 6-31G* are plotted in Figure 4 together with the MRCISD+Q benchmark results of Reference 27, and energies are relative to cis-PSB3.
From Figure 4 we see that the CC2, CCSD, and CC3 results for the DIR path show much of the same characteristics as the CT path.
CCSD displays too large vertical excitation energies, deviating between 0.22-0.28 eV, and surfaces, which deviate from the shapes of the MRCISD+Q surfaces (see Table 1).

| Basis set considerations
Most results presented for PSB3 in the literature are generated using 6-31G*. In the previous sections, we presented the performance of CC2, CCSD, and CC3 relative to MRCISD+Q, all using 6-31G*. We note that basis set effects for MRCISD+Q was explored in Reference 27 for the BLA path using Pople basis sets. In this section, we study the basis set dependence for coupled cluster calculations. We choose to extend from using 6-31G* to the correlation consistent basis set aug-cc-pVDZ rather than staying within the Pople basis set hierarchy, since Dunning's correlation consistent basis sets are seen to converge faster with respect to correlation energy than the Pople basis sets. 67 In Figures 5-7, we present the 6-31G* coupled cluster results from  Table 2, the nonparallelity values of the augcc-pVDZ calculations relative to the 6-31G* coupled cluster calculations are presented and demonstrate how much the shape of the surface is affected by choice of basis set. We see from

| Exploring vertical excitation energies for ground state equilibrium cis-and trans-azobenzene
In this section we investigate how choice of coupled cluster method and basis set, as well as accounting for solvent, affects computed S 0 !S 1 and S 0 !S 2 vertical excitation energies for the azobenzene molecule (see Figure 8) CC2 results are computed using cc-pVDZ, aug-cc-pVDZ, cc-pVTZ and aug-cc-pVTZ, whereas CC3 results are computed using cc-pVDZ and aug-cc-pVDZ. For calculations using aug-cc-pVDZ and aug-cc-pVTZ augmentation is only added to non-hydrogen atoms (using cc-pVDZ and cc-pVTZ on hydrogen). The vertical excitation energies are presented in Table 3 together with experimental UV-VIS results. 60 A discussion in connection to the experimental results is provided at the end of this section.
We first compare CC2 and CC3 results. The absolute difference between CC2 and CC3 S 0 ! S 1 excitation energies when the same basis set is used is 0.00-0.02 eV, whereas the absolute difference between CC2 and CC3 S 0 !S 2 excitation energies is 0.04-0.10 eV.
The slightly larger difference between CC2 and CC3 for the S 0 !S 2 transition relative to the S 0 !S 1 transition, may be due to the S 2 state having a lower fraction of single amplitudes and hence CC3 provides a better description. The fraction of singles in the CC3 calculations for cis a and trans a are listed in Table 4 (together with fraction of singles for the rotational path discussed in Section 2.4). However, the differences between CC2 and CC3 are still small for the structures studied in this section, and this is consistent with available benchmark literature. 62 From Table 3, we see that the choice of basis set impacts excitation energies for S 0 !S 1 differently than S 0 !S 2 for both sets of structures. Improving the basis set from cc-pVDZ to aug-cc-pVDZ lowers  Note: The basis set used for CC2 are cc-pVDZ, cc-pVTZ, aug-cc-pVDZ and aug-cc-pVTZ (cc-pVDZ and aug-cc-pVDZ for CC3). In addition, energies from UV-VIS spectroscopy performed in methanol 60 for cisand trans-azobenzene are presented.
From the results in Table 3 we can also discuss excitation energies for two different cis-and trans-geometries, that is, 5SA-CASSCF (6,6)/6-31G ground state optimized structures from Reference 32 (cis a and trans a ) and MP2/aug-cc-pVDZ ground state optimized structures (cis b and trans b ). A notable difference between the cis-structures is the CNNC dihedral angles, which is 9.5 for cis a and 6.0 for cis b .
For the trans structures a notable difference is the NNC bond angles, which are 117.9 and 118.2 for trans a and 113.3 for both angles for trans b . For simplicity, we only discuss CC2/aug-cc-pVTZ results, but we note that the trends observed for CC2/aug-cc-pVTZ is seen for all basis sets and also for CC3 results. For cis a and cis b , the differences in excitation energies are 0.38 eV for S 0 !S 1 and 0.10 eV for S 0 !S 2 . For the trans a and trans b structures, the differences in excitation energies are 0.20 eV for S 0 !S 1 and 0.36 eV for S 0 !S 2 . The differences between the sets of cis-and trans-structures are, therefore, significant. When we proceed to discuss and compare against experimental results we will only consider the MP2/aug-cc-pVDZ structures (cis b and trans b ). We expect MP2/aug-cc-pVDZ optimized structures to be more representative gas phase ground state structures than the ones generated using 5SA-CASSCF(6,6)/6-31G (cis a and trans a ) due to the inclusion of dynamical correlation and the use of a larger basis set.
Here, we compare computed results against absorption maxima from UV-VIS spectroscopy in methanol taken from Reference 60, which are also presented in Table 3.  68 showed an average deviation of 0.11 ± 0.08 eV between CC2 vertical excitations and absorption maxima computed using a nuclear ensemble approach at the CC2 level for a set of 28 molecules (60 transitions) to account for nuclear rearrangement.
Based on this, as well as errors related to the methods used, 62  Another factor which must be considered when comparing vertical excitation energies with UV-VIS absorption maxima is the solvent effect. In Table 5 for these methods and basis set, is that the S 0 !S 1 energies increase by 0.04-0.08 eV, while the S 0 !S 2 transitions are lowered by 0.04-0.09 eV.
The inclusion of solvent effects through PCM, therefore, does not appreciably improve the agreement between the computed vertical excitation energies and absorption maxima.

| Coupled cluster results for a rotational path of azobenzene including an excited state with doubly excited character
In this section, we discuss CC3 results for a cis to trans CNNCrotational path obtained from Zhu et al. 32 The path consists of the cis and trans ground state structures, a twisted minimum energy structure optimized for the second excited state (rot-S2), all optimized with 5SA-CASSCF(6,6)/6-31G, and structures generated by linear interpolation of internal coordinates (LIIC) with respect to these three structures (for details see Reference 32). In Figure 9, CC3 results using the cc-pVDZ basis set are presented for selected structures along the described rotational path for the ground state and the states corresponding to the two lowest singlet vertical excitation energies.
The maximum of S 0 corresponds to the aforementioned rot-S2 structure. To describe these states accurately, either cc-pVTZ, aug-cc-pVDZ or a bigger basis set should ideally be used, as can be seen from the basis set study of cis-and trans-azobenzene presented in Section 2.3. However, cc-pVDZ is sufficient to evaluate whether CC3 can describe these states.
Comparing to the results of Reference 32, it is seen that CC3 and 5SA-CASSCF (6,6) are in qualitative agreement for the ground and lowest excited state. However, for S 2 , we see from Figure 9  dihedral angle) around rot-S2, the n 2 π *2 state is the second excited state. The significant doubly excited character of the second excited state in this region is seen from the relative importance of singles amplitudes in the CC3 calculation. In Table 4 the fraction of singles amplitudes (jR 1 j/jRj) is shown for the first and second excited state throughout the rotational path. As ππ * and n 2 π *2 switch ordering,

| CONCLUSIONS
The results for PSB3 showed that vertical excitation energies computed with CC2 and CCSD deviated from MRCISD+Q results by 0.02-0.27 eV and 0.13-0.28 eV, respectively. The nonparallelity values also showed that CCSD generally gave deviation in the shape of the surfaces relative to MRCISD+Q, and while the nonparallelity values were small for the CC2 ground states, larger nonparallelity values were found for the CC2 excited states. Contrary, CC3 provided excellent agreement with MRCISD+Q, with vertical excitation energies differing by 0.01-0.04 eV.
As a consequence, the nonparallelity values for all ground and excited states were low. The performance of CC3 was similar to EOM-SF-CCSD (dT) presented in Reference 54, where it was seen that triples corrections are necessary. However, it appears that the use of a well-behaved triplet reference state was not required for the PSB3 surfaces explored here.
For equilibrium geometries of cis-and trans-azobenzene, a study of S 0 !S 1 and S 0 !S 2 vertical excitation energies computed by coupled cluster models was conducted. It was shown that the molecular geometry had the largest impact, yielding significant variations for vertical excitation energies and shifting the energy of S 1 and S 2 relative to each other. Further, a basis set of sufficient quality was required to describe these states. Both cc-pVTZ and aug-cc-pVDZ were shown to significantly change the vertical excitation energies relative to cc-pVDZ. However, aug-cc-pVDZ was shown to be the better choice of the two, since the inclusion of diffuse functions was found to be important to correctly describe S 0 !S 2 transitions.
F I G U R E 9 Ground state (S 0 ) and the two lowest excited singlet states (S 1 , S 2 ) calculated with CC3 using cc-pVDZ, for a selected set of azobenzene structures along the rotational path from Reference 32 Conversely, the differences between CC2 and CC3 results were shown to be small. By comparing computed vertical excitation energies to UV-VIS absorption maxima, CC2/aug-cc-pVTZ results were found to agree well for all excitations, except the S 0 !S 1 transition of cis-azobenzene. The computed vertical excitation energies were seen to be sensitive to molecular geometry, and this could in part explain the discrepancy between experimental and computational results. In the case of the CNNC rotational path from Reference 32, we presented the S 0 , S 1 , and S 2 surfaces using CC3. The CC3 S 0 and S 1 surfaces were in qualitative agreement with the CASSCF surfaces, while the S 2 surface differed significantly since ππ * and n 2 π *2 switched ordering. This was seen from the onset of significant doubly excited character of the S 2 state in regions of the path. Hence, CC2 and CCSD would not be able to describe the region of the S 2 surface surrounding rot-S2, where n 2 π *2 was the second excited state.
Thus, the article has demonstrated that the coupled cluster models, despite being predominantly single-reference, constitute valuable tools for preliminary and simple black box studies of photoswitchable systems.