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

Electronic structure methods enable highly accurate calculations of molecular structures, energetics, and properties; however, the domain of applicability of ab initio methods in terms of the size of a system is 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 a 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 (e.g., polarizable continuum), they are capable of describing specific interactions such as solvent-solute hydrogen bonding and salt bridges.

Numerous hybrid multi-scale methods¹ 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.² This basic QM/MM strategy has been successfully used for describing both the ground- and excited-state energetics and properties of large systems.^3–6 The QM/MM approach can be improved by using polarizable force fields in which the charge distributions in the MM part can respond to the changes in the 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 description can be achieved by using force fields derived from first principles, such as in the effective fragment potential (EFP) approach.^7–11 In the EFP method, 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. 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, and redox potentials in a variety of solvated systems including protein-bound chromophores.^12–19 QM/EFP is similar to the polarizable embedding (PE) approach.^20–23 The description of electrostatics and polarization is nearly identical in QM/EFP and PE; QM/EFP also includes (non-empirical) terms describing dispersion and exchange-repulsion interactions.

Higher-order response properties such as two-photon absorption (2PA) cross sections are particularly sensitive to the environment.²⁴ 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.²⁵ As pointed out in a recent validation study, the inclusion of solvent effects on 2PA cross sections is critically important for unambiguous comparison between theory and experiment.^26,27

Nonlinear 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 the one-photon transition moments expressed as the matrix elements between the initial and the final states only, the quality of the description of nonlinear 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.

The 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²⁸ or charge-transfer (CT) states delocalized over several solvent molecules.²⁹ This issue can be partially remedied by a careful choice of the QM subsystem, for example, by including some solvent molecules. 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 electronic correlation included in the model Hamiltonian.

In the context of 2PA calculations, the EFP method has been employed in combination with time-dependent density functional theory (TDDFT)³⁰ but not with high-level wave-function methods such as the equation-of-motion coupled-cluster singles and doubles method for excitation energies (EOM-EE-CCSD).^31–33 Coupled-cluster methods have been employed within the polarizable embedding scheme in the calculations of 2PA cross sections.²² In this article, we extend our EOM-CCSD implementation of the 2PA cross sections^34–36 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 the 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: (a) EOM-EE-CCSD/(EFP without polarization effects) and (b) EOM-EE-CCSD/MM, wherein the solvent molecules are represented by CHARMM22 point charges^37,38 (we denote these calculations as QM/CHARMM). 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 in the QM subsystem.

In the QM/EFP embedding scheme, the Hamiltonian of the whole QM/MM system is written as

where H^QM and H^MM describe the QM and MM parts, respectively, and H^QM/MM describes the interaction between these two subsystems. In the QM/EFP approach, the MM system is broken into effective fragments and its energy, E^MM, is computed as the sum of many-body electrostatic, polarization, dispersion, and exchange-repulsion interactions between the fragments

In calculations of solvated species, each solvent molecule is treated as an EFP fragment. The EFP method can also be applied to macromolecules via a fragmentation scheme that can deal with breaking covalent bonds.¹⁷ 

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^39–41 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 the induced dipoles on other fragments are included in the calculation of E^pol. A 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 distances between the permanent and induced multipoles can be found elsewhere.^9,10,42–44

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,

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. 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/MM, 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 the polarizable environment to different electron distributions of each state can be computed perturbatively¹² and were shown to improve the quality of computed energy differences.^{13,14,16–19} Here, we omit these corrections because this state-specific treatment cannot be easily incorporated into the response equations solved in 2PA calculations. In the current implementation,¹² the dispersion and exchange-repulsion interactions between the QM and MM systems (E^disp,QM/MM and E^exch,QM/MM) are treated approximately, as interactions between the EFP fragments, and are added to the state energies; thus, 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

and

where Ψ_gr and Ψ_ex are the ground- and excited-state wave functions. Note that in CIS (configuration interaction singles), TDDFT, and EOM-CC calculations, excitation energies are computed directly and not as energy differences from state-specific calculations.

Our choice of the electronic-structure method for the QM system is EOM-EE-CCSD^31,45,46 in which the excited-state wave functions are given by

Here, Φ₀ is the reference determinant and e^T|Φ₀⟩ is the CCSD wave function. The CCSD cluster operator T and the EOM-CCSD operator R comprise single and double excitation operators

The amplitudes of the operator $T^$ satisfy the following equations:

