Exploring water radiolysis in proton cancer therapy: Time-dependent, non-adiabatic simulations of H+ + (H2O)1-6

To elucidate microscopic details of proton cancer therapy (PCT), we apply the simplest-level electron nuclear dynamics (SLEND) method to H+ + (H2O)1-6 at ELab = 100 keV. These systems are computationally tractable prototypes to simulate water radiolysis reactions—i.e. the PCT processes that generate the DNA-damaging species against cancerous cells. To capture incipient bulk-water effects, ten (H2O)1-6 isomers are considered, ranging from quasi-planar/multiplanar (H2O)1-6 to “smallest-drop” prism and cage (H2O)6 structures. SLEND is a time-dependent, variational, non-adiabatic and direct method that adopts a nuclear classical-mechanics description and an electronic single-determinantal wavefunction in the Thouless representation. Short-time SLEND/6-31G* (n = 1–6) and /6-31G** (n = 1–5) simulations render cluster-to-projectile 1-electron-transfer (1-ET) total integral cross sections (ICSs) and 1-ET probabilities. In absolute quantitative terms, SLEND/6-31G* 1-ET ICS compares satisfactorily with alternative experimental and theoretical results only available for n = 1 and exhibits almost the same accuracy of the best alternative theoretical result. SLEND/6-31G** overestimates 1-ET ICS for n = 1, but a comparable overestimation is also observed with another theoretical method. An investigation on H+ + H indicates that electron direct ionization (DI) becomes significant with the large virtual-space quasi-continuum in large basis sets; thus, SLEND/6-31G** 1-ET ICS is overestimated by DI contributions. The solution to this problem is discussed. In relative quantitative terms, both SLEND/6-31* and /6-31G** 1-ET ICSs precisely fit into physically justified scaling formulae as a function of the cluster size; this indicates SLEND’s suitability for predicting properties of water clusters with varying size. Long-time SLEND/6-31G* (n = 1–4) simulations predict the formation of the DNA-damaging radicals H, OH, O and H3O. While “smallest-drop” isomers are included, no early manifestations of bulk water PCT properties are observed and simulations with larger water clusters will be needed to capture those effects. This study is the largest SLEND investigation on water radiolysis to date.

The selected (H 2 O) 1-6 series contains the initial terms of the long progression from molecular H 2 O to bulk water. Specifically, the first terms in this progression, (H 2 O) 1-5, have monoand di-cyclic quasi-planar/multiplanar structures [33][34][35] exhibiting no "bulky" shapes (cf. Fig  1), but two isomers in the last term-the prism and cage isomers of (H 2 O) 6 -have multi-cyclic, three-dimensional structures exhibiting drop-like shapes [34,35] (cf. Fig 1). In fact, these two (H 2 O) 6 isomers, particularly the prism, are considered the smallest possible drops of water [34,35]-i.e. the minimum water structures that manifest the three-dimensional hydrogen-bond structure [35] and solubility properties [36,37] of bulk water. Based on those facts, we expected that this series would reveal the earliest manifestations of bulk water effects on PCT properties; however, the present results do not display such manifestations and suggest that even larger water clusters should be considered (cf. Results and Discussion Section). Despite that outcome, all the predicted properties and reactions of H + + (H 2 O) 2-6 have never been measured or calculated before; therefore, the present results are truly predictive and fill a gap in the medical physics literature [16]. Furthermore, these results are important to understand PCT more thoroughly and to model water radiolysis processes [15][16][17][18][19][20] and radiation dosages [12][13][14][16][17][18][19][20] with MC methods. Finally, it should be noticed that our simulated phenomena are the first processes occurring upon a short-time direct collision of a proton with moderate size clusters. Other post-collision phenomena contributing to PCT such as local temperature increases [2][3][4] and hypothesized shock waves [38] in water require for their modelling longer simulation times, larger clusters and even different theoretical methods; therefore, those phenomena are not reproduced by the current simulations.

