Impact of solvation on the electronic resonances in uracil

Abstract

Interactions of low-energy electrons with the DNA and RNA nucleobases are known to form metastable states, known as electronic resonances. In this work, we study electron attachment to solvated uracil, an RNA nucleobase, using the orbital stabilization method at the Equation of Motion-Coupled Cluster for Electron Affinities with Singles and Doubles (EOM-EA-CCSD) level of theory with the Effective Fragment Potential (EFP) solvation method. We benchmarked the approach using multireference methods, as well as by comparing EFP and full quantum calculations. The impact of solvation on the first one particle (1p) shape resonance, formed by electron attachment to the π* LUMO orbital, as well as the first two particle one hole (2p1h) resonance, formed by electron attachment to neutral uracil's π–π* excited state, was investigated. We used molecular dynamics simulations for solvent configurations and applied charge stabilization technique-based biased sampling to procure configurations adequate to cover the entire range of the electron attachment energy distribution. The electron attachment energy in solution is found to be distributed over a wide range of energies, between 4.6 eV to 6.8 eV for the 2p1h resonance, and between −0.1 eV to 2 eV for the 1p resonance. The solvent effects were similar for the two resonances, indicating that the exact electron density of the state is not as important as the solvent configurations. Multireference calculations extended the findings showing that solvation effects are similar for the lowest four resonances, further indicating that the specific solute electron density is not as important, but rather the water configurations play the most important role in solvation effects. Finally, by comparing bulk solvation to clusters of uracil with a few water molecules around it, we find that the impact of microsolvation is very different from that of bulk solvation.

---

1 Introduction

DNA damage by ionizing radiation is caused to a large extent by indirect effects.¹ For example, interaction of ionizing radiation with the cellular environment produces highly reactive secondary species that can attack DNA. These reactive molecules can lead to lesions due to base damage and/or backbone breaking.^2–5 While DNA repair enzymes typically fix this damage without any harm to the cell,^6,7 improper repair may sometimes lead to cell death or carcinogenic events. Radiation therapy, a common cancer treatment, targets tumors with radiation, effectively shrinking them by inducing cell death.⁸ Among the reactive secondary species, low energy electrons, with kinetic energy below 20 eV, have been of special interest more recently,^9–11 since Sanche and coworkers¹² demonstrated, through their experiment on plasmid DNA, that low energy electrons, even those with sub-ionization energy, can lead to single and double strand breaks in DNA. The strand breaks in DNA have been attributed to the formation and reactivity of transient negative ions (TNIs).²

TNIs, formed by electron attachment to neutral species with negative electron affinity, are electronically metastable, and form electronic resonances. The formation of TNIs is the first step in the DNA strand breaks via dissociative electron attachment mechanism.^12–19 The second step involves electron transfer to an anti-bonding orbital, leading to the cleavage of the corresponding bond.¹⁶ For instance, electron transfer to the anti-bonding orbital of the glycosidic C–N bond leads to base release, while transfer to the anti-bonding orbital of the sugar–phosphate C–O bond results in cleavage of the backbone leading to DNA strand breaks.²⁰ Much work has been focused on electrons with energy below 3 eV, which are known to cleave the sugar–phosphate C–O bond.²¹ Only one-particle (1p) shape resonances are formed at this low energy range by attachment to a LUMO orbital, and there is fairly good understanding of the mechanism leading to strand breaks via the formation of 1p shape resonances of DNA nucleobases in the gas phase.^22–24

Strand breaks are even more prominent for higher energy electrons.^25,26 At these energies, the involvement of core-excited or two particle-one hole (2p1h) resonances is expected. These resonances are generated when the energy of the electrons is higher than the energy needed to excite the neutral molecule. So, the electron scattering induces an electronic excitation, and the electron is attached to the excited electronic configuration. While it is expected that 2p1h resonances are primarily important for the strand breaks above 4 eV, limited information exists about them, because of the difficulty to study them theoretically or experimentally. 2p1h resonances are challenging to simulate since rearrangement of more than one electron is needed. Experimentally, resonances reported in the energy range of 4–6 eV and 9–11 eV have been interpreted as 2p1h (core-excited shape or Feshbach) resonances.¹⁶

The electronic resonances of uracil have been extensively studied both theoretically and experimentally in the gas phase, since uracil is a prototype for pyrimidine bases and the simplest nucleobase.^27–38 Gianturco and Lucchese reported three π* resonances between 2.27 eV and 6.5 eV of energy with two σ* resonances located at 0.012 eV and 10.37 eV based on their scattering calculations.³² Later, Dora et al.³³ identified three shape resonances below 4.95 eV and three Feshbach resonances in the energy range of 6.17 eV to 8.12 eV using the R-matrix method. Cheng and Chen³⁴ performed Density Functional Theory (DFT) based stabilization method combined with analytic continuation to obtain resonance parameters, suggested the presence of three π* shape resonances below 5 eV and two σ* resonances in the energy region 5–8 eV. Kossoski and coworkers also observed three π* type shape resonances in uracil.³⁵ Similar results were found by Ehara et al. using CAP/SAC-CI (Complex Absorbing Potential/Symmetry-Adapted Cluster-Configuration Interaction) method.³⁶ Later, Fennimore and Matsika^37,38 reported three π* type shape resonances and two 2p1h type resonances. They characterized the resonance at energy around 5 eV as the first 2p1h resonance, formed by electron attachment to neutral uracil's π–π* triplet excited state which is approximately 4 eV higher than the ground state of neutral uracil. Another resonance around 6.5 eV of energy was attributed to electron attachment to an n–π* excited state of the neutral.

Since water is the most abundant molecule in the cellular environment, it is important to understand how aqueous solvent affects these electronic resonances. Few studies on the electronic resonances of DNA bases have been performed in microsolvated environment^39–46 compared to the gas phase. In the gas phase, DNA bases form dipole bound anions in the ground state with adiabatic electron affinities close to zero.^27,47–49 However, upon microhydration, the valence bound state, which appears as resonance state in the gas phase, becomes adiabatically bound, and even a single water molecule is enough for this transformation.^48,50,51 The geometry of the valence bound anion differs from that of the neutral species. Smyth and coworkers,⁴² based on their theoretical study on microsolvated uracil and thymine, showed that microsolvation leads to increased lifetime and reduced energy of the resonance state. This has also been seen for cytosine.⁵² In contrast, Kočišek and coworkers⁴³ based on their experimental results, found that microhydration suppresses the dissociative channel for the fragmentation of the N–H bond, conflicting with Smyth's theoretical prediction. Later, Sieradzka and coworkers⁴⁴ used R-matrix method to study resonances in thymine (methylated uracil) solvated with five water molecules. They concluded that hydrogen bonding affects the stability of resonance states based on donor and acceptor ability of water molecules. Varella and coworkers⁴⁵ performed Monte Carlo simulations to obtain uracil water clusters with six solvent molecules based on minimal distribution function. Their study revealed that microsolvation generally shifts the π* shape resonance states to lower energies. In a recent study, Verlet and coworkers⁴⁶ examined clusters of uracil with different number of water molecules using two-dimensional photoelectron spectroscopy and identified the two lowest resonances in microsolvated environment, suggesting a linear variation with the number of waters.

A few recent studies have applied molecular dynamics simulations^53,54 and Monte Carlo simulations⁴⁵ to go beyond microsolvation and better account for aqueous phase statistics. Mukherjee et al.⁵⁴ used hybrid quantum mechanics/molecular mechanics (QM/MM) simulations to study the electron attachment to solvated nucleobases. They concluded that the initial electron attachment is localized on water and acts as a doorway for the formation of the final nucleobase bound state. This electron transfer is facilitated by mixing of electronic and nuclear degrees of freedom. They reported that the ground anionic state is stabilized by the presence of water and becomes bound. Likewise, Anstöter et al.⁵³ using DFT based ab initio molecular dynamics simulations and orbital stabilization at the Equation of Motion-Coupled Cluster for Electron Affinities with Singles and Doubles (EOM-EA-CCSD) combined with effective fragment potentials (EFP), concluded that while the initial electron attachment is not localized on uracil, within 15 fs, the excess electron localizes on the uracil nuclear framework on the valence π* state. This π* state becomes stable within 1 ps. They found that both the solvent effect (solvent reorganization) and geometry relaxation (uracil core) lead to stabilization of the valence bound state.