where $H¯=e−THeT$ 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 $Rk$ and the corresponding energies $Ek$ for EOM-CCSD states are obtained by diagonalizing $H¯$ in the space of the reference, singly, and doubly excited determinants

Since $H¯$ is non-Hermitian, its left (⟨Ψ_k| = ⟨Φ₀L_k|) and right (|Ψ_k⟩ = |R_kΦ₀⟩) eigenstates are not Hermitian conjugates, but form a biorthonormal set

An important feature of the EOM-CCSD/EFP scheme is that it preserves the 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. This biorthogonality is important for computing transition properties such as 1PA and 2PA transition moments.

We compute the 2PA transition moments and cross sections for EOM-CCSD wave functions following the approach described in Refs. 34 and 36. Formally, the 2PA transition moments are given by the following sum-over-states expressions:

and

Here, ω₁ and ω₂ are the frequencies of the two absorbed photons satisfying the 2PA resonance condition ω₁ + ω₂ = Ω_fi with Ω_ni = E_n − E_i. To derive expressions for the EOM-CCSD 2PA transition moments, we replace the dipole-moment operator, μ^a, with the similarity-transformed dipole-moment operator, $μ¯a=e−TμaeT$⁠, and the wave functions, ⟨Ψ_k|s and |Ψ_k⟩s, with the EOM-CCSD target states, ⟨Φ₀L_k|s and |R_kΦ₀⟩s, in the above expressions. Then, we recast these expressions into the numerically and formally equivalent form

and

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

and

where Φ_Is denote the Slater determinants from the target-state manifold (e.g., the reference, singly excited, and doubly excited determinants for EOM-EE-CCSD).

The EOM-EE-CCSD/EFP method is implemented in the Q-Chem quantum-chemistry package.^47,48 The procedure for calculating 2PA cross sections involves the following steps:

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.

2. An iterative self-consistent field 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.

3. The CCSD T- and Λ-amplitudes 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. Using these amplitudes, the CCSD and EOM-EE-CCSD state and transition densities are formed for calculating the state and transition dipole moments.

4. In the 2PA step, the excitation energies without the state-specific corrections (due to 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 2PA cross sections is described in Refs. 34 and 36.

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 (1PDMs), γ,

where $p^†$ and $q^$ create and annihilate electrons in molecular orbitals ϕ_p and ϕ_q, respectively. Importantly, γ is related to the experimental observables, e.g., the 1PA transition moment is given by

The elements of transition 1PDMs can be interpreted as the amplitudes of an excitation operator that generates the final state from the initial correlated state

The squared norm of γ,

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.^49,50

Transition 1PDM can be used to define an exciton’s wave function (Ψ_exc), which is a two-particle quantity that contains essential information about the changes in electron density upon the transition^51–54 

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).