Methodology
The END theory and its SLEND version have been reviewed in detail in Refs. [26,29,39]; therefore, we provide a brief outline of them. END is a variational, time-dependent, direct, and non-adiabatic dynamical method [26,29]. END prescribes a total trial wavefunction jC END Total i ¼ jC END N i jC END e i, which consists of nuclear jC END N i and electronic jC END e i wavefunctions, and treats jC END Total i under the time-dependent variational principle (TDVP) [40]. The various versions of END differ in the kind of descriptions for the nuclei and electrons (e.g., density functional theory [30] [26] for electrons). In SLEND, the nuclear wavefunction jC SLEND N i for a system having N N nuclei is the product  6 ] optimized geometries of the water cluster isomers selected for this study. When two or more isomers of a given water cluster are considered, they are depicted/listed in the order of increasing energy; the first depicted/listed isomer is the lowest-energy isomer in its whole known series with the present theory [e.g. 6a is the lowest-energy (H 2 O) 6 structure out of 12 known (H 2 O) 6 isomers [34]]. 6a and 6b are the prism and cage isomers of (H 2 O) 6 with average positions R A (t), average momenta P A (t) and widths ΔR A . To lower computational cost, the zero-width limit, ΔR A ! 0, is applied to the nuclear wave packets after constructing the SLEND quantum Lagrangian (cf. next paragraph). That procedure generates a classical nuclear dynamics but with full retention of the nucleus-electron non-adiabatic coupling terms (cf. Eq 5).
As proven previously [6,26,27] with complex-valued molecular coefficients {z ph (t)}. The MSOs are obtained at initial time at the UHF level. The MSOs are constructed with travelling atomic basis set functions 0 A i ðr i ; R A ; P A Þi.e., Slater-type orbitals in terms of contracted Gaussian-type orbitals on the moving nuclear centers R A (t) and augmented with electron translation factors (ETFs) [41] to include explicit nuclear momenta P A (t) effects. MSOs and DSOs are spin-unrestricted; therefore, the unrestricted determinant |z(t);R(t),P(t)i can reasonably describe bond-breaking/-forming processes. The Thouless representation provides a non-redundant and singularity-free parameterization of an evolving single-determinantal state [26,29,40]. The SLEND dynamical equations are obtained according to the TDVP [40]. First, the quantum Lagrangian [26,29,40] L SLEND ¼ hC SLEND Total j ði=2Þ ð@ * =@t À @ =@tÞ ÀĤ Total j C SLEND Total i= hC SLEND Total j C SLEND Total i À 1 is constructed and then the zero-width limit is applied to the nuclear wave packets. Subsequently, the stationary condition is imposed to the quantum action [26,29,40] The resulting SLEND dynamical equations are [26,29]: where E total is the total (nuclear and electronic) energy and the dynamic metric matrices are [26,29] ðC XY Þ ik; jl ¼ À 2Im where X ik and Y jl denote either R A = i,k or P A = j,l and S = hz(t), R 0 (t), P 0 (t)|z(t), R(t), P(t)i. C R and C RR are equivalent to the ordinary non-adiabatic coupling terms [42]. Neglect of these terms seriously impairs the accuracy of the SLEND non-adiabatic dynamics [43]; therefore, those terms are kept in the present calculations. However, to lower computational cost, ETFs are not included in the present basis set so that C P = C PR = 0. The accelerated SLEND equations thus obtained have been successfully applied to various reactive and non-adiabatic processes, from a few eV [44] to the keV regime [45][46][47][48] and up to 900 keV [48] (cf. also the seminal END Ref. [29], especially page 948, where this approximation is applied exhaustively to high-energy non-adiabatic processes). At initial time, the reactants are prepared with positions fR i A g, momenta fP i A g and Thouless electronic state |z (i) = 0, R (i) i = |0i, i.e. the UHF ground state of the reactants' super-molecule. When the reaction accesses the non-adiabatic regime, z(t) 6 [1][2][3][4][5] ] are adopted. These basis sets provide good water clusters' geometries and energies (cf. Refs. [33,34]; for these basis sets' dynamical performance, cf. Electron Transfer Properties and Further Analysis Sub-Sections). Numerous theoretical [33,34,[54][55][56][57][58][59][60][61][62][63] and experimental [35,[64][65][66] studies have been devoted to determine the geometries and energies of water clusters because these systems are prototypes to study the structural [35] and solubility properties [36,37] of bulk water. In fact, the scientific literature about water clusters is vast and growing. Therefore, for brevity's sake, we limit ourselves to cite herein only the water clusters studies closely related to this investigation. It is well-known that the (H 2 O) 2-3 clusters present one isomer each [34,55]; however, the (H 2 O) n , n ! 4, clusters present a variable number of isomers that rapidly increase with the cluster size n [34,55]. (H 2 O) 1-5 present mono-and di-cyclic quasi-planar/multiplanar structures [33][34][35], whereas (H 2 O) n , n ! 6, present multi-cyclic, three-dimensional structures in addition to quasi-planar/multiplanar ones [34,35] (cf. Fig 1). The three-dimensional (H 2 O) 6 isomers (e.g. the prism and cage isomers named hexamer a and b in Fig 1) are considered the smallest drops of water [34,35]. In general, theory and experiments agree in regard to the structures and relative energies of the (H 2 O) 1-5 isomers, but discrepancies arise in the energy orders of the (H 2 O) n , n ! 6, isomers [35,56]. Recent spectroscopy experiments have identified the cage isomer as the lowest-energy (H 2 O) 6 structure [35,66]. However, ab initio calculations at various levels of accuracy have disagreed on whether the prism [34,[57][58][59], the cage [60,61], both of them [56], or the chair [62] isomer(s) is(are) the lowest-energy (H 2 O) 6 structure(s). Ultimately, the most accurate calculations with the coupled-cluster with singles, doubles and perturbative triples [CCSD(T)] method have identified the prism isomer as the lowest-energy (H 2 O) 6 structure [58].

Fig 2. H + + (H 2 O) 1-6 initial conditions.
A given water cluster (not depicted for clarity's sake) is placed at rest with its center of mass at the origin of the coordinate axes and with its major (pseudo-)plane of symmetry with maximum coincidence with the x-y plane. The H + projectile is initially prepared with position and momentum R 0 H þ and P 0 H þ and with impact parameter b (Panel I). Different projectile-target relative orientations Ω = (α, β, γ) are generated by rotating R 0 H þ and P 0 H þ through the extrinsic Euler angles γ (Panel I), β (Panel II), and α (Panel III) around the space-fixed z, y, and z axes, respectively. The definite initial conditions of the H + projectile to start the simulations, R i H þ and P i H þ , are shown in Panel IV (cf. text for more details). https://doi.org/10.1371/journal.pone.0174456.g002 and hexamer c since E Cage Hexamer b À E Prism Hexamer a = 4.2 kJ/mol and E Hexamer c À E Prism Hexamer a = 15.4 kJ/mol with UHF/6-31G Ã ). Notice that the UHF prism isomer is the lowest-energy (H 2 O) 6 structure as predicted by CCSD(T) [58]. The first depicted/listed isomer is also the lowest-energy structure in its whole set of isomers [34]. The present UHF calculations of (H 2 O) [1][2][3][4][5][6] are not intended to contribute to the resolution of the discrepancies about the (H 2 O) 6 energies because these calculations do not attain the accuracy of CCSD(T) [58,59]. Instead, these optimizations provide a good description of the water clusters [33,34] for subsequent SLEND/6-31G Ã and /6-31G ÃÃ simulations.

Initial states preparation and simulation times
Once the water clusters are optimized at the UHF/6-31G Ã and /6-31G ÃÃ levels, the super-molecular systems H + + (H 2 O) 1-6 are assembled for the initial conditions (cf. Fig 2). The (H 2 O) 1-6 targets are prepared at rest in their equilibrium geometries with their centers of mass placed at the origin of the laboratory-frame coordinate axes; the (H 2 O) 1-6 major (pseudo-)planes of symmetry are placed with maximum coincidence with the x-y plane. The H + projectile is first prepared with position R 0 H þ ¼ ðb ! 0; 0; þ 30 a:u:Þ and momentum P 0 where b ! 0 is the projectile impact parameter measured from the (H 2 O) 1-6 centers of mass, and p z H þ corresponds to E Lab = 100 keV (cf. Fig 2, panel I). Having set R 0 H þ and P 0 H þ , various projectile-target relative orientations can be generated by rotating a (H 2 O) 1-6 target according to ordinary Euler angles prescriptions [67]. However, such a procedure involves the electronic re-optimization of the (H 2 O) 1-6 targets at each new orientation. Therefore, since the H + bare ion requires no electronic optimization, we adopted the easier and equivalent procedure of keeping a (H 2 O) 1-6 target fixed while rotating the H + projectile around. The definite initial conditions of the H + projectile R i H þ and P i H þ are therefore obtained by rotating R 0 H þ and P 0 H þ through the extrinsic Euler angles [67] in the order: 1 st , 0 0 γ < 360 0 , 2 nd , 0 0 β 180 0 , and 3 rd , 0 0 α < 360 0 , Table 1 Computer simulations of water radiolysis around the space-fixed z, y, and z axes [67], respectively (cf. Fig 2, panels I, II and III for each angle rotation). α and γ are measured from the +x axis employing the R 0 H þ and P 0 H þ projections on the x-y plane and β is measured directly from the +z axis (cf. Fig 2). These rotations define projectile-target relative orientations O i = (α i , β i , γ i ). The definite initial conditions of the H + projectile for the simulations, R i H þ and P i H þ , are shown in Fig 2, panel IV. α and β determine the direction of an axis of incidence whereby an incoming H + trajectory runs parallel to that axis with a lateral separation b ! 0; γ is the polar angle of that trajectory around the axis of incidence. In all simulations, the selected values of the Euler angles correspond to a 60-point grid: O i , 1 i 60, developed in Ref. [68]. This grid displays a uniform sampling of the orientation space and provides a numerical quadrature that ensures the invariance of Euler-angle integrals under several rotation operations (e.g. Wigner D-matrices satisfy For a given orientation O i , b is varied according to the grids defined in Table 1. Simulations for the calculation of 1-ET ICSs utilize the grids denoted as [b] 1 in Table 1 Table 1 and run for a total time of 1,000 a.u. (24.19 fs); this much longer simulation time permits the manifestation of post-collision fragmentations that have longer time scales than those of the 1-ET processes. The described initial conditions generate a total of 25,020 trajectories to complete the H + + (H 2 O) 1-6 study.

Final states analysis and properties calculations
By the end of a simulation, various auxiliary codes in the PACE package identify and analyze the final products and calculate dynamical properties. The most important properties calculated herein are the cluster-to-proton 1-ET total ICSs, σ 1−ET , corresponding to H + + (H 2 O) 1-6 ! H + different cluster products [69] s 1À ET ¼ 1 4p where P 1−ET (α, β, γ, b) is the probability of a 1-ET process from a bound electronic state of the cluster to a bound electronic state of the projectile, henceforth named bound-to-bound 1-ET for brevity's sake. Notice that for atom-atom collisions involving spherically symmetric poten- , and Eq (7) simplifies to the familiar expression: [42]. In the present systems, an outgoing projectile can in practice capture up to two electrons: H + + A ! H 1−n + A −n , 0 n 2, because the probability of forming unstable H 1−n with n ! 3 is negligible. Under these conditions, P n−ET (α, β, γ, b), 0 n 2, are [69] P 0À ET ða; b; g; bÞ ¼ ð1 À N a Þð1 À N b Þ;
A to D [70][71][72][73], respectively. In addition, Table 2 includes results from two alternative theories, namely, Theory A: the basis generator method [23] (BGM), and Theory B: the continuum distorted wave-eikonal initial state (CDW-EIS) approximation [22]. From Table 2, one finds that the average experimental ICS, s Exp: 1À ET , and its average relative error, e Exp. , are 1.27 Å 2 and ± 10.62%, respectively. The theoretical s Theo: 1À ET 's and their average relative deviations D Theo: from the experimental values are: 1.54 Å 2 and +21.8% for SLEND/6-31G Ã , 1.00 Å 2 and +21.0% for BGM, and 0.589 Å 2 and -53.4% for CDW-EIS. Only the BGM result is within the error bars of one experiment, Exp. D [73], with Δ Theo. = 11.5%, but it lies on the lowest fringe part of the error bar range. The SLEND/6-31G Ã result is very close to the result from Exp. C [72] with Δ Theo. = 11.6% and not far from entering the upper part of the error bar range. In absolute Table 2 quantitative terms, the BGM and SLEND/6-31G Ã results are at the same level of accuracy and their agreement with the experimental data should be deemed satisfactory given the difficulty to both measure and predict the present 1-ET processes. Deviations of the obtained magnitude are not uncommon in measurements and predictions of similar complex processes (cf. Ref. [27] for the case of one experiment and four different theories including SLEND, where even higher deviations are observed). The CDW-EIS result compares less favorably with the experimental ones, being roughly a half of its experimental counterparts (D Theo: = -53.4%). The SLEND/6-31 ÃÃ result also compares less favorably with the experimental ones, but, opposite to CDW-EIS, its value is roughly twice as much as the experimental one (D Theo: = 63.0%). The reason and remediation of the SLEND/6-31 ÃÃ σ 1−ET overestimation will be discussed in detail in the Further Analysis Subsection. It suffices to say here that this overestimation results from the σ 1−ET contamination with electron transfers to the continuum of unbound states. The SLEND σ 1−ET for the polymeric systems (n = 2-6) are listed in Table 3. Unfortunately, to the best of our knowledge, there are no alternative experimental or theoretical ICSs for H + + (H 2 O) 2-6 ; therefore, current SLEND σ 1−ET for these systems are predictive. To facilitate the comparison among all the considered σ 1−ET , we plot them as a function of the cluster size n in Fig 3. There, each set of SLEND/6-31G Ã and /6-31G ÃÃ σ 1−ET is fit to the scaling formulae σ 1−ET (n) = cn 2/3 , where c are fitting coefficients reported in Fig 3. These formulae are by no means arbitrary because they can be justified on physical grounds as follows. The volume V(n) of the (H 2 O) n clusters should be approximately proportional to their size n, V(n) / n. If the clusters are represented by the minimal spheres enclosing all their atoms, then their volume V(n) and external area A(n) are V ðnÞ ¼ ð4=3Þ p R 3 n and AðnÞ ¼ 4p R 2 n , respectively, where R n is the radius of the (H 2 O) n sphere; therefore, A(n) / V (n) 2/3 / n 2/3 . In turn, the σ 1−ET are proportional to the effective external area A(n) of the (H 2 O) n exposed to the incident H + for ET processes [42]; therefore, σ 1−ET / A(n) / n 2/3 ) σ 1−ET (n) = cn 2/3 . In relative quantitative terms, the SLEND/6-31G Ã and /6-31G ÃÃ σ 1−ET fit remarkably well into the physically justified formulae σ 1−ET (n) = cn 2/3 with correlation factors R 2 = 0.983 in both cases (cf. Fig 3). This indicates that regardless of their absolute quantitative performance, the SLEND σ 1−ET scale correctly Table 3. SLEND cluster-to-proton 1-electron-transfer integral cross sections σ 1−ET for H + + (H 2 O) n , n = 2-6, at E Lab = 100 keV. Cf. Fig 1 for the structures of the water cluster isomers. The error of the SLEND integral cross sections is ± 0.005 Ä . Computer simulations of water radiolysis with the number of water molecules in the clusters. Therefore, with these fitting formulae, one can estimate the σ 1−ET of the immediately larger clusters: e.g., s 6À 31G Ã 1À ET = 5.96 and 6.51 Å for n = 7 and 8, respectively; these estimated values should be interpreted as average σ 1−ET over the various (H 2 O) 7 and (H 2 O) 8 [34] isomers, respectively. Inspection of Tables 2 and 3 and Fig 3  reveals that the SLEND/6-31G ÃÃ σ 1−ET are always higher in value than the SLEND/6-31G Ã ones for each cluster as was the case with the H 2 O monomer. Inspection of Table 3 and Fig 3 reveals that the SLEND σ 1−ET for the various isomers appearing at a given n ! 4 do not significantly differ in their values; this implies that these σ 1−ET are rather insensitive to the varying isomers' structures as targets. A similar finding was observed in the σ 1−ET of H + + DNA/RNA bases at E Lab = 80 keV [27], where similar bases differing in their structure even more than isomers exhibited close values of σ 1−ET . We expected that various σ 1−ET (n) = cn 2/3 formulae differing in their coefficient c would exclusively connect values from different sets of clusters-e.g. a single  [71], C [72] and D [73], Theory A: basis generator method (BGM) [23], Theory B: continuum distorted wave-eikonal initial state (CDW-EIS) approximation [22]]. SLEND values are fit to the scaling formula σ 1−ET (x) = cn 2/3 . The error of the SLEND ICSs is ± 0.005 Ä . The errors from the Theory A and B results compared herein were not reported.

Water Cluster Basis Set 1-electron ICS (Å 2 )
https://doi.org/10.1371/journal.pone.0174456.g003 Computer simulations of water radiolysis σ 1−ET (n) = cn 2/3 formula might have only fit well with results from the quasi-planar clusters (monomer, dimer, trimer, tetramer a, pentamer a and hexamer c), another single formula with other type of clusters, etc., but that is not case. In fact, the uniformity among the SLEND σ 1−ET (n) values for isomers at a given n led us to fit all of them with a single formula per basis set. Furthermore, we expected that the σ 1−ET of the drop-like prism and cage (H 2 O) 6 isomers would differ sharply from the σ 1−ET of the quasi-planar/multiplanar isomers so that it would be impossible to fit drop-like and non-drop-like σ 1−ET with a single formulae. Such a hypothetical fitting failure would manifest as a "phase transition" discontinuity from non-drop-like to drop-like σ 1−ET . However, no such "phase transition" is observed in the selected series of (H 2 O) 1-6 isomers. Thus, unlike the case of structural and solubility properties [35][36][37], the "magic number" of six waters in the prism and cage (H 2 O) 6 isomers do not bring about any hint of water radiolysis processes in bulk water. Likely, an extension of the current (H 2 O) 1-6 series with the (H 2 O) 7-20 isomers may bring about some type of bulk-water manifestations.
Total 1-ET ICSs σ 1−ET are not very detailed properties and cannot reveal some dynamical details of 1-ET processes. Therefore, Figs 4 and 5 show the orientation-averaged 1-ET probabilities P 1À ET ðbÞ and their b -weighted counterparts bP 1À ET ðbÞ vs. b for the considered (H 2 O) 1-6 isomers. With both basis sets, the P 1À ET ðbÞ are high in value at small impact parameters (roughly, b 6 a.u.) corresponding to close projectile-cluster encounters but they decrease rapidly at larger impact parameters. The P 1À ET ðbÞ vs. b plots show a variable number of maximum peaks (from one to three) depending on the considered water cluster. The bP 1À ET ðbÞ vs. b plots exhibit similar patterns to those of the P 1À ET ðbÞ but modulated by the b value. As the integrands of the σ 1−ET , the bP 1À ET ðbÞ plots indicate that the most important contributions to the σ 1−ET come from 1-ET processes arising from intermediate impact parameters (roughly, 2 a.u. b 9 a.u.).

Fragmentation reactions
Unlike time-independent scattering methods [21,22,74] applicable to PCT, SLEND simulations can reveal the reactants-to-products time evolution of the chemical reactions underlying  Fig 1 (1, 2 . . . 6b, 6c). https://doi.org/10.1371/journal.pone.0174456.g004 Computer simulations of water radiolysis the 1-ET processes. To study those reactions, the SLEND/6-31G Ã simulations to calculate σ 1−ET of H + + (H 2 O) 1-4 (only tetramer a for n = 4) were prolonged from their simulation times of 30.00 a.u. (0.7257 fs) to 1,000 a.u. (24.19 fs) using the impact parameter grids [b] 2 in Table 1. This much longer simulation time allows for the manifestation of fragmentation reactions that occur at longer time scales. In fact, no fragmentation was observed within the original time of 30.00 a.u. As Eq 6 shows, the final SLEND electronic wavefunction |z(t),R(t)i is a superposition of various UHF states corresponding to various products' channels [69,75] where H q 1 ¼0;þ1;À 1 proj: is the incoming/outgoing projectile; these channel states occur with different probabilities [69,75]. Therefore, when Eq 9 is applied to all the well-separated fragments at final time, N α and N β and their corresponding charges q i (e.g. q 1 = 1−N α −N β for the final H q 1 proj: ) are not necessarily integer numbers corresponding to canonical chemical species (e.g. H q 1 ¼þ1 proj , H q 1 ¼0 proj and H q 1 ¼À 1 proj ) but fractional numbers as the averages of the number of electrons and charges over the channels' probabilities. This was always the case in all the present simulations not leading to clusters fragmentations (i.e. H þ1 proj: þ ðH 2 OÞ n ! H q 1 proj: þ ðH 2 OÞ 1À q 1 n , with q i continuously varying in the range −1 q i +1). However, as in previous SLEND studies of H + + (H 2 O) 1-4 at E Lab = 1 keV [6], the present simulations leading to clusters fragmentations always bring about outgoing projectiles H q 1 proj: and clusters fragments A q 2 with q 1 = +1 and q 2 = 0, respectively (cf. Fig 6 caption). Thus, the present SLEND simulations predict that the fragmentation channel leading to outgoing H + and neutral fragments predominates over the others. However, it is known experimentally that proton-water collisions lead to fragmentations into ions [76][77][78]. SLEND can properly describe fragmentations into ions as shown in previous studies [e.g. cf. Refs. [79,80]]. Therefore, to allow the manifestation of those types of fragmentations here, it will be necessary to prolong even further the simulation time of the present calculations or, more likely, increase Water cluster isomers are denoted with the number code in Fig 1 (1, 2 . . . 6b, 6c). To illustrate the predicted fragmentations, we present a few animation stills from some representative simulations. Computer simulations of water radiolysis studies [6,75], Mulliken populations were good predictors for the time-evolution of the electrons over atoms and fragments. Furthermore, when all the fragments are well-separated at final time, the Mulliken populations converge to the unequivocal N α and N β in Eq 9 (cf. Fig 6  caption). In Fig 6, the colored spheres represent the classical nuclei (white = H and red = O), and the colored clouds represent selected electron density iso-surfaces (from red = lowest density to blue = highest density Further analysis and improvements SLEND/6-31G ÃÃ σ 1−ET compares unsatisfactorily with experiments in contrast to SLEND/6-31G Ã σ 1−ET . Indeed, it is surprising that SLEND/6-31G ÃÃ performs worse than SLEND/6-31G Ã in these Computer simulations of water radiolysis calculations since common knowledge dictates that the 6-31G ÃÃ basis set is better than the 6-31G Ã one. 6-31G ÃÃ is constructed from 6-31G Ã by augmenting the latter with p-type basis functions on the hydrogen atoms; thanks to this augmentation, 6-31G ÃÃ provides better time-independent molecular properties than 6-31G Ã . For instance, the UHF/6-31G ÃÃ energies and geometries of (H 2 O) 1-6 are more accurate than the UHF/6-31G Ã ones. However, this comparative time-independent performance does not necessarily extend to time-dependent calculations since these basis sets were designed to calculate static properties. To explain the SLEND/6-31G ÃÃ σ 1−ET overestimation, one should remember that the main component in the bound-to-bound SLEND σ 1−ET is the 1-ET probability P 1−ET (cf. Eqs 7 and 8). However, as derived in Ref. [69] and supposed in previous SLEND studies [29,75], the P n−ET , 0 n 2, in Eq 8 assume that the probabilities of ETs from the target A to the projectile H + are dominated by transitions with electrons transferring into the localized, discrete, bound states of H: H + + A: ! HÁ + AÁ + [pure charge-transfer (CT) processes]; instead, transitions with electrons scattering into the delocalized, continuous, unbound states of H are considered negligible: H + + A: ! (H + + e − ) + AÁ + [direct ionization (DI) processes]. Typical quantum chemistry basis sets, such as 6-31G Ã and 6-31G ÃÃ , are ultimately based on localized primitive Gaussian functions so that occupied spin-orbitals {ψ h } below the Fermi level represent localized, bound states in the discrete part of the spectrum. However, as a by-product of the UHF procedure, diffuse virtual spin-orbitals {ψ p } above the Fermi level approximately represent some of the delocalized, unbound states in the continuous part of the spectrum. Therefore, the virtual space constitutes the so-called quasi-continuum that may accommodate DI processes. However, if the basis set is not large, the DI contributions of a small quasi-continuum to the ET processes become negligible in comparison to the CT contributions of the occupied space. Under those conditions, the ET probabilities P n−ET in Eq 8 basically correspond to boundto-bound (occupied-space-to-occupied-space) CT processes. For that reason, with relatively small basis sets, those P n−ET consistently rendered correct bound-to-bound CT σ 1−ET in various systems [26,27,29]. However, if the basis set is large, the DI contributions of an enlarged quasi-continuum may become substantial. If so, the P n−ET in Eq 8 and resulting σ n−ET no longer correspond to pure bound-to-bound CT processes since they get contaminated with bound-to-quasi-continuum DI contributions. Therefore, one can hypothesize that SLEND with the smaller 6-31G Ã basis set can predict genuine CT σ n−ET via Eq 8 but not with the larger 6-31G ÃÃ one.
To verify the above hypothesis, we performed a series of SLEND simulations on the simple model system: H + + H at 40 keV E Lab 100 keV with the 6-31 ++ G ÃÃ basis set. The latter produces the best DI results for H + + H as shown shortly. However, instead of using P 1−ET in Eq 8 for σ 1−ET , the final-time electronic wavefunction jC SLEND e ðt f Þi, Eq 6, was projected onto the ground |0i and excited states |. . . h ! p . . .i of the target and the projectile. In this way, the evaluation of ET probabilities could distinguish the cases with the electron transferring into bound or unbound states of H. In the 40 E Lab 80 keV range, where experimental results are available [81], the DI ICSs, σ DI , deviate less than 10% from experimental data [81], with the best agreement at E Lab = 60 keV: s SLEND DI = 13.9Å and s EXPT: DI = 13.8 Å ) deviation Δ Theo. = 0.7%. Notably, these calculations produced accurate results even though ETFs were neglected as in Eq 4; this gives extra support to the ETFs' neglect in the H + + (H 2 O) n simulations and ruled it out as a source of the SLEND/6-31G ÃÃ σ 1−ET overestimation. For H + + H at E Lab = 100 keV, SLEND σ DI is 0.92 Å 2 . If P 1−ET in Eq 8 is used to calculate CT ICSs, σ CT , a σ DI contribution of 0.92 Å will be spuriously added to the σ CT making it overestimated. If this DI contribution is assumed to be similar to that in H + + H 2 O, the overestimated SLEND/6-31G ÃÃ s SLEND the SLEND/6-31G Ã s SLEND CT but it will be far smaller due to a smaller virtual space as suggested by previous calculations with comparable basis sets [26,27,29]. The calculation of the CT σ CT in H + + H 2 O is more complicated than that of H + + H because the former has more than one electron. For H + + H 2 O, numerous excited states |. . . h ! p . . .i from |0i forming a full CI expansion should be generated so that jC SLEND e ðt f Þi in Eq 6 can be projected on all those states. This more demanding capability is not currently available in PACE but is under development.