To the best of our knowledge, there have been no theoretical studies to date to account for the effect of bulk solvation on the 2p1h resonances, despite their importance in DNA strand breaks.²⁵ The main focus of this work is to investigate the distribution of electron attachment energies of the first 2p1h resonance when solvation is taken into account. Our calculations also provide the electron attachment energies for the 1p resonance, and we compare the solvation effects on these two resonances. We also compared the solvent effect in clusters to investigate the effect of solvation in bulk. This work examines only vertical attachment, so we do not incorporate the stabilization from nuclear rearrangement or solvent reorganization.

2 Theory and computational details

2.1 Classical molecular dynamics

In this work, snapshots of previous classical molecular dynamics (CMD)⁵⁴ performed by some of us are used. A summary of how these were performed is below. The CMD simulations were performed to obtain the configurations that mimic the bulk aqueous environment around uracil. Uracil is frozen in the geometry obtained at the MP2/6-31G* level of theory.^55,56 For the water molecules, the TIP3P model⁵⁷ was used, while uracil was treated with the CHARMM compatible force field⁵⁸ using the NAMD15 package.⁵⁹ In the minimization step, uracil was kept fixed in a cubic box of 40 Å edge with 2460 water molecules. The minimized structure was then heated to 300 K using Langevin thermostat. The equilibration run was carried out for 500 ps using periodic boundary conditions. During the production run, constant temperature at 300 K and constant pressure (1 bar) conditions were maintained with Nose–Hoover Langevin piston pressure control. We have studied 300 different configurations that generated in equal time interval of 25 picosecond from a 10 ns production run trajectory.

2.2 Electronic structure methods

The EOM-EA-CCSD method is used to obtain the electron attachment energies.^60,61 This method can produce accurate electron attachment (EA) energies using an operator to attach an electron to the neutral reference, thus balancing neutral and anion correlation. In the most common cases, the reference is the neutral closed shell ground state of a molecule. This however fails to properly treat anionic states that differ by more than one electron from the ground state, which is the case for 2p1h resonances. To overcome this deficiency, one can choose a different reference. If we choose as reference the neutral triplet ground state, then we can generate 2p1h resonances which have that state as their parent state. This is the approach we have chosen here. There are however some serious limitations with this approach. It can be used only if the ground triplet state is the correct parent state for the resonance we are interested in, so in most cases it is limited to being able to produce only one 2p1h resonance. Furthermore, we have found in previous work³⁸ that there is mixing between the 1p and 2p1h resonances if they have similar energies. This mixing cannot be described when the 1p and 2p1h resonances are generated from different calculations. Nevertheless, we have used this approach here since it is enabling us to run calculations for 300 snapshots to account for solvation effects in a statistical manner. We will compare with multireference methods to examine its validity.

In order to check the results against multireference methods, we have carried out complete active space self-consistent field (CASSCF)⁶² calculations for representative solvated configurations. An active space of 11 electrons in 9 orbitals was used, CAS(11,9). Since the focus of this work is to consider only resonances initiated from π orbitals, we included all the valence π(a′′) orbitals (5π and 3π*). In addition, we included an extra diffuse π character orbital in order to have avoided crossings in the stabilization curves to obtain the widths. The orbitals included in the active space are shown in ESI,† Fig. S1. 5 doublet states of A′′ symmetry were averaged in the CASSCF calculations. The energy of neutral uracil was obtained with a CASSCF(10,9) with the same active space and an average of 5 singlet A′ states. The energy of the neutral at α = 1 is used as reference to obtain the EAs. The (11,9) active space was also used in QM/MM calculations using classical point charges to describe the solvent (TIP3P point charges of −0.834 and 0.417 atomic units are used for oxygen and hydrogen, respectively). The MOLRPO^63–65 software was used for the CASSCF calculations.⁶⁶

The solvation effects in EOM-EA-CCSD have been considered using the EFP method.^67–70 EFP is a quantum mechanical potential that can either be viewed as a fragmentation model or as a polarizable force field model with the parameters obtained from quantum calculations. Like force fields, the total non-covalent interaction energy can be written as the sum of several interaction terms, but unlike force fields, EFP uses fragments rather than atoms as the basic units. This method allows an accurate representation of the interactions between the QM region and the solvent fragments, and the solvent-solvent fragments. To account for solvent response to electron rearrangement in the EOM target states, a perturbative non-iterative correction is computed for each EOM root using the one-electron density of that state which is used to re-polarize the environment.^71,72 The structures used for the electronic structure calculations are taken from the CMD trajectory. EOM-EA-CCSD/EFP single point calculations were carried out using the Q-Chem computational software (version 5.4).⁷³

For one configuration, the solvent effect is also evaluated using polarizable continuum model (PCM)⁷⁴ as well as QM/MM with water as point charges. PCM is performed using a non-equilibrium conductor like-PCM (C-PCM) model as implemented in Q-Chem,⁷³ with dielectric constant of 78.39 for the solvent. In the QM/MM method, point charges of −0.834 and 0.417 atomic units are used for oxygen and hydrogen, respectively, in water (TIP3P charges). In these approaches, the solvent effect is incorporated at the Hartree–Fock level, and EA-EOM-CCSD/cc-pVDZ+1p is used to obtain the EA.

2.3 Charge stabilization method

Using the charge stabilization (CS) method,^75,76 the electron attachment energy of the resonance states can be obtained at affordable cost because we can apply it using a smaller number of ab initio points and smaller basis set than other more demanding methods for electronic resonances. As will be discussed in detail in the results, the qualitative trends are reproduced, since the error is consistent for all configurations. In this technique, a small positive charge (z) is added at the center of the molecule. This charge is increased in small steps to generate the potential that efficiently converts a metastable state into a bound anionic state. Then the energy of the state of interest is calculated as a function of the positive charge. The obtained energies vary linearly with the extra positive charge because of the linear dependence of the Coulomb interactions, and extrapolation for z = 0 provides the electron attachment energy for the desired state. A representative plot for one snapshot is presented in Fig. 1.

Fig. 1 Electron attachment energies plotted against charge in charge stabilization method. Extrapolation to zero provides the resonance position.

In this study, the electron attachment energies to build the CS plots for the first 1p and 2p1h states are obtained at the combined EOM-EA-CCSD/EFP level of theory with the cc-pVDZ basis set.⁷⁷ For all 300 snapshots from the CMD trajectory, the positive charge is varied from 0.6 to 0.8 with a difference of 0.1. This optimum charge range was tested by examining the linear behavior of the energy for a wide range of charges on 25 snapshots. Results for both 3 and 5 points are shown in ESI† (Table S3 and Fig. S5). For lower values of charge, the resonance does not become bound. The CS method allows for the calculation of electron attachment energies of the resonance states, but it does not provide information about the resonance widths. Since this is a much cheaper approach, we can apply it to many configurations and then choose only appropriate configurations for the more expensive orbital stabilization approach, which will provide more accurate positions, as well as widths.

2.4 Orbital stabilization method

Hazi and Taylor first introduced the stabilization method,⁷⁸ which utilizes stabilization graphs plotted using the excited energies of an anion. In the orbital stabilization method (OSM), the spatial extent of the Gaussian basis functions is varied using a scaling parameter (α). Resonances can be identified as stationary states exhibiting avoided crossings with discretized continuum states. The real energies at the avoided crossings are analytically continued into the complex plane using the Generalized Padè Approximant (GPA) approach. The quadratic polynomial GPA is implemented in this work, where the following equation is used

E²P + EQ + R = 0. (1)