Equations (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

As discussed in Subsection II F, this reasoning can be extended to quantify the extent of charge-transfer character.^55,56

Further simplification in exciton’s representation can be achieved by applying singular-value decomposition (SVD) to γ. The SVD procedure reduces the description of electronic excitation to one-electron transitions between pairs of orbitals called natural transition orbitals^54,56–59 (NTOs). Using NTOs, the exciton’s wave function can be written as

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 σ_Ks are significant, giving rise to a 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 many-body wave functions in terms of matrix elements between orbitals, facilitating the connection between molecular orbital theory and experimental observables

A detailed description of the NTO analysis and its implementation in Q-Chem is given in Refs. 54 and 55.

We recently extended the concepts of transition 1PDMs and NTOs to 2PA transitions.³⁶ 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 transition.

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 the 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

where $ϕ̃K,Pa$ and $ϕ̃K,Ha$ are the Kth particle and hole NTOs for the 2PA a-component (a, b = x, y, z) transition 1PDM with a singular value of $σKa$⁠.

Since the NTOs are linear combinations of atomic orbitals (χs), i.e.,

Eq. (26) can be expressed in the basis of χs as

The corresponding 2PA a-component transition 1PDM in the atomic orbitals basis is

By partitioning the QM system into the chromophore (A) and solvent (B), we can partition the above 1PDM into the four parts, $γf←iaA←A$⁠, $γf←iaA←B$⁠, $γf←iaB←A$⁠, and $γf←iaB←B$⁠, where ^C←Dγ represents the contribution of the hole residing on fragment D and a particle residing on fragment C

The probability of finding the hole on fragment D and the particle on fragment C is now given by

The four probabilities corresponding to the four fragment 1PDMs give a measure of the charge-transfer contributions within and across fragments into the 2PA transition moment. This is similar to the charge-transfer (CT) numbers used to quantify the extent of charge transfer in electronic transitions.^54,55 Breaking the system into a QM solute and an 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

In other words, the performance of fragment approaches such as QM/EFP deteriorates with increasing contributions of the solvent to the 2PA transition 1PDMs.

We computed these CT numbers for the 2PA transitions in microhydrated clusters studied here as well as for the representative snapshots from solvated PYP simulations. The analysis reveals that the charge-transfer probabilities between the chromophore and the solvent are low (⁠$Qf←iaA←B+Qf←iaB←A$ < 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 a noticeable charge-transfer character in the virtual states.

All calculations were carried out using the Q-Chem quantum-chemistry package.^47,48 We used the B3LYP/6-31+G^* optimized geometry for pNA, the RI-MP2/cc-pVTZ geometry for thymine, and the CCSD/6-31+G^* geometry for PYPb in the calculations of bare chromophores. The structures used in the QM/EFP calculations for the microsolvated clusters of pNA, thymine, and PYPb were taken from Refs. 12, 60, and 61, respectively; they are shown in Fig. 1. 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 (1) EFP fragments, (2) EFP fragments without polarization, and (3) CHARMM22 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⁻¹ (3.1492 eV). We used the TheoDORE package⁵³ for computing the charge-transfer probabilities following the 2PA NTO calculations with Q-Chem.

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 structures 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 that 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 supplementary material.

To model bulk solvation, we performed molecular dynamics (MD) simulations with a PYPb chromophore and one sodium cation in a periodic box of 1717 water molecules using the NAMD software.⁶² We used the CHARMM22 parameters for water and the sodium cation; the parameters for PYPb are given in the supplementary material. 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 the supplementary material). Along 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.

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.

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 solute-solvent interactions by analyzing the orbital character of the transitions. We calculated the 2PA NTOs using our wave-function analysis toolkit.³⁶ The leading NTO pairs and the molecular orientations are shown in Fig. 2.

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. 2) 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. A more detailed analysis of this transition can be found in Ref. 36. The charge-transfer character is further supported by analyzing the participation ratios (PR_NTO) computed for the z-component 1PDM and its component ωDMs (see Ref. 36 for definitions). Whereas the PR_NTO for the z-component 1PDM is ∼1-2 for both degenerate and non-degenerate cases, the PR_NTOs for its component ωDMs are >20. This is a signature of intramolecular charge-transfer transitions. Intramolecular charge-transfer 2PA transitions are driven by the change in the dipole moment upon excitation (⟨Ψ_f|μ^a|Ψ_f⟩ − ⟨Ψ_i|μ^a|Ψ_i⟩), as given by an approximate expression for the 2PA transition moment

which reduces to

for degenerate photons and diagonal elements of the 2PA transition-moment tensor.^36,63,64 Equation (34) 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 (Ω_fi = ω₁ + ω₂). Thus, understanding the difference in the permanent dipole moments of the ground and excited states, transition dipole moment, and excitation energy in the microsolvated clusters is the key to understanding how solvent impacts 2PA cross sections. 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.

In our calculations, thymine (C_s) is in the xy plane. The lowest A′ transition with degenerate photons has a small microscopic cross section dominated by the M_xx and M_yy components. 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 PR_NTOs 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 the dipole moment of the molecule in the excited state, which suggests red shift in aqueous solution.

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 x-component 2PA 1PDM gives two dominant NTO pairs, which represent the 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 the supplementary material (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.⁶⁵ By 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 PR_NTOs 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. Equation (34), therefore, explains the large 2PA cross section for this transition.

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. This is expected for weakly interacting solute-solvent systems. The change in the oscillator strength is also small. Table S1 in the supplementary material 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 increases the 2PA cross section by up to 8% for degenerate photons and 15% for non-degenerate photons in accordance with Eq. (33). Across the three clusters, the full QM cross section increases from dihydrates to tetrahydrates and decreases slightly from tetrahydrates to hexahydrates.

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 (Table I). Omitting the polarization contributions of the EFP fragments 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 the 2PA transition with non-degenerate photons for the hexahydrate than the dihydrate. Similarly, the QM/CHARMM approach 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 the tetrahydrate and an increase from tetrahydrate to hexahydrate. Yet, the small differences between QM/EFP, QM/EFP with polarization turned off, and QM/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.

We also report the QM/EFP results for the 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 was reported for the core-ionization energies of solvated glycine.¹⁹ 

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 supplementary material (Table S1). We note that the errors in the excitation energies for these three QM/MM approaches relative to the 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).

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%, the oscillator strength by up to 14%, and the 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–T2 relative to the bare thymine (see Table II), but in the T3 cluster, the excitation energy is higher than in bare thymine. The difference in trends can be attributed to the following structural feature: only in the T3 cluster, the water molecule is near the carbonyl group that is across from the methyl group, thereby stabilizing the hole orbital shown in Fig. 2(b). 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.

