The effect of polarizable environment on two-photon absorption cross sections characterized by the equation-of-motion coupled-cluster singles and doubles method combined with the effective fragment potential approach
Solvent can significantly affect nonlinear optical properties of solvated chromophores. We report an extension of a hybrid polarizable embedding method incorporating solvent effects in the electronic structure calculations of two-photon absorption (2PA) cross sections. The approach uses the equation-of-motion coupled-cluster singles and doubles method for excitation energies (EOM-EE-CCSD) for the quantum region (QM) and the effective fragment potential (EFP) method for the classical region (MM). The impact of the environment on the 2PA cross sections is investigated in microhydrated clusters of para-nitroaniline, thymine, and deprotonated anionic chromophore of photoactive yellow protein (PYPb). The performance of EOM-EE-CCSD/EFP is assessed by comparing the 2PA cross sections against full quantum calculations as well as against non-polarizable QM/MM electrostatic embedding approach. The EOM-EE-CCSD/EFP approach reproduces well the main trends in the 2PA cross sections across different clusters and its performance improves when few water molecules are included in the QM subsystem. For the studied transitions with the strong charge-transfer character, the errors in the 2PA cross sections correlate with the errors in the excitation energies, transition moments, and dipole-moment differences between the initial and final states. Analysis of the 2PA transitions in terms of natural transition orbitals is also reported. The analysis of 2PA transition densities allows us to introduce a a rigorous metric for assessing the domain of applicability of QM/MM methods for 2PA cross sections calculations. Finally, the calculations of 2PA cross section of PYPb in bulk water are used to assess the convergence of the 2PA cross sections with respect to the number of water molecules included in the QM subsystem. Relative to the small QM (chromophore only) calculations, the inclusion of water molecules from the first solvation shell leads to a large change (∼22%) in the calculated 2PA cross section. The convergence of the key properties, such as the excitation energies, transition moments, and dipole-moment differences with the size of the QM subsystem, is investigated using configuration interaction singles method. I. INTRODUCTIONAb initio electronic structure methods enable highly accurate calculations of molecular structures, energetics, and properties. The domain of applicability of these approaches in terms of the size of a system is, however, limited by their steep scaling and large computational costs. Despite significant progress in computer hardware and algorithms, full quantum-chemical modeling of condensed-phase systems remains elusive. This limitation can be circumvented by using multi-scale approaches combining high-level quantum mechanical (QM) description of a system (i.e., a solvated chromophore) with a more approximate treatment, such as classical molecular mechanics (MM), of the environment. Among many different flavors, explicit solvent approaches are particularly attractive, because, in contrast to implicit solvent models (polarizable continuum), they are capable of describing specific interactions such as solvent-solute hydrogen bonding and salt bridges.Numerous hybrid multi-scale methods 1 differ in their description of the interactions between the QM subsystem (treated at a high level of theory) and the environment. In the most commonly used QM/MM approach, classical force fields are used to describe the environment (hence, MM) and the interaction between the QM and MM parts is described by electrostatic embedding 2 . This basic QM/MM strategy has been successfully used for describing not only the ground-state energetics and properties of large systems, but also photochemical properties 3-6 . The basic QM/MM scheme can be improved by using polarizable force fields in which the charge distributions in the MM part can respond to the changes in electronic distribution in the QM part.The quality of force fields (whether polarizable or not) is essential for obtaining reliable electronic properties of the QM subsystem. Force fields are commonly derived using empirical parameterization aiming at reproducing particular observables. A more rigorous
<!>II. THEORY<!>B. The EOM-EE-CCSD/EFP method<!>C. 2PA within the EOM-EE-CCSD/EFP scheme<!>D. Summary of the algorithm<!>34.<!>III. COMPUTATIONAL DETAILS<!>IV. RESULTS AND DISCUSSION<!>B. Microsolvated p-nitroaniline clusters<!>C. Microsolvated thymine clusters<!>D. Microsolvated clusters of the anionic PYP chromophore<!>V. CONCLUSIONS
<p>description can be achieved by using force fields derived from first principles, such as in the effective fragment potential (EFP) method [7][8][9][10][11] , wherein the MM subsystem is fragmented into smaller subunits (fragments) and the force-field parameters of each of these fragments are obtained from ab initio calculations. In the EFP approach, the fragments are frozen (i.e., their structures do not change), which enables using a fixed set of parameters computed once. In principle, these parameters can be computed "on the fly", but this, of course, leads to increased computational costs. The QM/EFP scheme has been shown to yield accurate solvatochromic shifts, ionization and detachment energies, electron-attachment energies, as well as redox potentials in a variety of solvated systems including protein-bound chromophores [12][13][14][15][16][17][18][19] . EFP is similar to polarizable embedding (PE) approach [20][21][22][23] . The description of electrostatics and polarization is nearly identical in EFP and PE; however, EFP features more rigorous treatment of dispersion and exchange repulsion.</p><p>Higher-order response properties such as two-photon absorption (2PA) cross sections are particularly sensitive to the environment 24 . For example, the 2PA cross sections of red fluorescent proteins with identical DsRed-like chromophores (DsRed, mCherry, and mStrawberry) differ by a factor of five 25 . Non-linear response properties are formally given by sum-over-states expressions, which depend on the wave functions and excitation energies of all electronic states of the system. Thus, in contrast to one-photon transition moments expressed as matrix elements between the initial and the final states only, the quality of description of non-linear properties, such as 2PA cross sections, depends on the quality of the representation of the entire spectrum by a model Hamiltonian. This is the primary reason why modeling 2PA cross sections is challenging. QM/MM approach has a number of limitations that make reliable calculations of 2PA cross sections difficult. First, QM/MM schemes with only one chromophore in the QM part cannot describe delocalized states, such as those in molecular aggregates 26 or charge-transfer states delocalized over several solvent molecules 27 . This issue can be partially remedied by a careful choice of the QM subsystem, for example, by including some solvent molecules.</p><p>Second, small errors in the QM/MM excitation energies (present in the denominator of the formal expression) and transition dipole moments can accumulate leading to large errors in the calculated 2PA cross sections. Third, the quality of computed 2PA cross sections is sensitive to the degree of electron correlation included in the model Hamiltonian.</p><p>In the context of 2PA calculations, the EFP method has been employed in combination with time-dependent density functional theory (TDDFT) 28 , but not with high-level wavefunction methods, such as equation-of-motion coupled-cluster singles and doubles method for excitation energies (EOM-EE-CCSD) [29][30][31] . Coupled-cluster methods have been employed within polarizable embedding scheme in calculations of 2PA cross sections 22 . In this article, we extend our EOM-CCSD implementation of the 2PA cross sections [32][33][34] to include the effect of solvent environment via QM/EFP embedding. We assess the quality of the EOM-EE-CCSD/EFP 2PA cross sections (i.e., computed with the EOM-EE-CCSD wave functions embedded in the polarizable environment of EFP fragments) relative to the full EOM-EE-CCSD calculations and identify the main sources of errors. We consider low-lying transitions in small microsolvated clusters of three chromophores, para-nitroaniline (pNA), thymine, and deprotonated anionic chromophore of photoactive yellow protein (PYPb). To analyze the role of solvent polarization, we compare the EOM-EE-CCSD/EFP 2PA cross sections for the microhydrated clusters with the results from two non-polarizable QM/MM schemes:</p><p>(a) EOM-EE-CCSD/EFP without polarization effects and (b) EOM-EE-CCSD/CHARMM, wherein the solvent molecules are represented by CHARMM point charges 35 . To assess the effect of bulk solvation, we calculate the 2PA cross sections of the PYPb chromophore in water and investigate the convergence of the results with respect to the QM size by including an increased number of water molecules into the QM subsystem.</p><!><p>A. The QM/EFP scheme In the QM/EFP embedding scheme, the Hamiltonian of the whole QM/MM system is written as</p><p>where H QM and H M M describe the QM and MM parts, respectively, and H QM/M M describes the interaction between these two subsystems. In the QM/EFP approach, the MM system is broken into effective fragments and its energy, E M M , is computed as the sum of manybody electrostatic, polarization, dispersion, and exchange-repulsion interactions between the fragments:</p><p>In calculations of solvated species, each solvent molecule is treated as an EFP fragment.</p><p>The EFP method can also be applied to macromolecules via a fragmentation scheme that can deal with breaking covalent bonds 17 .</p><p>In polar and hydrogen-bonded solvents such as the ones studied here, electrostatics and polarization are the leading terms for inter-fragment interactions. The electrostatic interaction energy, E elec , is computed using Stone's distributed multipole analysis [36][37][38] by placing permanent multipoles (up to octopoles) at the atoms and bond midpoints. Only the interactions of the induced dipoles with the permanent multipoles and induced dipoles on other fragments are included in the calculation of E pol . Detailed description of the dispersion and exchange-repulsion terms as well as the damping functions needed to correct the artificial charge-penetration effects and to prevent "polarization catastrophe" at close separations of the permanent and induced multipoles can be found elsewhere 9,10,39-41 .</p><p>The QM system is embedded in the polarizable environment such that the electrostatic and polarization interactions between the QM and MM subsystems are included explicitly in the one-electron part of H QM :</p><p>where φs are the molecular orbitals. V elec is the Coulomb potential due to the nuclear charges and permanent multipoles of the fragments, whereas V pol is the potential due to the induced multipoles. The induced dipoles on each fragment are calculated self-consistently with each other and with the Hartree-Fock (HF) wave function using a two-level iterative scheme.</p><p>The converged induced dipoles of the reference HF determinant are fixed during subsequent calculations of the ground-and excited-state wave-function amplitudes (e.g., T amplitudes in CC and excitation amplitudes in EOM-CC). The contribution to the polarization energy due to the interactions of the induced dipoles with the QM electronic density, E pol,QM/M M , is included in the state energies. Since the induced multipoles calculated at the HF level are kept fixed in the calculations of the excited-state wave functions, their contribution to the excitation energies is the same for all excited states. In other words, beyond the HF level the polarizable environment is not allowed to respond to the change in the electronic wave function. State-specific energy corrections due to the response of polarizable environment to different electron distributions of each state can be computed perturbatively 12 and were shown to improve the quality of computed energy differences 13,14,[16][17][18][19] . Here,we omit these corrections because this state-specific treatment cannot be easily incorporated in the response equations solved in 2PA calculations. In the current implementation 12 , the dispersion and exchange-repulsion interactions between the QM and MM systems (E disp,QM/M M and E exch,QM/M M ) are treated approximately, as interactions between the EFP fragments, and are added to the state energies, so they do not contribute to wave functions, excitation energies, and other transition properties. The total energies of the ground and excited states and the excitation energy can be written as</p><p>where Ψ gr and Ψ ex are the ground-and excited-state wave functions (in CIS, TDDFT, and EOM-CC calculations, excitation energies are computed directly, and not as energy differences from state-specific calculations).</p><!><p>Our choice of the electronic structure method for the QM system is EOM-EE-CCSD 29,42,43 in which the excited-state wave functions are given by</p><p>Here, Φ 0 is the reference determinant and e T |Φ 0 > is the CCSD wave function. The CCSD cluster operator T and the EOM-CCSD operator R comprise single and double excitation operators:</p><p>The amplitudes of the operator T satisfy the following equations:</p><p>where H = e −T He T is the similarity-transformed Hamiltonian and Φ µ denote singly and doubly excited determinants. In the EOM/EFP scheme, the one-particle component of the Hamiltonian H is given by Eq. ( 3). The EOM amplitudes R and the corresponding energies for EOM-CCSD states are obtained by diagonalizing H in the space of the reference, singly, and doubly excited determinants:</p><p>Since H is non-Hermitian, its left (</p><p>are not Hermitian conjugates, but form a biorthonormal set</p><p>An important feature of the EOM-CCSD/EFP scheme is that it preserves biorthonormality of the target EOM states, because the induced multipoles are computed self-consistently at the HF step and kept fixed during the calculations of the CCSD and EOM amplitudes.</p><p>This biorthogonality is important for computing transition properties such as 1PA and 2PA transition moments.</p><!><p>We compute the 2PA transition moments and cross sections for EOM-CCSD wave functions following the approach described in Refs. 32 and 34. Formally, the 2PA transition moments are given by the following sum-over-states expressions:</p><p>Here, ω 1 and ω 2 are the frequencies of the two absorbed photons satisfying the 2PA resonance</p><p>To derive expressions for the EOM-CCSD 2PA transition moments, we replace the dipole-moment operator, µ a , with the similaritytransformed dipole-moment operator, μa = e −T µ a e T and the wave functions, Ψ k |s and |Ψ k s, with the EOM-CCSD target states, Φ 0 L k |s and |R k Φ 0 s, in the above expressions.</p><p>Then, we recast these expressions into the numerically and formally equivalent form as follows:</p><p>where X and X are the first-order response wave functions computed by solving auxiliary response equations. This reformulation of the sum-over-states expressions leads to practical expressions of the 2PA moments in terms of just the zeroth-and first-order wave functions of the initial and the final states. The response equations are:</p><p>where Φ I s denote Slater determinants from the target-state manifold (e.g., the reference, singly excited, and doubly excited determinants for EOM-EE-CCSD).</p><!><p>The EOM-EE-CCSD/EFP method is implemented in Q-Chem quantum chemistry package 44,45 . The procedure for calculating 2PA cross sections involves the following steps:</p><p>1. First, using a standard procedure 9,10 , all MM solvent molecules are replaced by rigid effective fragments placed at the same position and orientation as the original solvent molecules. The QM subsystem is now in the field of these effective fragments.</p><p>2. The two-level iterative HF procedure is then employed such that the induced dipoles on each EFP fragments are self-consistently iterated with the HF wave function and with induced dipoles of other EFP fragments. The converged induced dipoles are then kept fixed in the subsequent steps of the calculation.</p><p>3. The CCSD T -and Λ-and EOM-EE-CCSD R-and L-amplitudes and energies are then computed using the HF wave function and orbitals calculated in the second step.</p><p>Using these amplitudes, the CCSD and EOM-EE-CCSD state and transition densities are formed for calculating the state and transition dipole moments.</p><p>4. In the 2PA step, the excitation energies without the polarization response of the environment to the wave functions of the ground and excited states are used to determine the photon frequencies according to the 2PA resonance condition. The procedure for calculating the EOM-EE-CCSD/EFP 2PA cross sections is described in Refs. 32 and</p><!><p>E. Characterization of 1PA and 2PA transitions in terms of natural transition orbitals</p><p>Valuable insight into solvent effects can be gained by understanding electronic transitions (1PA, 2PA, etc.) in terms of molecular orbitals. This entails mapping the changes in electronic density induced by excitation onto pairs of molecular orbitals. For 1PA transitions, such maps are provided by the reduced transition one-particle density matrices (1PDM), γ:</p><p>where p † and q create and annihilate electrons in molecular orbitals φ p and φ q , respectively.</p><p>Importantly, γ is related to the experimental observables, e.g., the 1PA transition moment is given by</p><p>The elements of transition 1PDM can be interpreted as amplitudes of an excitation operator that generates the final state from the initial correlated state:</p><p>The squared norm of γ:</p><p>gives the weight of one-electron character of the transition and provides a bound to the respective expectation values by virtue of the Cauchy-Schwarz inequality 46,47 .</p><p>Transition 1PDM can be used to define an exciton's wave function (Ψ exc ), a two-particle quantity that contains essential information about the changes in electron density upon the transition [48][49][50][51] :</p><p>where r H and r P are the hole and particle (electron) coordinates (using r H = r P = r, Ψ exc is reduced to the transition density).</p><p>Eqs, ( 20) and ( 22) allow one to interpret the individual elements of γ as weights of particular configurations. For example, using localized orbitals, one can express the probability of finding the exciton in a spatial domain D as:</p><p>As discussed in the next subsection, this reasoning can be extended to quantify the extent of charge-transfer character 52,53 .</p><p>Further simplification in exciton's representation can be achieved by applying singularvalue decomposition (SVD) to γ. The SVD procedure reduces the description of electronic excitation to one-electron transitions between pairs of orbitals called Natural Transition Orbitals 51,53-56 (NTOs). Using NTOs, the exciton's wave function can be written as:</p><p>where φK,P and φK,H are the particle and hole orbitals obtained by SVD of γ and σ K are the corresponding singular values. Usually, only a few σ K s are significant, giving rise to simple interpretation of excited-state characters in terms of excitations between hole and particle orbitals. By using NTOs, one can express the inter-state matrix elements between manybody wavefunctions in terms of matrix elements between orbitals, facilitating the connection between molecular orbital theory and experimental observables:</p><p>A detailed description of NTO analysis and its implementation in Q-Chem is given in Refs.</p><p>51 and 52.</p><p>We recently extended the concepts of transition 1PDMs and NTOs to 2PA transitions 34 .</p><p>In addition to enabling visualization of a 2PA transition in terms of molecular orbitals, the 2PA NTO analysis provides insight into the character of the virtual state involved and the type of 2PA transitions.</p><p>F. Conditions necessary for using fragment approach for describing 2PA transitions in condensed phase</p><p>Let us now discuss the domain of applicability of the fragment approaches for modeling 2PA. With an aim to provide a metric quantifying whether the conditions needed for QM/EFP treatment are met by a particular 2PA transition, we invoke the NTO analysis of the respective perturbed 1PDM. As per Eq. ( 25), the 2PA moments can be expressed in terms of matrix elements between the respective hole and particle orbitals:</p><p>where φa K,P and φa K,H are the K th particle and hole NTOs for the 2PA a-component (a, b = x, y, z) transition 1PDM with singular value of σ a K . Since the NTOs are linear combinations of atomic orbitals (χs)</p><p>Eq. ( 26) can be expressed in the basis of χs as</p><p>The corresponding 2PA a-component transition 1PDM in the atomic orbitals basis is:</p><p>By partitioning the QM system into the chromophore (A) and solvent (B), we can partition the above 1PDM into the four parts, A←A γ a f ←i , A←B γ a f ←i , B←A γ a f ←i , and B←B γ a f ←i , where C←D γ represents the contribution of the hole residing on fragment D and a particle residing on fragment C:</p><p>Then the probability of finding the hole on fragment D and the particle on fragment C is given by</p><p>The four probabilities corresponding to the four fragment 1PDMs give a measure of the charge transfer within and across fragments to the 2PA transition. This is similar to the charge-transfer numbers used to quantify the extent of charge-transfer in electronic transition 51,52 . Breaking the system into a QM solute and MM solvent does not introduce large errors only when the probability of finding both the hole and particle is predominantly on the solute, i.e., when</p><p>The performance of fragment approaches such as QM/EFP deteriorates with increasing contributions of solvent to the 2PA transition 1PDMs.</p><p>We computed these CT numbers for the 2PA transitions in microhydrated clusters studied here as well as for representative snapshots from solvated PYP simulations. The analysis reveal that that the charge-transfer probabilities between the chromophore and the solvent are low ( A←B Q a f ←i + B←A Q a f ←i <0.01). The probability of finding both the hole and particle NTOs on the solvent molecules is even lower. This indicates that the QM/MM fragmenting of these systems should not introduce large errors into the 2PA cross sections, because the 2PA transition is predominantly localized on the chromophore. This encouraging finding is far from obvious, as, in contrast to 1PA transitions for which only the initial and final states need to be localized on the solute, one could expect noticeable charge-transfer character in the virtual states.</p><!><p>All calculations were carried out using the Q-Chem quantum-chemistry package 44,45 . We used the B3LYP/6-31+G* optimized geometry for pNA, RI-MP2/cc-pVTZ geometry for thymine, and CCSD/6-31+G* geometry for PYPb in the calculations of bare chromophores.</p><p>The structures used in the QM/EFP calculations for the microsolvated clusters of pNA, thymine, and PYPb were taken from Refs. 12, 57, and 58, respectively. We used the EOM-EE-CCSD/6-31+G* level of theory within three different QM/MM schemes in which the MM subsystem was represented by the EFP fragments, EFP fragments without polarization, and CHARMM point charges. Core electrons were kept frozen in all calculations. In each of these systems, we characterized the lowest symmetric 2PA transition considering degenerate photons. For pNA, we also considered non-degenerate photons, the first photon having a frequency of 25,400 cm −1 (3.1492 eV). We used TheoDORE package 50 for computing the charge-transfer probabilities following the 2PA NTO calculations with Q-Chem.</p><p>In the QM/EFP calculations, all solvent water molecules were replaced by the standard effective fragments placed at the same positions and orientations as the initial solvent molecules 9,10 . This means that the Cartesian geometries of the initial and the QM/EFP model clusters are slightly different. To eliminate the ambiguity due to slight geometry change, we used these QM/EFP geometries in all calculations of the microsolvated clusters, meaning the structures of the water molecules in extended QM subsystems are exactly the same as in the effective fragments. All Cartesian geometries are provided in the SI.</p><p>To model bulk solvation, we performed molecular dynamics simulations with a PYPb chromophore and one sodium cation in a periodic box of 1717 water molecules using NAMD 59 . We used the CHARMM22 parameters for water and sodium cation; the parameters for PYPb are given in the SI. After the initial minimization for 20 ps (with 2 fs time step), we equilibrated the system for 2 ns (with 1 fs time step) in an NPT ensemble at 298 K and 1 atm and followed up with a production run of 1 ns. From this trajectory, we picked 21 snapshots separated by 50 ps for the 2PA calculations using EOM-EE-CCSD/6-31+G*//EFP (the pdb files for the snapshots are given in SI). In the equilibrium trajectory, the sodium cation is more than 10 Å from the PYPb chromophore. In the 2PA calculations, we included only the chromophore and water molecules.</p><p>To investigate the effect of the QM size on the computed 2PA cross section, we systematically extended our QM subsystem by including explicit water molecules at the phenolate and carboxylic acid ends of PYPb. In these calculations, the structures of QM water molecules were identical to the structures in the MD snapshots. Using these snapshots, we also performed CIS/6-31+G*//EFP calculations to investigate the convergence of the key components -excitation energies, transition dipole moments, and dipole-moment differences between the ground and excited states. In these CIS/EFP calculations, the structures of the model systems were obtained by using the standard EFP procedure 9,10 on the MD snapshots, meaning that the structures of water molecules in the QM system are slightly different from the structures in the MD snapshots.</p><!><p>A. Molecular orbital framework for 2PA transitions in pNA, thymine, and PYP chromophore To understand the effect of microsolvation on the 2PA transitions of the three benchmark molecules, we begin by characterizing these transitions in the bare chromophores. We aim to understand the nature of the solute-solvent interactions by analyzing the orbital character of the transitions. We calculated the 2PA NTOs using our wave-function analysis toolkit 34 .</p><p>The leading NTO pairs and the molecular orientations are shown in Fig. 1.</p><p>In our calculations, the pNA molecule is in the xz plane, with the z axis along the dipole moment. The lowest symmetric 2PA transition has a large microscopic cross section (see Table I) and is dominated by the zz component of the 2PA transition moment. The NTO analysis of the z component 1PDM shows that this transition is well described by just one NTO pair with a large singular value. This NTO pair (shown in Fig. 1) reveals the intramolecular charge-transfer character: upon this transition the dipole moment increases due to the shift in the electronic density towards the electron-withdrawing nitro group.</p><p>A more detailed analysis of this transition can be found in Ref. 34</p><p>which reduces to</p><p>for degenerate photons and diagonal elements of the 2PA transition moment tensor 34,60,61 .</p><p>Eq. ( 33) also shows that in this case, the 2PA transition moment is proportional to the transition dipole moment ( Ψ f |µ a |Ψ i ) and inversely proportional to the excitation energy</p><p>. Thus, understanding the difference in the permanent dipole moment in the ground and excited state, transition dipole moment, and excitation energy in the microsolvated clusters is the key to understanding how solvent impacts 2PA cross sections.</p><p>For example, bulk solvation stabilizes the excited state more than the ground state because of a larger excited-state dipole moment, resulting in the lowering of the excitation energy relative to the gas phase. However, in microsolvated clusters, the change in the excitation energy depends on the relative orientation of the solute molecules and the extent to which they stabilize or destabilize the hole and particle orbitals.</p><p>In our calculations, thymine (C s ) is in xy plane. The lowest A transition with degenerate photons has a small microscopic cross section dominated by the M xx and M yy components.</p><p>The NTO analyses of the x-and y-components of the respective 1PDMs show that each of these component transitions can be described by a single pair of NTOs, which indicates that the transition is accompanied by some intramolecular charge transfer, away from the oxygen of the carbonyl group across the methyl group. This is further validated by the P R N T O s for these 1PDMs (∼1-2) and their component ωDMs (>20). So, this 2PA transition can also be explained by Eq. ( 34) and the changes in its 2PA cross section in microhydrated clusters can be explained by understanding the corresponding changes in the transition moments, differences in permanent dipoles, and excitation energy. The shift in electronic density increases dipole moment of the molecule in the excited state, which suggests red shift in aqueous solution.</p><p>PYPb also has C s symmetry and is in the xy plane in our calculations. For this chromophore, the lowest symmetric 2PA transition with degenerate photons has a large microscopic 2PA cross section with a dominant M xx component. The 2PA NTO analysis for the</p><p>x-component 2PA 1PDM gives two dominant NTO pairs, which represent opposing channels of electronic density flow upon the 2PA transition. The dominant channel reveals charge transfer from the carboxylic acid to the phenolate and the second channel reveals the charge transfer from the phenolate to the carboxylic acid. The latter channel corresponds to the orbital character of the 1PA transition between these states shown in SI (Fig. S1). Thus, the dominant NTO pair indicates that the virtual state is dominated by the initial state, which has a larger electronic dipole than the excited state 62 . In contrast, in the studied transition in pNA, the character of the virtual state is dominated by the final state with a larger dipole moment than the ground state and only one NTO pair. Quantitative analysis of the P R N T O s for this 1PDM (∼1-2) and its component ωDMs (>30) also shows that this 2PA transition in PYPb has some intramolecular charge-transfer character. Moreover, the respective transition dipole moment is also large. Eq. ( 34), therefore, explains the large 2PA cross section for this transition.</p><!><p>In all microhydrated clusters of pNA, the excitation energy calculated with the full QM approach is slightly lower than that of the bare chromophore. This indicates that the water molecules, which are clustered around the nitro group, stabilize the particle orbital of the lowest symmetric 2PA transition more than the hole orbital. Moreover, the character of the dominant 2PA NTOs for this transition does not change upon microsolvation, as expected for weakly interacting solute-solvent systems. The change in the oscillator strength is also small. Table S1 in the SI indicates that the difference in permanent dipole moments of the ground and excited states increases, which along with the slight decrease in the excitation energy -according to Eq. ( 33)-increases the 2PA cross section by up to 8% for degenerate photons and 15% for non-degenerate photons. Across the three clusters, the full QM cross section increases from dihydrate to tetrahydrate and decreases slightly from tetrahydrate to hexahydrate.</p><p>The QM/EFP scheme with just the pNA molecule in QM reproduces this trend, but slightly exaggerates the changes across these clusters, with errors of about 2-9% relative to the respective full QM cross sections, as shown in Table I. Omitting the polarization contribution in the QM/EFP calculations results in smaller errors (2-5%) relative to the full QM cross sections. This approach, however, yields slightly smaller cross sections for 2PA transition with non-degenerate photons for the hexahydrate than the dihydrate. Similarly, QM/MM(CHARMM) also yields smaller errors relative to the full QM cross sections across the three clusters, but fails to reproduce the trend obtained with the full QM calculation. In fact, QM/CHARMM shows a decrease in the cross sections for both sets of photons from dihydrate to tetrahydrate and an increase from tetrahydrate to hexahydrate.</p><p>Yet, the small differences between QM/EFP, QM/EFP with polarization turned off, and QM/MM(CHARMM) preclude us from arguing in favor of any of these approaches. Such small differences also indicate that for these systems the impact of intermolecular polarization on 2PA cross sections is not large.</p><p>We also report QM/EFP results for dihydrate and tetrahydrate clusters with the QM subsystem comprising the pNA and its nearest hydrogen-bonded water. These calculations result in smaller errors of 2-5% relative to the full QM cross sections. Similar reduction of the QM/EFP errors upon inclusion of the hydrogen-bonded waters in the QM system were reported for core ionization energies of solvated glycine 19 .</p><p>The errors in the 2PA cross sections for the three QM/MM schemes relative to the full QM results can be traced back to the slightly smaller errors in the excitation energies, transition moments, and dipole-moment differences with the QM/EFP scheme as tabulated in the SI (Table S1). We note that the errors in the excitation energies for these three QM/MM approaches relative to full QM results are very small (<1%), but the errors in the transition dipole moments and dipole-moment differences are slightly larger (up to 2-3% for both).</p><!><p>The orbital character of the lowest A 2PA transition in thymine does not change upon microsolvation. In the full QM calculation, the microsolvation changes the excitation energy by up to 3%, oscillator strength by up to 14%, and 2PA cross sections by up to 20% (we note that the magnitude of the microscopic cross section is small). In the full QM calculation, the excitation energy of this transition decreases in T1 to T2 relative to the bare thymine (see Table II), but in the T3 cluster, the excitation energy is higher than in bare thymine.</p><p>The difference in trends can be attributed to the following structural feature: only in the T3</p><p>cluster the water molecule is near the carbonyl group that is across from the methyl group, thereby stabilizing the hole orbital shown in Fig. 1b. All three QM/MM schemes reproduce this trend in the excitation energy across the clusters with negligible errors (<1%) relative to the full QM calculation. The 2PA cross section decreases from T1 to T2 to gas-phase thymine to T3. All three QM/MM schemes, however, yield smaller 2PA cross sections for the T1 cluster than for the T2 cluster, with comparable errors. For the T3 cluster, QM/EFP (2% error) does slightly better than the other two schemes (6-7% error). The errors in the 2PA cross sections relative to the full QM calculations can be traced back to even smaller errors in the transition moments and dipole-moment differences between the ground and excited states (Table S2). For thymine clusters, these differences across the three QM/MM schemes are small. The larger excitation energies and smaller dipole-moment differences between the ground and excited states of thymine compared to pNA and PYPb explains why these thymine systems have smaller 2PA cross sections.</p><p>The results for the thymine dihydrates and trihydrate computed with the three QM/MM schemes are also comparable, with errors ranging between 11-20% relative to the full QM calculation. This indicates that polarization has negligible effect on the 2PA cross sections in thymine clusters. Similar to our observation for the pNA clusters, inclusion of explicit water molecules in the QM subsystem for thymine reduces the errors of the QM/EFP calculation. In terms of the structure, the T12 dihydrate can be viewed as superposition of T1 and T2 monohydrate clusters. Considering the large QM/EFP errors in the 2PA cross section for the T1 cluster, the inclusion of the T1 fragment (the tilde indicates that T1 fragment and T1 cluster are structurally similar) in the QM subsystem reduces the error in the 2PA cross sections from 14% to about 3% for the T12 cluster. Similarly, the structure of T112 trihydrate can be described as a superposition of the T11 dihydrate and the T2 monohydrate. The QM/EFP errors in the 2PA cross section relative to the full QM calculation are reduced from 19% to 3% when the T11 fragment (structurally similar to the T11 cluster) of T112 is included in QM.</p><!><p>The lowest excited state of the phenolate chromophore of PYP (PYPb) has a large 2PA cross section (Table III). In the microsolvated clusters, 2PA cross section for this transition increases significantly (up to 110% in the PYPb-W P W P 1 cluster). We also note that the variations in the excitation energy for this transition are large across the three clusters. In the PYPb-W P 1 and PYPb-W P W P 1 clusters, the water molecules are hydrogen-bonded with the phenolate oxygen. This interaction stabilizes the hole NTO of the transition (see SI Fig. S1) more than the particle NTO and results in an increase in the excitation energy. In the PYPb-W C W P 1 cluster, water molecules make hydrogen bonds with both the phenolate and the carboxylic acid group. This stabilizes both the hole and particle NTOs to a similar extent and the excitation energy does not change significantly in this cluster relative to the bare PYPb chromophore.</p><p>The three QM/MM schemes show small errors (<2%) in the excitation energy relative to the full QM calculation. The QM/EFP errors are slightly smaller than those of the other two non-polarizable QM/MM schemes. This is expected for a charged chromophore for which mutual polarization is more pronounced. Interestingly, QM/EFP overestimates the excitation energies for all clusters, whereas the other two QM/MM schemes underestimate the excitation energies.</p><p>In the PYPb-W P 1 and PYPb-W C W P 1 clusters, the QM/EFP errors in the 2PA cross sections relative to the full QM calculation are smaller than the errors of the other two QM/MM schemes, which again signifies that solvent polarization plays an important role in these clusters. The errors for the PYPb-W P W P 1 cluster are comparable across all three schemes. The errors in the 2PA cross sections can be traced back to smaller errors in the transition moments and electronic dipole differences between the ground and excited states given in Table S3. Interestingly, QM/EFP yields more accurate dipole-moment differences than the other two schemes but performs poorly for the transition moments.</p><p>In the next section, we investigate the effect of bulk solvation on PYPb and show that the accuracy of the 2PA cross sections can be significantly improved by including nearest water molecules in the QM subsystem.</p><p>E. Anionic PYP chromophore in aqueous solution Fig. 2 and SI Table S4 show the excitation energies and microscopic 2PA cross sections with degenerate photons for the 21 MD snapshots of PYPb in solution computed with QM/EFP. We use the following notations for calculations with different QM subsystems:  With just the PYPb chromophore in QM, the average 2PA cross section across the snap-shots shows about 28% increase relative to the bare chromophore (Table III) with a large standard deviation. Except snapshot #15, the 2PA cross sections for all other snapshots are higher than the reference gas-phase value. As expected, the excitation energies across all snapshots increase relative to the bare chromophore. The transition moments and dipolemoment differences between the initial and the final states are given in SI Table S4. We note that the transition moments, in general, decrease while the dipole-moment differences increase, the latter overcompensating the former resulting in higher 2PA cross sections in solution relative to the gas phase.</p><p>The QM/EFP results for the microsolvated clusters show that 2PA cross sections are sensitive to the choice of the QM subsystem and that the inclusion of solvent molecules in QM causes noticeable change in the cross section. We carried out a similar analysis for the PYPb chromophore in solution, wherein we extend the QM subsystem by including more water molecules in the QM part. First, we included the first solvation shell (defined by the cutoff distance of 2.5 Å) around the phenolate oxygen (denoted by PYPb+). We see that the PYPb+ 2PA cross sections increase significantly (increase in the average is ∼15%) and excitation energies decrease (by up to 0.11 eV) across the snapshots relative to the QM/EFP results with just the chromophore in QM. On the other hand, the inclusion of the first solvation shell of water molecules around the carboxylic acid group (denoted by +PYPb) increases the 2PA cross sections to a smaller extent than the PYPb+ calculations in five of the six snapshots for which we conducted such calculations. In general, the excitation energies increase slightly (up to 0.02 eV) in the +PYPb calculations relative to the QM/EFP calculations with just the chromophore in QM.</p><p>We note that the 2PA cross sections for the +PYPb+ QM/EFP calculations, wherein the QM subsystem includes the first solvation shells around both the phenolate oxygen and the carboxylic acid group, can be approximated by adding the differences between PYPb+ and PYPb results to the differences between +PYPb and PYPb results:</p><p>The increase in the average 2PA cross section across snapshots for the +PYPb+ calculations is ∼22%. These changes in the 2PA cross sections between QM/EFP calculations with different QM subsystems can be traced back to the corresponding changes in the transition moments and dipole-moment differences between the ground and excited states given in SI Table S4. These calculations clearly suggest that the convergence of the 2PA cross sections with the size of QM subsystem depends on the convergence of the excitation energies, transition moments, and dipole-moment differences.</p><p>Given the above results, a crucial question is then How fast do these properties converge with the size of the QM subsystem? In other words, how large the QM solvation shell around the chromophore should be for the converged bulk results? While we could not carry out such convergence studies with EOM-EE-CCSD, we can investigate the convergence of the key electronic properties contributing to the 2PA cross sections using a lower-level method.</p><p>Towards this goal, we investigated the convergence of the excitation energies, transition moments, and dipole-moment differences by computing them with the CIS/EFP scheme for snapshots #13 and #17 (Fig. 3). For this analysis, we gradually extend the QM subsystem by including more and more water molecules in QM using the EFP-modified geometries (as in the case of the microhydrated clusters, this strategy allows us to compare properties of different embedded QM subsystems using identical structures). In all model structures, the chromophore is oriented in the same way, such that the x-component of the 1PA transition moment and dipole-moment difference always dominate. We approximate the corresponding M xx moment according to Eq. ( 34) and then estimate the error in the 2PA cross section in terms of the error in this 2PA transition moment as a function of the QM size. The use of EFP-modified extended QM subsystem does not introduce significant differences in the behavior of the excitation energies, total transition moments, and dipole-moment differences in our analysis. The differences in these properties for these two EFP-modified snapshots across different embedded QM subsystems (given in SI Table S5) are small when compared to the corresponding values for embedded QM subsystems constructed with the geometries taken from the MD snapshots (given in SI Table S6). We expand the QM subsystem by adding waters within a cutoff distance from either the phenolate oxygen or the carboxylic acid group of the chromophore. We vary the cutoff distance from 0 Å (chromophore only) to 5 Å. For snapshots #13 and #17, the maximum number of water molecules included in QM was 54 and 43, respectively (Fig. 3(a)). The cutoff distance of 2.5 Å gives the structure corresponding to the +PYPb+ calculations discussed above. for snapshot #17 is well-converged after the inclusion of 20 water molecules in the QM subsystem. However, the convergence is slower for snapshot #13 for which convergence is reached only after ∼45 molecules are included in QM. Whereas the absolute errors relative to the calculation with the largest QM subsystem are small for snapshot #17 (<2.5%), the corresponding errors are slightly larger for snapshot #13 (up to 6.2%). Including the first solvation shell decreases the absolute errors for snapshot #17 from 0.9% to 0.1%, while the absolute errors increase for snapshot #13 from 4.3% to 5.5%. On the other hand, the change is larger and convergence is slower (after inclusion of ∼45 water molecules in QM) for snapshot #13. This is highlighted in that the absolute errors in +PYPb+ estimates for this 2PA transition moment are 6.2% for snapshot #13 as against 0.3% for snapshot #17. Since the approximate error in the 2PA cross sections is twice the error in the dominant M xx component, for snapshot #13, we estimate that our best EOM-EE-CCSD 2PA cross section (+PYPb+ calculation) is within ∼12.5% from the converged bulk value.</p><p>Fig. S2, Table S5, and the discussion that follows in the SI provide the CIS/EFP error analysis for all 21 snapshots as a function of the cutoff distance. The errors in the average microscopic 2PA cross sections for PYP and +PYP+ calculations are expressed in terms of the errors in the M xx transition moments relative to the converged values (Eq. (S1)).</p><p>These errors are -12.6% and -3.3%, respectively. This clearly shows that the errors in the average 2PA cross sections relative to the converged values drop significantly when the water molecules in the first solvation shell are treated at the QM level. Thus, we estimate that our best bulk-averaged value of EOM-CCSD 2PA cross section (+PYPb+ calculation) is about -3.3% from the converged result.</p><!><p>We extended the EOM-CCSD method of calculating 2PA cross sections to condensedphase calculations, wherein the quantum system is embedded in the polarizable environment represented by EFP fragments. We also presented a metric based on transition 1PDMs that helps to assess the applicability of a particular QM/MM embedding scheme to studying 2PA transitions.</p><p>We analyzed the solvent effects on 2PA cross sections in microhydrated clusters of pNA, thymine, and PYPb chromophore. We benchmarked the EOM-EE-CCSD/EFP results against full EOM-EE-CCSD calculations for these systems. We observe errors of <20% in the EOM-EE-CCSD/EFP 2PA cross sections relative to the full EOM-EE-CCSD treatment when only the chromophore is included in the QM subsystem. The EOM-EE-CCSD/EFP method captures the main trends in the 2PA cross sections as the full QM calculations. By systematic inclusion of explicit solvent molecules in QM, we find that the calculated 2PA cross sections are sensitive to the QM size. The EOM-EE-CCSD/EFP errors are reduced to ∼5% when nearest water molecule(s) are added to the QM subsystem.</p><p>In the EOM-EE-CCSD/EFP method, we turned off the contribution of the polarization effects and analyzed its impact on the 2PA cross sections. We find that polarization has a large impact on the EOM-EE-CCSD/EFP 2PA cross sections of (anionic) PYPb clusters but not for clusters of thymine and pNA. We also compared our results with EOM-EE-CCSD/MM approach, wherein the solvent molecules are represented by CHARMM charges.</p><p>While we do not see significant differences in the results with these three QM/MM schemes for thymine and pNA clusters, the performance of the EOM-EE-CCSD/EFP results is better for the PYPb clusters.</p><p>We also performed NTO analyses of the 2PA transitions in order to understand their orbital character and to obtain insight into the impact of the solvent. The analyses revealed the intramolecular charge-transfer character of these 2PA transitions, which then allowed us to evaluate the 2PA transition moments from three key components: excitation energies, 1PA transition moments, and dipole-moment differences between the ground and excited states.</p><p>The errors in the EOM-EE-CCSD/EFP 2PA cross sections relative to the full EOM-EE-CCSD results were then traced back to the errors in the excitation energies, 1PA transition moments, and dipole-moment differences.</p><p>We have also reported EOM-EE-CCSD/EFP 2PA cross section of the lowest symmetric transition in PYPb in solution. Our analysis shows that inclusion of few water molecules in the QM subsystem reduces the errors in the excitation energies, 1PA transition moments, and dipole-moment differences significantly relative to the converged values with extended QM subsystems. Consequently, the inclusion of water molecules in the QM improves the calculated 2PA cross section significantly.</p>