P, Q, and R are the coefficients which are polynomials of the scaling parameter (α) and are given as follows:
The GPAs are denoted by (n_i, n_j, n_k) where n_i, n_j, and n_k are the number of coefficients in P, Q, and R polynomials. In this work, (4, 4, 4), (5, 5, 5) and (7, 7, 7) GPAs are used to calculate resonance parameters. The (7, 7, 7) GPA is usually reported unless this did not converge in which case either (4, 4, 4) or (5, 5, 5) GPA was used. These cases are listed in Table S4 of the ESI.† The stabilization plot is generated at the QM/EFP level for the solvated uracil where EOM-EA-CCSD/cc-pVDZ+1p is used for the QM region. This basis set is generated by adding an extra 1p diffuse function to the cc-pVDZ basis set, the exponent of which is half of the most diffuse p function of the standard basis set. For isolated uracil, a stabilization plot is also generated with a larger basis set (aug-cc-pVDZ+1p)⁷⁷ at the EOM-EA-CCSD level. Here also the extra p function is added in a similar fashion. The energies and α values from the avoided crossing are substituted in the GPA equation providing linear equations in terms of unknown coefficients. The unknown coefficients are obtained by solving linear equations using standard matrix techniques. These known coefficients are substituted back in the initial GPA and the complex stationary point corresponding to the resonance state is found by dE/dα = 0. After substituting the computed complex stationary point in the polynomial, complex energies are obtained.^37,79

Since the OSM uses a larger basis set than CS and many more points are needed, it is not possible to apply it to all 300 snapshots. For this reason, we used the distribution obtained from the CS results on the 300 snapshots to obtain representative snapshots that span the whole distribution, and we carried out OSM on those representative snapshots. Fig. 2 shows the 13 snapshots that were chosen. Five of them were chosen to have an average value of EA (chosen from the peak of the distribution), and the others were chosen from either side of the distribution taken at steps of 0.5 standard deviations. The figure only shows the structure of uracil and water molecules in the first solvation shell.

Fig. 2 Structure of 13 snapshots used for GPA calculations from the CS distribution. The figure only shows the structure of uracil and water molecules in the first solvation shell.

3 Results and discussion

3.1 Description of 1p and 2p1h resonances in uracil

Uracil is known to form 1p and 2p1h resonances upon electron attachment.³⁷ The 1p resonances are formed by an electron attachment to unoccupied molecular orbitals. There are three low lying π* resonances, and the lowest one is formed by attachment of an electron to the LUMO, 1π* orbital (see Fig. 3). Alternatively, when the electron has enough energy to excite the molecule first, inelastic scattering can lead to an excited state, followed by attachment of an electron, leading to the formation of 2p1h or core-excited resonances. When the energy of the anionic state is below that of the neutral excited state, a core-excited Feshbach resonance is formed, whereas when it is above the energy of the neutral excited state, a core-excited shape resonance is formed. While the term core-excited resonances has been most often used in the literature, here we prefer the term 2p1h to describe these resonances, since it is more descriptive and avoids confusion as core electrons are not involved. Since 2p1h states involve electron attachment to neutral excited states, these resonances are expected to have energies similar to the neutral excited states. Singlet or triplet states can be the parent state for a resonance, but usually triplet states have lower energies. In the case of neutral uracil, the lowest triplet state is a ³(π)¹(π*)¹ around 3.5–4 eV, while a ³(n)¹(π*)¹ state is about 1 eV higher.⁸⁰ The corresponding 2p1h resonances can be formed by attachment to the anti-bonding π* orbital, leading to (π)¹(π*)² or (n)¹(π*)². Both of these 2p1h resonances have been studied in the gas phase,³⁷ and it was found that the (π)¹(π*)² is lower in energy, and since it has A′′ symmetry it can mix with the nearby 1p resonances that have the same symmetry. For this reason, we focus on the (π)¹(π*)² 2p1h resonance in this work.

Fig. 3 Energy level diagram defining all the states studied here. S₀ and T₁ represent the ground singlet and triplet states of neutral uracil, respectively. D₀(1p) and D_N(2p1h) represent the first 1p and 2p1h resonances, respectively. The energy gap between the singlet and triplet state is denoted by ΔE(S–T). EA-D₀ is the electron attachment energy for the 1p resonance. EA-D_N(T) and EA-D_N(S) represent the electron attachment energies to the triplet or singlet neutral state, respectively, to produce the 2p1h resonance. The corresponding orbital level diagrams for each resonance are shown on the right.

2p1h resonances present a theoretical challenge for the underlying electronic structure method, when trying to use extended quantum mechanical methods to study them. Multireference methods are important in order to treat properly the mixing, however they can be too expensive and complicated to use with explicit solvent. Alternatively, EOM-EA-CCSD has been used successfully for 1p resonances, but it is less straightforward to apply it to 2p1h resonances. We can however start from the triplet ground state of the neutral and attach an electron. This can only be used if we can obtain the appropriate target triplet state. In uracil the ground triplet state is the (π)¹(π*)¹, so we can apply this technique for the resonance of interest. Starting from the triplet state as reference we can also obtain the first shape resonance by attachment of an electron to the π orbital instead of the π*. In that case however, we cannot obtain the width since the shape resonance is stable compared to the triplet neutral state. Fig. 3 shows the energy level diagram for the electron attachment energies of the 1p and 2p1h resonances studied in this work, and the relevant neutral states. The first state S₀, represents the ground state of neutral uracil. Above it lies the first 1p resonance state D₀(1p), and the energy gap between the two states is defined as the electron attachment energy for shape resonance (EA-D₀). The 2p1h resonance (D_N(2p1h)) is obtained by electron attachment to the triplet T₁ state, and this EA is denoted as EA-D_N(T). In order to obtain the EA from the ground neutral state we add the gap between the singlet and triplet neutral states, ΔE(S–T), and we obtain EA-D_N(S). Here we will use EOM-EA-CCSD, while we will examine the effect of multireference methods in Section 3.7.

3.2 Benchmark

In this section, we examine the effect of our methodological choices.

3.2.1 Basis set effect. The basis set effect is evaluated by treating gas phase uracil with the chosen basis set (cc-pVDZ+1p) and comparing the results to the larger basis set (aug-cc-pVDZ+1p). It can be seen from Table 1 that increase in basis set size, results in decreased electron attachment energy and width by 0.16 eV and 0.12 eV, respectively. For the bulk solvated uracil at the QM/EFP level, we could not use the larger basis set because the system was not converged to the correct triplet state at the CCSD level, so we could not obtain the desired 2p1h resonance. Instead the calculations converged to a Rydberg state. Therefore, in this work, except for isolated uracil, all the calculations to generate orbital stabilization graphs are performed with the cc-pVDZ+1p basis set. A much more detailed study of the effect of basis sets on the 1p resonances of cytosine was performed by Verma et al.⁵² using aug-cc-pVDZ, aug-cc-pVTZ and aug-cc-pVQZ basis sets expanded with extra p functions. In that study it was found that the additional polarization is important and stabilizes the resonances by 0.2–1 eV. In the current work, these basis sets are prohibitive because of the electronic structure methods used. So it is clear that there will be an important effect by the basis set choice, and the exact positions and widths we calculate have some uncertainties. The solvation effects and trends that we observe however should still be valid.

Table 1 Electron attachment energies from the triplet state, EA-D_N(T), and corresponding widths (in parenthesis) of the first 2p1h resonance of gas phase uracil and uracil at different level of solvation obtained with the GPA method where the stabilization curves are generated at the EOM-EA-CCSD/cc-pVDZ+1p level of theory for one representative snapshot (snapshot 8) (see Fig. 4 for structures). All values are in eV

System	EA-DN(T)(Γ)
Uracil (QM)/aug-cc-pVDZ+1p	1.611 (0.046)
Uracil (QM)/cc-pVDZ+1p	1.778 (0.167)
Uracil (QM) + PCM	2.139 (0.029)
Uracil (QM) + 2460H2O (point charges)	1.885 (0.027)
Uracil (QM) + 2460H2O (EFP)	1.759 (0.027)
Uracil (QM) + 3H2O (QM) + 2457H2O (EFP)	1.749 (0.053)
Uracil (QM) + 7H2O (QM) + 2453H2O (EFP)	1.696 (0.068)
Uracil (QM) + 1aH2O (QM)	1.519 (0.035)
Uracil (QM) + 1bH2O (QM)	1.649 (0.014)
Uracil (QM) + 2H2O (QM)	1.279 (0.026)
Uracil (QM) + 3H2O (QM)	1.612 (0.017)
Uracil (QM) + 4H2O (QM)	1.366 (0.011)
Uracil (QM) + 7H2O (QM)	1.424 (0.066)

---