The results for the thymine dihydrates and trihydrates computed with the three QM/MM schemes are also comparable, with errors ranging between 11% and 20% relative to the full QM calculation. This indicates that polarization has a 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 polyhydrates reduces the errors of the QM/EFP calculation. In terms of the structure, the T12 dihydrate can be viewed as a 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 $T̃$1 fragment (the tilde indicates that the $T̃$1 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 the 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 $T̃$11 fragment (structurally similar to the T11 cluster) of T112 is included in QM.

The lowest excited state of the phenolate chromophore of PYP (PYPb) has a large 2PA cross section (Table III). In the microsolvated clusters, the 2PA cross section for this transition increases significantly (up to 110% in the PYPb-W_PW_P1 cluster). We also note that the variations in the excitation energy for this transition are large across the three clusters. In the PYPb-W_P1 and PYPb-W_PW_P1 clusters, the water molecules are hydrogen-bonded with the phenolate oxygen. This interaction stabilizes the hole NTO of the transition (see Fig. S1 of the supplementary material) more than the particle NTO and results in an increase in the excitation energy. In the PYPb-W_CW_P1 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.

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.

In the PYPb-W_P1 and PYPb-W_CW_P1 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_PW_P1 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.

In Sec. IV E, 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.

Table S4 in the supplementary material and Fig. 3 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: (a) PYPb means that only PYPb is included in QM, (b) PYPb+ means that PYPb plus all water molecules within 2.5 Å of the phenolate oxygen are included, (c) +PYPb includes PYPb plus all water molecules within 2.5 Å of the carboxylic acid group, (d) +PYPb+ includes PYPb plus all water molecules within 2.5 Å of the phenolate oxygen and the carboxylic acid group, (e) PYPb++ includes PYPb plus all water molecules within 3.5 Å of the phenolate oxygen, and (f) PYPb+++ includes PYPb plus all water molecules within 4.0 Å of the phenolate oxygen.

With just the PYPb chromophore in QM, the average 2PA cross section across the snapshots shows about 28% increase relative to the bare chromophore (Table III) with a large standard deviation. With an exception of snapshot #15, the 2PA cross sections 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 dipole-moment differences between the initial and the final states are given in Table S4 (supplementary material). 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.

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 a 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.

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 the PYPb+ and PYPb results to the differences between the +PYPb and PYPb results

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 Table S4 (supplementary material).

Since the changes due to the inclusion of water molecules in QM at the two ends of the PYPb chromophore are approximately additive, we test if the inclusion of the first solvation shell in QM gives sufficiently converged results by including water molecules beyond the first solvation shell only at the phenolate end. For snapshots #13 and #17, we include all water molecules within 4.0 Å from the phenolate oxygen. For snapshot #13, the excitation energy and 2PA cross section for the PYPb+++ calculations show less than a 1% change relative to the results for the PYPb+ calculation with just the first solvation shell in QM. The change in the transition moments from PYPb+ to PYPb+++ for this snapshot is positive, which cancels the negative change in the dipole-moment difference leading to an overall smaller change in the 2PA cross section. For snapshot #17, the change in the 2PA cross section from PYPb+ to PYPb++ (water molecules within 3.5 Å) is about 5%. However, the 2PA cross section change by <1% from PYPb++ to PYPb+++ calculations. The transition moments and dipole-moment differences also show negligible changes from PYPb++ to PYPb+++. 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.

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. 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 all snapshots (the results for snapshots #13 and #17 are presented in Fig. 4). 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 the 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 dominates. 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 the 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 Table S5 of the supplementary material) are small when compared to the corresponding values for embedded QM subsystems constructed with the geometries taken from the MD snapshots (given in Table S6 of the supplementary material for snapshots #13 and #17).

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. 4(a)]. The cutoff distance of 2.5 Å gives the structure corresponding to the +PYPb+ calculations discussed above.

Figure 4(b) shows the behavior of the excitation energy as a function of the QM water molecules. In both snapshots, we observe an initial steep drop from the calculation with the bare chromophore in QM to the +PYPb+ calculation. The excitation-energy curve for snapshot #17 plateaus after about 20 water molecules (∼4 Å cutoff) are included in QM. By contrast, this curve for snapshot #13 shows a slow gradual decrease and plateaus when about 42 water molecules (∼4.6 Å cutoff) are included in the QM. Including the first solvation shell (+PYPb+ calculations) around the two ends of PYPb reduces the error in the excitation energy from 3.5% to 2.6% and 1.6% to 0.8% for snapshots #13 and #17, respectively, relative to the calculations with the largest QM subsystems.

Figure 4(c) shows the behavior of the x-component transition moment as a function of QM water molecules. After an initial steep increase upon addition of the first solvation shell, the changes are much smaller for both snapshots. The errors in the +PYPb+ calculations are 1.8% and 0.4% for snapshots #13 and #17, respectively, relative to the calculations with the largest QM subsystems.

Figure 4(d) shows the behavior of the x-component dipole-moment difference as a function of QM water molecules. As for the excitation energies, the dipole-moment difference 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%.

Figure 4(e) shows the approximate M_xx transition moment, calculated using Eq. (34), as a function of the number of QM water molecules. Similar to the behavior of the excitation energy and dipole-moment difference, the xx-component 2PA transition moment for snapshot #17 changes by no more than 3.3% after the inclusion of ∼20 water molecules. 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.

Figure S2, Table S5, and the discussion in the supplementary material 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)]. 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 the EOM-CCSD 2PA cross section (+PYPb+ calculation) is about 3-4% below the converged result. The analysis of the convergence of 2PA cross sections with respect to the size of the QM subsystem shows that the first solvation shell provides a good balance between computational cost and converged results for the final averaged spectra, despite slower convergence for individual snapshots. A similar conclusion for another polarizable embedding QM/MM scheme has been reported in Ref. 66 for converged X-ray absorption spectra averaged over multiple snapshots.