Fig. 4 Clusters with various numbers of water for a representative snapshot (snapshot 8). This snapshot is used for results in Table 1. (a) and (b) are monohydrated structures of uracil where water is hydrogen bonded to the same atom of uracil with different orientation in space. (c) uracil-water cluster with two water molecules. (d) uracil-water cluster with three water molecules. (e) uracil-water cluster with four water molecules and (f) uracil-water cluster with seven water molecules.

3.2.2 Description of the solvent. Anions are very sensitive to solvation, so it is important to use the appropriate method for the description of solvent. A common way is to treat solvent molecules with EFP which is computationally affordable for large number of solvent molecules.^67–70 In previous work, we had examined the accuracy of EFP for excited states and for resonances in anions, and we found that it is quite accurate as long as the electron density is mostly localized on the solute rather than the water.⁸¹

Here we perform some additional benchmarking. For this purpose, we have considered uracil with three water molecules (uracil–(H₂O)₃) cluster. In order to be able to do full QM calculations for several snapshots, we restricted the number of waters to three. We chose five snapshots that have similar EA, at the average value (see snapshots 5–9 in Fig. 2). Table 2 presents the electron attachment energy from the neutral triplet state of uracil–(H₂O)₃ clusters for the first 2p1h shape resonance. To generate uracil water clusters, a minimal distance between a water molecule and uracil is used as a cutoff in the snapshots from the MD simulations. The cutoff is varied to include only three water molecules around uracil. The efficiency of EFP to account for the solvent effect is evaluated against a full QM approach. For three out of the five clusters under study, the electron attachment energy at the EFP level is overestimated, while for the other two clusters it is underestimated, with respect to water treated at the QM level. However, the error when using EFP is less than 0.06 eV, for both the energy and the width, with an average error of about 0.025 eV. So, EFP can be considered as a good alternative to treat solvent molecules in this case.

Table 2 Electron attachment energy from the triplet state, EA-D_N(T), and width (in parenthesis) for the first 2p1h resonance of uracil water clusters with 3 water molecules, uracil–(H₂O)₃, where the water molecules are treated at the QM level or by using EFP. Results from orbital stabilization method at the EOM-EA-CCSD/cc-pVDZ+1p level are shown. The differences between the two approaches are also shown as errors. All values are in eV

Snapshot	QM	EFP	Error	Error
	EA-DN(T)(Γ)	EA-DN(T)(Γ)	EA-DN(T)	Γ
5	1.656 (0.017)	1.716 (0.013)	0.060	−0.004
6	1.260 (0.014)	1.242 (0.017)	−0.018	0.003
7	1.443 (0.025)	1.419 (0.088)	−0.024	0.063
8	1.612 (0.017)	1.638 (0.044)	0.026	0.027
9	1.348 (0.024)	1.353 (0.072)	0.005	0.048

---

Other approaches for the solvent are compared in Table 1, which shows results of using a polarizable continuum model or point charges in a traditional QM/MM scheme or EFP. The comparison is made for a single snapshot using all 2460 water molecules for EFP and point charges. While the widths using the three approaches are very similar, the positions vary more. In particular, the EA value from PCM differs from EFP by 0.38 eV, while the point charges differ from EFP by 0.13 eV. More importantly, PCM predicts that the resonance is destabilized in solution compared to gas phase uracil by 0.36 eV while EFP predicts a very small stabilization of 0.02 eV. PCM at the simple level incorporated here is clearly incapable of accurately predicting the solvation effects.

We also examine whether the results change if we include some explicit QM water molecules in addition to the bulk EFP water molecules. Adding three QM water molecules to the EFP treatment changes the EA by only 0.01 eV, but it has a stronger effect for the width, which doubles (although in actual values the change is only 0.026 eV). Adding seven water in the quantum region changes the EA by another 0.05 eV and increases the width by another 0.015 eV. The results do not seem to converge by systematically adding water molecules, although we are not able to go beyond seven QM water molecules. In addition, there is no clear correlation between the EA and the number of hydrogen bonds between uracil and water. In general, we can conclude that explicit quantum water molecules have some effect, especially for the widths. However, given the substantial increase in computational cost, we can still expect semi-quantitative results with EFP without adding water molecules at the quantum level.

3.3 Charge stabilization results

Fig. 5 shows the distribution of electron attachment energies using 300 snapshots for both the 1p (Fig. 5a) and 2p1h (Fig. 5b) resonances. The EA for the 2p1h resonance is calculated using the triplet state as reference, so the obtained energy has to be adjusted by adding the energy gap between the ground singlet state and the triplet state ΔE(S–T). The distribution for ΔE(S–T) is also shown in Fig. 5c. The 1p distribution is centered around 1.5 eV and has a spread of about 2 eV, while the 2p1h distribution is centered around 6.5 eV and has a similar spread. While these results will be improved using OSM, we already see a qualitative picture. Similarly to what was found in the gas phase, the 2p1h resonance is not a Feshbach resonance, being energetically above the neutral triplet state on average by almost 2.5 eV according to these calculations. The most apparent observation is that while the distributions for the resonances are wide, with a range of 2 eV, the distribution for the neutral state ΔE(S–T) is much narrower, with an energy range of only 0.26 eV, starting from 3.84 eV to 4.10 eV, with an average of 3.97 eV. These differences in the distribution show very clearly how much more important the solvent effect is for the anion compared to the neutral.

Fig. 5 Electron attachment energy calculated using the charge stabilization method at the EOM-EA-CCSD/cc-pVDZ/EFP level, for (a) the first 1p resonance, EA-D₀, and (b) the first 2p1h resonance, EA-D_N(S). (c) Singlet–triplet gap, ΔE(S–T), distribution at the EOM-EA-CCSD/cc-pVDZ/EFP level.

3.4 Orbital stabilization results

While charge stabilization gives an estimate of the positions, it provides no information about the widths. So, we move to OSM in order to refine the positions and obtain more information about the widths. However, due to the high computational cost, we cannot perform OSM for all 300 configurations. Therefore, we have chosen 13 configurations from the entire distribution (Fig. 5b) at different positions from the extreme left to the extreme right as described in methodology and shown in Fig. 2.

Fig. 6 shows a representative orbital stabilization plot, indicating the points around the avoided crossing that were used in the GPA. Other plots are shown in ESI† (Fig. S6–S18). The (7, 7, 7) GPA is used in all cases, except when it did not work, either (4, 4, 4) or (5, 5, 5) GPA was used. We have included information about which GPA was used in each case in ESI† (Table S4). Since the reference is a triplet state, the lowest state is obtained by attaching an electron to the HOMO, which gives rise to the first 1p resonance. This has negative energy in the plot since it is below the neutral triplet state. We use these energies to obtain the position of the 1p resonance, but we do not get information about the widths. The 2p1h resonance is at positive energies, indicating that it is above the neutral triplet state, as was also seen in the CS results. Here, we can use the avoided crossings to obtain the widths. We have done this for 13 snapshots, and the results are shown in Table 3.

Fig. 6 Orbital stabilization plot for a representative snapshot using EOM-EA-CCSD/cc-pVDZ+1p/EFP. Blue points indicate the points that were used in GPA.

Table 3 Electron attachment energies, EA-D_N(T), and widths Γ in parenthesis, for different configurations chosen from the distribution in Fig. 5b picking them using the average and standard deviation. All values are calculated using EOM-EA-CCSD/cc-pVDZ+1p and shown in eV. The number of water molecules in the cluster calculations are also shown. Configurations are shown in Fig. 2

Snapshot	# Water	Clusters	Bulk
		EA-DN(T)(Γ)	EA-DN(T)(Γ)
1 (Average − 2.0 SD)	10	1.122 (0.029)	0.912 (0.030)
2 (Average − 1.5 SD)	12	1.838 (0.050)	1.148 (0.025)
3 (Average − 1.0 SD)	10	1.132 (0.048)	1.354 (0.021)
4 (Average − 0.5 SD)	9	0.917 (0.030)	1.511 (0.014)
5 (Average)	8	1.353 (0.102)	1.709 (0.038)
6 (Average)	9	1.807 (0.128)	1.716 (0.034)
7 (Average)	10	1.199 (0.073)	1.749 (0.038)
8 (Average)	7	1.539 (0.044)	1.759 (0.027)
9 (Average)	10	1.091 (0.033)	1.778 (0.025)
10 (Average + 0.5 SD)	10	1.653 (0.058)	1.981 (0.045)
11 (Average + 1.0 SD)	10	1.716 (0.082)	2.107 (0.018)
12 (Average + 1.5 SD)	10	2.156 (0.039)	2.356 (0.065)
13 (Average + 2.0 SD)	6	1.488 (0.069)	2.569 (0.012)

---

Similarly to the CS results, the EA-D_N(T) obtained from OSM follow the same trends, with snapshot 1 having the lowest EA and snapshot 13 the highest. The correlation is actually very strong if we plot the OSM vs. CS results. We have plotted results for both the 1p and 2p1h resonances in Fig. 7, and in both cases there is a very strong visual correlation. In order to examine the correlation more quantitatively, we carried out a linear regression of the electron attachment energies at the GPA versus charge stabilization level for the 13 snapshots. In both cases, the coefficient of determination, R², is very close to 1, an indication of the very strong correlation: R² = 0.9898 for the 1p resonance and R² = 0.9969 for the 2p1h resonance. This made us realize that we can use the CS results to extrapolate and obtain OSM (GPA) results for all 300 snapshots. To obtain the GPA electron attachment energies for all 300 snapshots, we use the expressions obtained from the linear fitting. The following equations show the expressions for 1p (eqn (2)) and 2p1h (eqn (3)) resonances.

E^1p_GPA = 0.9942(E^1p_CS) − 0.4909 (2)

E^2p1h_GPA = 0.9880(E^2p1h_CS) − 0.8103 (3)

Fig. 7 Regression correlating the GPA and CS energies for the 1p (top) and 2p1h (bottom) resonances obtained from EOM-EA-CCSD.

E _GPA denotes the EA obtained from the GPA calculations while E_CS denotes the EA from CS calculations. From the above two expressions, the obtained GPA EA energies distribution for all the 300 snapshots of the 1p and 2p1h resonances are presented in Fig. 8. The distributions are wide, similarly to the CS ones. The EA values have shifted to lower energies by about 0.5 eV for 1p resonance and around 0.8 eV for 2p1h resonance. The 1p resonance state is distributed over a wide range (of 2 eV) with average at 0.99 eV. The electron gets attached at energies ranging from −0.09 to 1.99 eV forming the lowest energy valence π* resonance, or in some cases at the lower energy tail a bound state is predicted. However, configurations with negative value are limited and more than 98% of the total configurations possess positive value, and are metastable.

Fig. 8 Electron attachment energies for (a) the first 1p resonance, EA-D₀, and (b) the first 2p1h resonance, EA-D_N(S), calculated using GPA method at the EOM-EA-CCSD/cc-pVDZ+1p/EFP level.

We want to put these values into context by comparing this to previous work on the stability of the first resonance in uracil. While here we see mostly metastable character, the state is eventually stabilized by about 3 eV when the solvent and the geometry are allowed to relax to the anionic electronic density.⁵³ So even though vertical electron attachment leads to a metastable state, after relaxation the uracil anion becomes stable in aqueous environment.

The 2p1h resonance electron attachment energies are also distributed over a wide range of 2.2 eV with an average at 5.7 eV. So, the electrons that can form 2p1h resonances can have energies in between 4.6 eV to 6.8 eV. In this case all EA are above the triplet state, so there is no chance for this to be Feshbach in character. Based on the Gaussian distribution for 300 snapshots, it can be concluded that the electron attachment energies are significantly affected by the distribution of water molecules in both resonances. The widths for the 13 snapshots are shown in Table 3. Their values range from 0.012–0.065 eV, without any obvious correlation to the EAs. This may be a reflection of the uncertainties that we have in the widths, because of the small size of the basis set.

Finally, we want to examine whether the solvation effects leading to the extended distributions are similar or different for the two types of resonances we have calculated. For this purpose, we plotted the EA values with respect to each other in Fig. 9. The figure shows a surprisingly strong linear correlation between the EA for 1p vs. EA for 2p1h resonance, with an R² = 0.9933. This implies that the solvation effects are similar for the two resonances, and whichever structure stabilizes the 1p resonance it will also stabilize the 2p1h resonance. This is a very interesting result, and can be potentially very useful in such studies.

Fig. 9 (a) Correlation plot between EA-D_N(S) and EA-D₀ at the EOM-EA-CCSD level. Linear regression gives: y = 1.0106x + 4.7116 with R² = 0.9848. (b) Correlation plots between EA of the three resonances 2²A′′, 3²A′′, 4²A′′ vs. 1²A′′ (EA-D₀) at the CASSCF level. Linear regression gives: 2²A′′: y = 1.3245x + 0.6306 with R² = 0.9948; 3²A′′: y = 1.2473x + 3.9387 with R² = 0.9915; 4²A′′: y = 1.2401x + 5.2333 with R² = 0.991. All values are obtained from the orbital stabilization method.

In order to check that this correlation is not artificial since we have used the same reference to do EOM-EA-CCSD for both resonances, we also calculated the 1p resonance using the neutral ground state as the reference for three snapshots with very different solvation effects. The results are shown in ESI† (Table S6) which shows that the EA-D₀ is similar regardless of the reference used. Thus, it is confirmed that this is not an artifact of the way we used to calculate the solvation effects.

3.5 Comparison of microsolvated and bulk solvated uracil

An important question we want to understand better is what causes the wide distribution, and specifically, whether local effects, such as hydrogen bonding, or effects of bulk solvation are responsible. In a broader context, it is also important to know whether we can extrapolate to bulk solvation by studying clusters. In order to address these questions, we compared clusters to bulk aqueous environment. This is achieved by performing the OSM calculations for bulk and microsolvated uracil for the 13 snapshots shown in Fig. 2. For bulk solvation, uracil is solvated with 2460 water molecules, while in microsolvated uracil, water is considered up to the first solvation shell (using a cutoff radius of 2.7 Å). This cutoff radius is chosen based on the work of Krylov and coworkers, where by using the first minima of radial distribution functions it is ensured that all water molecules that have a direct H-bond to uracil are included.⁸² So, the number of water molecules is not the same in all 13 snapshots, and it is listed in Table 3. These snapshots are arranged in increasing value of electron attachment energies obtained from the charge stabilization method. Table 3 and Fig. 10 show the comparison between the microsolvation and bulk solvation for these 13 snapshots.

Fig. 10 Comparison of EA-D_N(T) (top) and corresponding widths (bottom) between a cluster including only water molecules around the first solvation of uracil and bulk solvation for the 2p1h resonance, calculated using EOM-EA-CCSD/cc-pCDZ+1p. The thirteen snapshots used are tabulated in Table 3. Similar results for the 1p resonance are shown in ESI† (Fig. S19 and Table S5).

Fig. 10 makes the comparison very clear. Since the snapshots are ordered in order of increasing EA, the values for bulk increase from left to right. On the contrary, there is no similar increase for the EA obtained from the clusters. There is no correlation between the bulk-solvation values and the cluster values. In fact, the same configurations that are significantly stabilized by solvent in bulk-solvation compared to isolated uracil are destabilized when water only up to the first solvation shell is considered and vice versa. For example, the electron attachment energy in the case of the second configuration is destabilized in cluster, while it is significantly stabilized in bulk solvation. Similarly, the last snapshot which is destabilized upon bulk-solvation is stabilized in cluster and the difference in electron attachment energy between the two is more than 1 eV. The widths behave differently as well, with bulk-solvation widths distributed over a narrow range of 0.05 eV while for the clusters they distributed over a range of around 0.1 eV. The widths are also not correlated between the bulk and cluster values, although in this case there is no clear trend in either case. Finally, the 1p resonance shows similar effects, i.e. a difference between the bulk and microsolvation effects. Plots equivalent to Fig. 10 for the 1p resonance are shown in ESI† (Fig. S19 and Table S5).