We extended the EOM-CCSD method of calculating 2PA cross sections to condensed-phase 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 us to assess the applicability of a particular QM/MM embedding scheme for studying 2PA transitions.

We analyzed the solvent effects on 2PA cross sections in microhydrated clusters of pNA, thymine, and the PYPb chromophore. We benchmarked the EOM-EE-CCSD/EFP results against full EOM-EE-CCSD calculations for these systems. When only the chromophore is included in the QM subsystem, the errors in the EOM-EE-CCSD/EFP 2PA cross sections relative to the full EOM-EE-CCSD treatment are <20%. The EOM-EE-CCSD/EFP method captures the main trends in the 2PA cross sections as the full QM calculations. A systematic inclusion of explicit solvent molecules in QM reveals that the calculated 2PA cross sections are sensitive to the QM size. When nearest water molecule(s) are added to the QM subsystem, the EOM-EE-CCSD/EFP errors drop down to ∼5%.

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 the EOM-EE-CCSD/MM approach, wherein the solvent molecules are represented by CHARMM charges. While we do not see significant differences in the results with these three QM/MM schemes for thymine and pNA clusters, EOM-EE-CCSD/EFP performs better for the PYPb clusters.

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. 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.

We have also reported the 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.

See supplementary material for (1) 1PA NTOs for the studied transitions in pNA, thymine, and PYPb in the gas phase; (2) key quantities such as excitation energies, transition moments, dipole-moment differences, and 2PA cross sections for the microhydrated clusters and aqueous PYPb calculated with EOM-EE-CCSD/EFP and CIS/EFP methods; (3) force-field parameters for PYPb; (4) Cartesian coordinates; and (5) a separate .txt file containing pdbs for the 21 MD snapshots of PYPb in solution.

A.I.K. acknowledges support by the U.S. National Science Foundation (Grant No. CHE-1566428). K.D.N. acknowledges helpful discussions with Dr. Ilya Kaliman on EFP and Dr. Atanu Acharya and Tirthendu Sen on NAMD simulations.