From this comparison it can be concluded that water molecules beyond the first solvation shell that are not directly hydrogen bonded to uracil play an important role in the solvation effects of the resonances. Both the solvent–solvent interactions and long range solute–solvent interactions are crucial when considering the aqueous environment around uracil anion. Clusters with few water molecules can be misleading for drawing conclusions about the effect of bulk aqueous environment. Quite interestingly however, even though the microsolvation effects are different from bulk solvation, they also show a correlation between the 1p and 2p1h effects, similar to what we saw for bulk solvation in Fig. 9. An equivalent plot is shown in ESI† (Fig. S20) for the 13 microsolvated structures, exhibiting similar correlation.

Another comparison is made using a single snapshot and varying the number of water molecules included. Snapshot 8 is chosen because it is a representative snapshot close to the average of the CS distribution, as shown in Fig. 2. This way we assume our results will be more representative of the average. Fig. 4 shows clusters with various water molecules for this snapshot, while Table 1 presents the electron attachment energies and widths. Very interestingly, the EA-D_N(T) changes depending on the number of water molecules included in the cluster, and does not seem to converge to a specific value. The value for the cluster with three water molecules is 1.61 eV and it decreases to 1.37 eV for four water molecules, and somewhat increases again for seven water molecules. Overall, the 2p1h state, similar to bulk solvation, is also stabilized in uracil water cluster by the solvent effect but the extent to which the state is stabilized is different from that of bulk solvation. It should be noted that these different clusters show different hydrogen bonding patterns. An attempt was made to see if the hydrogen bonding pattern has a clear correlation to the EA, but we were not able to find one.

3.6 Comparison to isolated uracil

Table 1 can be used to examine the effect of solvation compared to isolated uracil on snapshot 8, which is the representative snapshot with average EA values. Study on the snapshot 8 at different level of solvation, reveals that in bulk solvation the energy is slightly decreased compared to isolated uracil while the widths are reduced by one order of magnitude, which leads to longer lifetime of the 2p1h state in bulk solvation. Microsolvation also shows similar effect on the width, while the energy is stabilized by a greater extent compared to bulk-solvation. To expand on these observations, we examined the effect of solvation compared to isolated uracil for four other snapshots close to the mean of the CS distribution (snapshots 5–9 in Fig. 2). The discussion here then only applies to the average.

Fig. 11 presents the solvent effect on the 1p and 2p1h resonances compared to isolated uracil for these five snapshots. The figure shows results for bulk solvation described only using EFP (black curves), bulk solvation described with EFP plus three water molecules in the QM region (red curves), as well as clusters with three water molecules only (blue curves). It is immediately apparent that the EA results for bulk solvation are similar qualitatively regardless of the description, but the cluster results are very different in all cases. In bulk, the energy is stabilized by a small amount compared to free uracil for both 1p and 2p1h resonances. The stabilization is less than 0.2 eV for the 1p resonance and less than 0.1 eV for the 2p1h resonance. The cluster effects for the energies have the same direction but larger magnitude in general for both resonances, and the stabilization can be up to 0.6 eV. The width of the 2p1h resonance decreases about 0.1–0.15 eV and there isn't a great difference (or trend) between the bulk and clusters. It is emphasized again that these effects on stabilization are only for snapshots very near the average of the distribution.

Fig. 11 Difference in electron affinity values for (solvated minus isolated) uracil for five different snapshots at three different level of calculations for (a) the first 1p resonance, (b) the first 2p1h resonance. (c) Difference in width for (solvated minus isolated) uracil for five different snapshots at three different level of calculations for the first 2p1h resonance. The actual values are shown in Table S7 (ESI†). Results at the EOM-EA-CCSD/cc-pVSZ+1p level.

3.7 Multireference description

A main issue with our calculations is that by using EOM-EA-CCSD from the triplet state to describe the 2p1h resonance, we are not able to describe its mixing with the third 1p shape resonance which exists in uracil at similar energies. This mixing has been identified in previous gas phase calculations.^37,38 Multireference methods are able to describe the mixing because they treat both the 1p and the 2p1h resonances equivalently in the same calculation. Performing OSM using multireference approaches for all the snapshots we did in this work, however, is prohibitive. So, as a test we have examined the performance of CASSCF on five representative snapshots. We chose the snapshots to represent the distribution, choosing snapshots 1, 3, 8, 11, and 13 from Fig. 2. Point charges were used for the water molecules. The results are compared to EOM-EA-CCSD in Table 4. CASSCF predicts all the four lowest A′′ resonances, three 1p resonances and the 2p1h resonance. However, the third and fourth resonances mix heavily. Using EOM-EA-CCSD we only calculated the first and fourth resonances.

Table 4 EA energies (and widths in parenthesis) in eV for the four A′′ resonances in solvated uracil using MM. All calculations used the cc-pVDZ+1p basis set. In CASSCF calculations 5 states were averaged

Snapshots	Method	1^2A′′	2^2A′′	3^2A′′	4^2A′′
		(1π*)^1	(2π*)^1	(3π*)^1	(π)^1(π*)^2
1 (Average − 2.0 SD)	EOM-EA-CCSD	0.097			4.815(0.021)
3 (Average − 1.0 SD)	EOM-EA-CCSD	0.484			5.283(0.021)
8 (Average)	EOM-EA-CCSD	0.942			5.659(0.027)
11 (Average + 1.0 SD)	EOM-EA-CCSD	1.331			6.139(0.018)
13 (Average + 2.0 SD)	EOM-EA-CCSD	1.729			6.486(0.012)
1 (Average − 2.0 SD)	CASSCF	1.751(0)	2.998(0.001)	6.140(0.089)	7.413(0.523)
3 (Average − 1.0 SD)	CASSCF	1.722(0)	2.907(0)	6.132(0.011)	7.311(0.366)
8 (Average)	CASSCF	2.373(0)	3.749(0)	6.802(0.113)	8.198(0.362)
11 (Average + 1.0 SD)	CASSCF	2.527(0)	3.903(0)	7.059(0.061)	8.463(0.462)
13 (Average + 2.0 SD)	CASSCF	3.065(0)	4.746(0)	7.827(0.067)	8.966(0.562)

---

The CASSCF energies show the same general trends as EOM-EA-CCSD, i.e. the same snapshots that are stabilized with EOM-EA-CCSD are stabilized at the CASSCF level as well. In order to demonstrate the trends quantitatively, we show the correlation plots between the EA of the 1p (1²A′′) resonance and the 2p1h (4²A′′) resonances in Fig. 9b. The CASSCF results show a strong linear correlation as well, similarly to the EOM-EA-CCSD. Using these 5 snapshots R² = 0.991 for CASSCF in comparison to R² = 0.985 for EOM-EA-CCSD in Fig. 9a. Furthermore, the linear equations relating the 1p and 2p1h resonances are

y = 1.2401x + 5.2333 (4)

for the CASSCF EA and

y = 1.0106x + 4.7116 (5)

for the EOM-EA-CCSD values. These equations indicate that the gap between the two resonances is mostly a constant at the EOM-EA-CCSD level, i.e. the 2p1h energies are 4.71 eV above the 1π* resonance regardless of solvation effects. For CASSCF the slope is not 1, so the gap is not always constant, but it is close to 5.23 eV.

CASSCF has the advantage that it can produce all the resonances between the 1p and 2p1h. So we have four resonances to examine how they are affected by solvation. The correlation plots in Fig. 9b include also the correlation between 1²A′′ vs. 2²A′′ and between 1²A′′ vs. 3²A′′. Quite interestingly, the correlation is strong in all cases. This is another strong indication that the solvation effects are not very sensitive to the actual electronic structure wavefunction of each state.

Comparison between the widths obtained at the CASSCF level with the widths obtained with EOM-EA-CCSD (Table 4) shows that they are affected more by mixing. The width of the fourth resonance is much larger at the CASSCF level compared to EOM-EA-CCSD. This may be an effect of the mixing with the 3π* character. However, it should be pointed out that the uncertainty of the values of widths at the CASSCF level is higher since it also depends on the number of diffuse orbitals in the active space which are needed to have continuum like states in the stabilization plots. Here we only have one diffuse orbital so there is a limited number of continuum like states. A detailed discussion on the challenges of CASSCF with OSM is given in ESI,† Section S2. In any case, it is reasonable to expect that the width will be affected by the wavefunction mixing.

The CASSCF results show similar trends to the EOM-EA-CCSD. Taken together with the correlation observed in Fig. 9 they suggest that solvation effects on the resonance position and width of nucleobases do not depend as much on the electronic distribution of the solute, but they depend strongly on the water configurations. So, even though the EOM-EA-CCSD calculations do not show mixing between the 1p and 2p1h resonances, the solvation effects predicted should still be valid.

A disadvantage of CASSCF is that it cannot accurately describe the imbalance between the neutral and the anion. As a result, the first 1p resonance is more than 1 eV higher than the corresponding value at the EOM-EA-CCSD level. If we focus on the gap between the first and fourth resonances, EOM-EA-CCSD predicts it to be around 4.7 eV, while CASSCF 5.2 eV. According to our previous work, dynamical correlation with addition of perturbation theory is needed to bring both values lower.³⁸ Comparing the values in Table 4, however, we can establish further whether this disadvantage is important. The EA values obtained at the EOM-EA-CCSD level were correlated to the CASSCF values. This correlation plot is shown in the ESI† (Fig. S4) and shows that there is indeed correlation and one could use one theory to predict the other. This nice correlation also suggests that point charges, which were used with CASSCF, can predict the same qualitative effects as EFP, which was used with EOM-EA-CCSD.

4 Conclusions

In this study we have examined electron attachment to solvated uracil, focusing on the first 2p1h resonance. Using snapshots from a molecular dynamics simulation, and the EOM-EA-CCSD method with EFP description for the water molecules combined with OSM, we were able to obtain a distribution of the electron attachment energies and the widths of the resonance. The distribution for the electron attachment energy is broad, ranging between 4.6 eV to 6.8 eV. The distribution for the first 1p resonance is also obtained, and it is found to be similarly wide, with energies in the range of −0.1 eV to 2.0 eV, so some of the configurations lead to a stable anion.

The 1p and 2p1h resonances studied show very similar solvation effects. A wide distribution of energies is predicted for both, and the solvation effects for the two resonances correlate very strongly, suggesting that the effect is driven mostly by the solvent distribution, while the solute electron density plays a minor role.

Using multireference CASSCF we were able to calculate the lowest four resonances for a small number of snapshots, which are representative of the whole distribution. The CASSCF results confirmed the EOM-EA-CCSD observation that solvation effects are not driven by the solute electron density, and they extended them even further to include all resonances. As a result, even though there is mixing between the third and fourth resonances, which is properly described at the CASSCF level but is missing at the EOM-EA-CCSD level, the solvation effects on EA are not affected by the mixing. The widths are affected by the mixing.

Comparisons between clusters of uracil with a small number of water molecules and uracil in bulk solvation show that microsolvation is very different from bulk solvation, leading to very different stabilization effects. So, we should not rely on microsolvation to extrapolate to bulk solvation in this case.

Since several compromises have been made in this work in order to be able to do calculations for 300 snapshots, there are some uncertainties in our calculated values. Nevertheless, we are still able to make significant observations on trends as summarized above.

Data availability

The data supporting this article have been included as part of the ESI.†

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under Award No. DE-SC0019394, as part of the Computational Chemical Sciences Program. Part of the computational work was performed using the Expanse at San Diego Supercomputer Center (SDSC) through allocation CHE140114 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program, which is supported by National Science Foundation grant numbers 2138259, 2138286, 2138307, 2137603, and 2138296.⁸³ AKD acknowledges support from the CRG project of SERB, India.

Notes and references

1. B. D. Michael and P. O'Neill, Science, 2000, 287, 1603–1604.
2. E. Alizadeh, T. M. Orlando and L. Sanche, Annu. Rev. Phys. Chem., 2015, 66, 379–398.
3. D. S. Frohlinde, in Comparison of Mechanisms for DNA Strand Break Formation by the Direct and Indirect effect of Radiation, ed. M. G. Simic, L. Grossman, A. C. Upton and D. S. Bergtold, Springer US, Boston, MA, 1986, pp. 19–27.
4. R. Téoule, Int. J. Radiat. Biol. Relat. Stud. Phys., Chem. Med., 1987, 51, 573–589.
5. P. Jeggo and M. Löbrich, Radiat. Prot. Dosim., 2007, 122, 124–127.
6. T. C. Evans and N. M. Nichols, Curr. Protoc. Mol. Biol., 2008, 84, 3–9.
7. M. R. Lieber, Annu. Rev. Biochem., 2010, 79, 181–211.
8. M. Lomax, L. Folkes and P. O'Neill, Clin. Oncol., 2013, 25, 578–585.
9. E. Alizadeh, A. G. Sanz, G. García and L. Sanche, J. Phys. Chem. Lett., 2013, 4, 820–825.
10. L. Sanche, Eur. Phys. J. D, 2005, 35, 367–390.
11. L. Sanche, in Low Energy Electron Damage To DNA, ed. M. K. Shukla and J. Leszczynski, Springer Netherlands, Dordrecht, 2008, pp. 531–575.
12. B. Boudaiffa, P. Cloutier, D. Hunting, M. A. Huels and L. Sanche, Science, 2000, 287, 1658–1660.
13. A. M. Stanislav, A. Pshenichnyuk and A. S. Komolov, Int. Rev. Phys. Chem., 2018, 37, 125–170.
14. S. Antonello and F. Maran, Chem. Soc. Rev., 2005, 34, 418–428.
15. A. Kumar and M. D. Sevilla, Handb. Comput. Chem., 2015, 1–63.
16. X. Luo, Y. Zheng and L. Sanche, J. Chem. Phys., 2014, 140, 155101.
17. I. Bald, R. Čurík, J. Kopyra and M. Tarana, in Dissociative Electron Attachment to Biomolecules, ed. A. V. Solov'yov, Springer International Publishing, Cham, 2017, pp. 159–207.
18. Z. Li, Y. Zheng, P. Cloutier, L. Sanche and J. R. Wagner, J. Am. Chem. Soc., 2008, 130, 5612–5613.
19. J. Simons, Acc. Chem. Res., 2006, 39, 772–779.
20. R. Barrios, P. Skurski and J. Simons, J. Phys. Chem. B, 2002, 106, 7991–7994.
21. J. Berdys, I. Anusiewicz, P. Skurski and J. Simons, J. Am. Chem. Soc., 2004, 126, 6441–6447.
22. F. Martin, P. D. Burrow, Z. Cai, P. Cloutier, D. Hunting and L. Sanche, Phys. Rev. Lett., 2004, 93, 068101.
23. B. Kumari, A. Huwaidi, G. Robert, P. Cloutier, A. D. Bass, L. Sanche and J. R. Wagner, J. Phys. Chem. B, 2022, 126, 5175–5184.
24. S. Tonzani and C. H. Greene, J. Chem. Phys., 2006, 124, 054312.
25. X. Luo, Y. Zheng and L. Sanche, J. Chem. Phys., 2014, 140, 155101.
26. R. Schürmann, T. Tsering, K. Tanzer, S. Denifl, S. V. K. Kumar and I. Bald, Angew. Chem., Int. Ed., 2017, 56, 10952–10955.
27. N. A. Oyler and L. Adamowicz, J. Phys. Chem., 1993, 97, 11122–11123.
28. S. Denifl, S. Ptasińska, G. Hanel, B. Gstir, M. Probst, P. Scheier and T. Märk, J. Chem. Phys., 2004, 120, 6557–6565.
29. T. Takayanagi, T. Asakura and H. Motegi, J. Phys. Chem. A, 2009, 113, 4795–4801.
30. R. Abouaf and H. Dunet, Eur. Phys. J. D, 2005, 35, 405–410.
31. Y. Kawarai, T. Weber, Y. Azuma, C. Winstead, V. McKoy, A. Belkacem and D. Slaughter, J. Phys. Chem. Lett., 2014, 5, 3854–3858.
32. F. A. Gianturco and R. Lucchese, J. Chem. Phys., 2004, 120, 7446–7455.
33. A. Dora, J. Tennyson, L. Bryjko and T. Van Mourik, J. Chem. Phys., 2009, 130, 164307.
34. H.-Y. Cheng and C.-W. Chen, J. Phys. Chem. A, 2011, 115, 10113–10121.
35. F. Kossoski, M. Bettega and M. D. N. Varella, J. Chem. Phys., 2014, 140, 024317.
36. Y. Kanazawa, M. Ehara and T. Sommerfeld, J. Phys. Chem. A, 2016, 120, 1545–1553.
37. M. A. Fennimore and S. Matsika, Phys. Chem. Chem. Phys., 2016, 18, 30536–30545.
38. M. A. Fennimore and S. Matsika, J. Phys. Chem. A, 2018, 122, 4048–4057.
39. S. Kim and H. F. Schaefer, J. Chem. Phys., 2006, 125, 144305.
40. S. M. Bachrach and M. W. Dzierlenga, J. Phys. Chem. A, 2011, 115, 5674–5683.
41. X. Bao, H. Sun, N.-B. Wong and J. Gu, J. Phys. Chem. B, 2006, 110, 5865–5874.
42. M. Smyth, J. Kohanoff and I. Fabrikant, J. Chem. Phys., 2014, 140, 184313.
43. J. Kocisek, A. Pysanenko, M. Fárnk and J. Fedor, J. Phys. Chem. Lett., 2016, 7, 3401–3405.
44. A. Sieradzka and J. D. Gorfinkiel, J. Chem. Phys., 2017, 147, 034303.
45. L. Cornetta, K. Coutinho and M. D. N. Varella, J. Chem. Phys., 2020, 152, 084301.
46. G. A. Cooper, C. J. Clarke and J. R. Verlet, J. Am. Chem. Soc., 2022, 145, 1319–1326.
47. C. Desfrançois, V. Periquet, Y. Bouteiller and J. Schermann, J. Phys. Chem. A, 1998, 102, 1274–1278.
48. J. Hendricks, S. Lyapustina, H. De Clercq and K. Bowen, J. Chem. Phys., 1998, 108, 8–11.
49. D. Tripathi and A. K. Dutta, Int. J. Quantum Chem., 2019, 119, e25875.
50. H. Motegi and T. Takayanagi, THEOCHEM, 2009, 907, 85–92.
51. I. Anusiewicz, P. Skurski and J. Simons, J. Phys. Chem. A, 2020, 124, 2064–2076.
52. P. Verma, M. Mukherjee, D. Bhattacharya, I. Haritan and A. K. Dutta, J. Chem. Phys., 2023, 159, 214303.
53. C. S. Anstoter, M. DelloStritto, M. L. Klein and S. Matsika, J. Phys. Chem. A, 2021, 125, 6995–7003.
54. M. Mukherjee, D. Tripathi and A. K. Dutta, J. Chem. Phys., 2020, 153, 044305.
55. M. Head-Gordon, J. A. Pople and M. J. Frisch, Chem. Phys. Lett., 1988, 153, 503–506.
56. P. C. Hariharan and J. A. Pople, Theor. Chim. Acta, 1973, 28, 213–222.
57. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, J. Chem. Phys., 1983, 79, 926–935.
58. A. D. MacKerell Jr, N. Banavali and N. Foloppe, Biopolymers, 2000, 56, 257–265.
59. J. C. Phillips, R. Braun, W. Wang, J. Gumbart, E. Tajkhorshid, E. Villa, C. Chipot, R. D. Skeel, L. Kale and K. Schulten, J. Comput. Chem., 2005, 26, 1781–1802.
60. M. Nooijen and R. J. Bartlett, J. Chem. Phys., 1995, 102, 3629–3647.
61. A. I. Krylov, Annu. Rev. Phys. Chem., 2008, 59, 433–462.
62. P.-Å. Malmqvist and B. O. Roos, Chem. Phys. Lett., 1989, 155, 189–194.
63. H.-J. Werner, P. J. Knowles, F. R. Manby, J. A. Black, K. Doll, A. Heßelmann, D. Kats, A. Köhn, T. Korona, D. A. Kreplin, Q. Ma, I. Miller, F. Thomas, A. Mitrushchenkov, K. A. Peterson, I. Polyak, G. Rauhut and M. Sibaev, J. Chem. Phys., 2020, 152, 144107.
64. H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby and M. Schütz, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 242–253.
65. H.-J. Werner, P. J. Knowles, P. Celani, W. Györffy, A. Hesselmann, D. Kats, G. Knizia, A. Köhn, T. Korona, D. Kreplin, R. Lindh, Q. Ma, F. R. Manby, A. Mitrushenkov, G. Rauhut, M. Schütz, K. R. Shamasundar, T. B. Adler, R. D. Amos, S. J. Bennie, A. Bernhardsson, A. Berning, J. A. Black, P. J. Bygrave, R. Cimiraglia, D. L. Cooper, D. Coughtrie, M. J. O. Deegan, A. J. Dobbyn, K. Doll, M. Dornbach, F. Eckert, S. Erfort, E. Goll, C. Hampel, G. Hetzer, J. G. Hill, M. Hodges, T. Hrenar, G. Jansen, C. Köppl, C. Kollmar, S. J. R. Lee, Y. Liu, A. W. Lloyd, R. A. Mata, A. J. May, B. Mussard, S. J. McNicholas, W. Meyer, T. F. Miller III, M. E. Mura, A. Nicklass, D. P. O'Neill, P. Palmieri, D. Peng, K. A. Peterson, K. Pflüger, R. Pitzer, I. Polyak, M. Reiher, J. O. Richardson, J. B. Robinson, B. Schröder, M. Schwilk, T. Shiozaki, M. Sibaev, H. Stoll, A. J. Stone, R. Tarroni, T. Thorsteinsson, J. Toulouse, M. Wang, M. Welborn and B. Ziegler, MOLPRO, version, a package of ab initio programs, see https://www.molpro.net.
66. B. O. Roos, Adv. Chem. Phys., 1987, 69, 399–445.
67. M. S. Gordon, L. Slipchenko, H. Li and J. H. Jensen, Annu. Rep. Comput. Chem., 2007, 3, 177–193.
68. M. S. Gordon, D. G. Fedorov, S. R. Pruitt and L. V. Slipchenko, Chem. Rev., 2012, 112, 632–672.
69. M. S. Gordon, Q. A. Smith, P. Xu and L. V. Slipchenko, Annu. Rev. Phys. Chem., 2013, 64, 553–578.
70. A. DeFusco, N. Minezawa, L. V. Slipchenko, F. Zahariev and M. S. Gordon, J. Phys. Chem. Lett., 2011, 2, 2184–2192.
71. L. V. Slipchenko, J. Phys. Chem. A, 2010, 114, 8824–8830.
72. D. Kosenkov and L. V. Slipchenko, J. Phys. Chem. A, 2011, 115, 392–401.
73. E. Epifanovsky, A. T. Gilbert, X. Feng, J. Lee, Y. Mao, N. Mardirossian, P. Pokhilko, A. F. White, M. P. Coons and A. L. Dempwolff, et al. , J. Chem. Phys., 2021, 155, 084801.
74. B. Mennucci, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 386–404.
75. B. Nestmann and S. D. Peyerimhoff, J. Phys. B: At. Mol. Phys., 1985, 18, 615.
76. A. Whitehead, R. Barrios and J. Simons, J. Chem. Phys., 2002, 116, 2848–2851.
77. T. H. Dunning Jr, J. Chem. Phys., 1989, 90, 1007–1023.
78. H. S. Taylor, Adv. Chem. Phys., 1970, 91–147.
79. M. Pyla and S. Matsika, J. Chem. Phys., 2024, 161, 154306.
80. E. Epifanovsky, K. Kowalski, P.-D. Fan, M. Valiev, S. Matsika and A. I. Krylov, J. Phys. Chem. A, 2008, 112, 9983–9992.
81. C. S. Anstöter, S. Abou-Hatab, M. Thodika and S. Matsika, J. Phys. Chem. A, 2022, 126, 8508–8518.
82. D. Ghosh, O. Isayev, L. V. Slipchenko and A. I. Krylov, J. Phys. Chem. A, 2011, 115, 6028–6038.
83. T. J. Boerner, S. Deems, T. R. Furlani, S. L. Knuth and J. Towns, Practice and Experience in Advanced Research Computing, 2023, pp. 173–176.

---

Footnote

† Electronic supplementary information (ESI) available: Coordinates of uracil, additional tables, active space orbitals involved in CASSCF calculation and stabilization plots. See DOI: https://doi.org/10.1039/d4cp04333b

---