Home | About | Journals | Submit | Contact Us | Français |

**|**PLoS One**|**v.12(4); 2017**|**PMC5380356

Formats

Article sections

- Abstract
- Introduction and background
- Methodology
- Computational details
- Results and discussion
- Conclusions
- References

Authors

Related links

PLoS One. 2017; 12(4): e0174456.

Published online 2017 April 4. doi: 10.1371/journal.pone.0174456

PMCID: PMC5380356

Austin J. Privett,^{¤a}^{} Erico S. Teixeira,^{} Christopher Stopera,^{¤b}^{} and Jorge A. Morales^{*}^{}

Ilia Solov'yov, Editor^{}

Department of Chemistry and Biochemistry, Texas Tech University, Lubbock, Texas, United States of America

Syddansk Universitet, DENMARK

**Conceptualization:**JAM.**Data curation:**AJP EST JAM.**Formal analysis:**AJP CS EST JAM.**Funding acquisition:**JAM.**Investigation:**AJP CS EST JAM.**Methodology:**AJP JAM.**Project administration:**JAM.**Resources:**JAM.**Software:**AJP EST JAM.**Supervision:**JAM.**Validation:**AJP CS EST.**Visualization:**AJP EST.**Writing – original draft:**AJP JAM.**Writing – review & editing:**JAM.

* E-mail: ude.utt@selarom.egroj

Received 2016 October 31; Accepted 2017 March 9.

Copyright © 2017 Privett et al

This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

To elucidate microscopic details of proton cancer therapy (PCT), we apply the simplest-level electron nuclear dynamics (SLEND) method to H^{+} + (H_{2}O)_{1-6} at E_{Lab} = 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 (H_{2}O)_{1-6} isomers are considered, ranging from quasi-planar/multiplanar (H_{2}O)_{1-6} to “smallest-drop” prism and cage (H_{2}O)_{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 H_{3}O. 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.

Proton cancer therapy (PCT) is an approved type of radiotherapy that utilizes high-energy H^{+} projectiles to fight cancer [1–7]. The ultimate effect of this radiation is to damage the DNA of cancerous cells [1–7]. If left unrepaired, this damage produces various anomalies in cancerous cells that eventually lead to their death (apoptosis) [1–7]. While PCT damages both healthy and cancerous cells, the latter have a high rate of division and a reduced ability to repair damaged DNA. Thus, cancer cells are particular vulnerable to radiation attacks[8].

PCT radiation is applied as collimated beams of H^{+} projectiles at an initial kinetic energy of 70–250 MeV [1–7]. As they penetrate the patient’s body, these projectiles lose their energy through molecular interactions until they reach thermal energy in deep tissues. The radiation energy deposited in the tissues is a measure of its potential for DNA damage. Typically, a plot of the radiation energy loss *vs*. the radiation travelled distance exhibits a maximum of energy deposition [2]. Conventional X-ray therapy exhibits a broad deposition maximum not far after the photons’ penetration into the body, followed by a gradual energy loss at deeper penetrations [2]. In contrast, PCT (and other ion therapies) exhibits a sharp maximum peak known as the Bragg peak; that peak occurs just before the H^{+} projectiles are stopped in deep tissues [2]. Thus, by focusing Bragg peaks on a deep tumor, PCT inflicts a maximum DNA damage on that region and a minimum DNA damage on the surrounding healthy tissues [1–7].

In all radiotherapies, the radiation predominantly interacts with cellular H_{2}O because the latter constitutes about 70% of the human cell mass. This interaction triggers water radiolysis reactions—i.e. a series of cascade hydrolytic reactions producing DNA-damaging species. In PCT, water radiolysis generates various secondary species: (a) free radicals (e.g. H^{+} + H_{2}O → H^{+} + H·+ OH·)[2–4, 6, 9], (b) secondary ions (e.g. H^{+} + H_{2}O → 2H^{+}
*+* OH^{-})[2–4, 6, 9], (c) reactive molecules (e.g. H^{+} + 2H_{2}O → H^{+} + H_{2} + H_{2}O_{2}) [10], and (d) localized heat in the medium [2–4, 6, 9]. These species can react with H_{2}O and generate similar tertiary species and so forth. All these reactive products eventually reach cellular DNA and cause various types of damage: DNA bases’ fragmentations and deletions, and single-, double- and clustered-strand DNA breaks [2–5, 11].

PCT comprises various processes spanning different space (*l* = 10^{−10}–10^{0} m) and time (*t* = 10−21–10^{5} s) scales [2–4]. For instance, water radiolysis reactions, DNA damage at the genome level [12] and tumor remissions lie at the microscopic (10^{−10} ≤*l*≤10^{−8} m), mesoscopic (10^{−8} ≤*l*≤10^{−3} m) and macroscopic (10^{−3} ≤*l*≤10^{0} m) scales, respectively. The scale of a process determines the appropriate methods for its study. Thus, in theoretical/computational studies, microscopic water radiolysis reactions can be feasibly simulated with *ab initio* quantum-mechanics methods. In contrast, mesoscopic energy-loss and Bragg peak processes are only tractable with Monte Carlo (MC) models [12–14]. Quantum-mechanics and MC methods act in synergy: the former predict properties (e.g. reaction cross sections) required as input data for MC simulations [9, 12–20], and the latter calculate proper radiation doses for treatments [18].

Although PCT is clinically approved, various PCT details at the microscopic scale remain uncertain [1–7]. Knowledge of those details is essential for a rational design of PCT seeking to maximize its therapeutic power and minimize its side effects [1–7, 18]. While the predominant paradigm for cancer research is experimental/clinical, theoretical/computational methods can reveal microscopic details of PCT more exhaustively than experimental/clinical techniques and without putting human subjects at risk. Therefore, time-independent scattering [21–23] and time-dependent dynamics methods [6, 24–27] have been applied to computationally feasible prototypes of PCT reactions (e.g. H^{+} + (H_{2}O)_{1–4} to model water radiolysis reactions [6, 21–26, 28] and H^{+} + DNA/RNA bases to model DNA proton damage [6, 21, 22, 26, 27]).

Among quantum-mechanics methods for PCT [6, 24–27], the electron nuclear dynamics (END) theory [26, 29] offers distinctive capabilities to study PCT reactions. END is a (1) time-dependent, (2) direct and (3) non-adiabatic method to simulate chemical reactions. These attributes are valuable for PCT simulations because they afford: (1) time-dependent detail, (2) independence from predetermined potential energy surfaces, and (3) capability of describing high-energy non-adiabatic processes [e.g. electron transfers (ETs)]. Among different END versions [26, 29, 30], the simplest-level (SL) END (SLEND) describes the nuclei and electrons in terms of classical mechanics and a Thouless single-determinantal wavefunction, respectively [26, 29, 31]. Thus, SLEND possesses a suitable balance between accuracy and computational feasibility to simulate large PCT prototypes (cf. previous SLEND studies of H^{+} + (H_{2}O)_{n}, *n* = 1 [24, 28], 2 [25], and 3–4 [6, 26] and of H^{+} + DNA/RNA bases and DNA base pairs [6, 26, 27]).

Following the aforesaid precedents [6, 21–26, 28], we present herein an exploratory SLEND study of PCT water radiolysis reactions with the H^{+} + (H_{2}O)_{1-6} prototypes at *E*_{Lab} = 100 keV. This energy is selected because it corresponds to the average Bragg peak energy in bulk water [32], the medium where water radiolysis occurs. However, no *ab initio* quantum-mechanics methods can simulate bulk water due to prohibitive computational cost and, therefore, those methods treat the above-mentioned water-clusters prototypes. Surprisingly, most quantum-mechanics studies of water radiolysis have utilized the smallest prototype: H^{+} + H_{2}O [21–24, 28], whose “cluster” is the farthest from being bulk water. While these H^{+} + H_{2}O studies were indeed useful for investigating radiolysis processes [21–24, 28], they could not completely capture the processes occurring in bulk water. For that reason, previous SLEND studies explored H^{+} + (H_{2}O)_{n} prototypes with n = 2 [25] and 3–4 [6, 26]. However, these studies still concentrated on the smallest possible clusters and involved a limited number of proton-cluster orientations and simulations. Therefore, to overcome all the discussed limitations, we study herein the H^{+} + (H_{2}O)_{1-6} prototypes that include ten isomers in a larger series of clusters (H_{2}O)_{1-6} (cf. Fig 1) and involve a larger number of proton-cluster orientations (60) and simulations (25,020).

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 mono- and di-cyclic quasi-planar/multiplanar structures [33–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–20] and radiation dosages [12–14, 16–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–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.

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 $|{\mathrm{\Psi}}_{Total}^{END}\u27e9=|{\mathrm{\Psi}}_{N}^{END}\u27e9\phantom{\rule{0.25em}{0ex}}|{\mathrm{\Psi}}_{e}^{END}\u27e9$, which consists of nuclear $|{\mathrm{\Psi}}_{N}^{END}\u27e9$ and electronic $|{\mathrm{\Psi}}_{e}^{END}\u27e9$ wavefunctions, and treats $|{\mathrm{\Psi}}_{Total}^{END}\u27e9$ 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 $|{\mathrm{\Psi}}_{N}^{SLEND}\u27e9$ for a system having *N*_{N} nuclei is the product of 3*N*_{N} frozen, narrow Gaussian wave packets:

$$\u27e8\mathbf{X}|{\mathrm{\Psi}}_{N}^{SLEND}\left(t\right)\u27e9=\u27e8\mathbf{X}|\mathbf{R}\left(t\right),\phantom{\rule{0.25em}{0ex}}\mathbf{P}\left(t\right)\u27e9={\displaystyle \prod _{A=1}^{3{N}_{N}}\mathrm{exp}\left\{-{\left[\frac{{X}_{A}-{R}_{A}\left(t\right)}{2\mathrm{\Delta}{R}_{A}}\right]}^{2}+i{P}_{A}\left(t\right)\left[{X}_{A}-{R}_{A}\left(t\right)\right]\right\}}$$

(1)

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], classical nuclear dynamics does not impair the accuracy of PCT simulations because they happen at high collision energies. The SLEND electronic wavefunction $|{\mathrm{\Psi}}_{e}^{SLEND}\u27e9$ is a spin-unrestricted, single-determinantal wavefunction in the Thouless representation [31]. Specifically, taking *N*_{e} occupied {*ψ*_{h}}, and *K* – *N*_{e} virtual {*ψ*_{p}} molecular spin-orbitals (MSOs), $|{\mathrm{\Psi}}_{e}^{SLEND}\u27e9$ is [31]:

$$|{\mathrm{\Psi}}_{e}^{SLEND}\u27e9=|\mathbf{z}\left(t\right);\mathbf{R}\left(t\right),\mathbf{P}\left(t\right)\u27e9=\mathrm{det}\left\{{\chi}_{h}\left[{\mathbf{x}}_{h};\phantom{\rule{0.25em}{0ex}}\mathbf{z}\left(t\right),\mathbf{R}\left(t\right),\mathbf{P}\left(t\right)\right]\right\}=\text{exp}\left[{\displaystyle \sum _{h=1,\phantom{\rule{0.25em}{0ex}}p={N}_{e}+1}^{{N}_{e},\phantom{\rule{0.25em}{0ex}}K}{z}_{ph}\left(t\right){b}_{p}^{\u2020}{b}_{h}}\right]|0\u27e9;$$

(2)

where |0⟩ = |*ψ*_{Ne}…*ψ*_{1}⟩ is an unrestricted Hartree-Fock (UHF) reference state and {*χ*_{h}} are dynamical spin-orbitals (DSOs)

$${\chi}_{h}\left[\mathbf{x};\phantom{\rule{0.25em}{0ex}}\mathbf{z}\left(t\right),\mathbf{R}\left(t\right),\mathbf{P}\left(t\right)\right]={\psi}_{h}\left[\mathbf{x};\phantom{\rule{0.25em}{0ex}}\mathbf{R}\left(t\right),\mathbf{P}\left(t\right)\right]+{\displaystyle \sum _{p={N}_{e}+1}^{K}{z}_{ph}\left(t\right){\psi}_{p}}\left[\mathbf{x};\phantom{\rule{0.25em}{0ex}}\mathbf{R}\left(t\right),\mathbf{P}\left(t\right)\right];\phantom{\rule{0.25em}{0ex}}\left(1\le h\le {N}_{e}\right)$$

(3)

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 *ϕ*_{Ai}(**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*) 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}=\u27e8{\mathrm{\Psi}}_{Total}^{SLEND}\phantom{\rule{0.25em}{0ex}}|\phantom{\rule{0.25em}{0ex}}\left(i/2\right)\phantom{\rule{0.25em}{0ex}}(\stackrel{\rightharpoonup}{\partial}/\partial t-\overleftarrow{\partial}/\partial t)-{\widehat{H}}_{Total}\phantom{\rule{0.25em}{0ex}}|\phantom{\rule{0.25em}{0ex}}{\mathrm{\Psi}}_{Total}^{SLEND}\u27e9/\u27e8{\mathrm{\Psi}}_{Total}^{SLEND}\phantom{\rule{0.25em}{0ex}}|\phantom{\rule{0.25em}{0ex}}{{\mathrm{\Psi}}_{Total}^{SLEND}\u27e9}^{-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] *A*_{SLEND}, $\delta {A}_{SLEND}=\delta {\displaystyle {\int}_{{t}_{1}}^{{t}_{2}}{L}_{SLEND}(t)dt}=0$. The resulting SLEND dynamical equations are [26, 29]:

$$\left[\begin{array}{cccc}i\mathbf{C}& \mathbf{0}& i{\mathbf{C}}_{\mathbf{R}}& i{\mathbf{C}}_{\mathbf{P}}\\ & & & \\ \mathbf{0}& -i{\mathbf{C}}^{*}& -i{\mathbf{C}}_{\mathbf{R}}^{*}& -i{\mathbf{C}}_{\mathbf{P}}^{*}\\ & & & \\ i{\mathbf{C}}_{\mathbf{R}}^{\u2020}& -i{\mathbf{C}}_{\mathbf{R}}^{T}& {\mathbf{C}}_{\mathbf{R}\mathbf{R}}& -\mathbf{I}+{\mathbf{C}}_{\mathbf{R}\mathbf{P}}\\ & & & \\ i{\mathbf{C}}_{\mathbf{P}}^{\u2020}& -i{\mathbf{C}}_{\mathbf{P}}^{T}& \mathbf{I}+{\mathbf{C}}_{\mathbf{P}\mathbf{R}}& {\mathbf{C}}_{\mathbf{R}\mathbf{P}}\end{array}\right]\left[\begin{array}{c}\frac{d\mathbf{z}}{dt}\\ \frac{d{\mathbf{z}}^{*}}{dt}\\ \frac{d\mathbf{R}}{dt}\\ \frac{d\mathbf{P}}{dt}\end{array}\right]=\left[\begin{array}{c}\frac{\partial {E}_{Total}}{\partial {\mathbf{z}}^{*}}\\ \frac{\partial {E}_{Total}}{\partial {\mathbf{z}}^{}}\\ \frac{\partial {E}_{Total}}{\partial \mathbf{R}}\\ \frac{\partial {E}_{Total}}{\partial \mathbf{P}}\end{array}\right]$$

(4)

where *E*_{total} is the total (nuclear and electronic) energy and the dynamic metric matrices are [26, 29]

$${\left({\mathbf{C}}_{XY}\right)}_{ik,\phantom{\rule{0.25em}{0ex}}jl}={-2\mathrm{Im}\frac{{\partial}^{2}\mathrm{ln}S}{\partial {X}_{ik}^{\text{'}}\partial {Y}_{jl}}|}_{\begin{array}{l}{\mathbf{R}}^{\prime}=\mathbf{R}\\ \mathbf{P}\text{'}=\mathbf{P}\end{array}};\phantom{\rule{0.5em}{0ex}}{\left({\mathbf{C}}_{{X}_{ik}}\right)}_{ph}={\frac{{\partial}^{2}\mathrm{ln}S}{\partial {z}_{ph}^{*}\partial {X}_{ik}}|}_{\begin{array}{l}{\mathbf{R}}^{\prime}=\mathbf{R}\\ \mathbf{P}\text{'}=\mathbf{P}\end{array}};\phantom{\rule{0.5em}{0ex}}{\mathbf{C}}_{ph,\phantom{\rule{0.25em}{0ex}}qg}={\frac{{\partial}^{2}\mathrm{ln}S}{\partial {z}_{ph}^{*}\partial {z}_{qg}}|}_{\begin{array}{l}{\mathbf{R}}^{\prime}=\mathbf{R}\\ \mathbf{P}\text{'}=\mathbf{P}\end{array}}$$

(5)

where *X*_{ik} and *Y*_{jl} denote either **R**_{A = i,k} or **P**_{A = j,l} and *S* = **z**(*t*), **R**′(*t*), **P**′(*t*)|**z**(*t*), **R**(*t*), **P**(*t*). **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–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 $\left\{{\mathbf{R}}_{A}^{i}\right\}$, momenta $\left\{{\mathbf{P}}_{A}^{i}\right\}$ and Thouless electronic state |**z**^{(i)} = **0**, **R**^{(i)} = |0, i.e. the UHF ground state of the reactants’ super-molecule. When the reaction accesses the non-adiabatic regime, **z**(*t*) ≠ **0** and |**z**(*t*),**R**(*t*) becomes a superposition of ground |0 and excited |… *h* → *p* … UHF states (cf. Eq 2)

$$|{\mathrm{\Psi}}_{e}^{SLEND}\left(t\right)\u27e9=|\mathbf{z}(t);\mathbf{R}\left(t\right),\mathbf{P}\left(t\right)\u27e9={\displaystyle {\sum}_{h=1,\phantom{\rule{0.25em}{0ex}}p={N}_{e}+1}^{K,\phantom{\rule{0.25em}{0ex}}{N}_{e}+1}{z}_{ph}\left(t\right)|\dots \phantom{\rule{0.25em}{0ex}}h\to p\phantom{\rule{0.25em}{0ex}}\dots \u27e9}+\phantom{\rule{0.25em}{0ex}}\dots $$

(6)

All the present SLEND simulations were computed with our END program PACE (Python-Accelerated Coherent states Electron-nuclear dynamics, T. V. Grimes, E. S. Teixeira and J. A. Morales, Texas Tech University, 2010–2016; cf. Ref. [26], Sect. 4). PACE combines several advanced computer science techniques such as a mixed programming language (Python for logic flow and Fortran and C++ for calculations), intra- and internode parallelization, and the fast OED/ERD atomic integral package [49] from the ACES III/IV [50] program. In addition, the water clusters’ geometries at initial time were computed with the NWChem [51] and GAMESS [52] programs.

Present SLEND simulations start with the super-molecular systems H^{+} + (H_{2}O)_{1-6} optimized at the UHF level and having projectile-to-target [H^{+}-to-(H_{2}O)_{1-6}] separations ≥ 30.00 a.u. (cf. Fig 2). Integration of the SLEND Eq 4 requires the evaluation of their various terms at numerous adaptive time steps on a total of 25,020 trajectories. Thus, for feasibility’s sake, the medium size basis sets 6-31G* [53] [for H^{+} + (H_{2}O)_{1-6}] and 6-31G** [53] [for H^{+} + (H_{2}O)_{1-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–63] and experimental [35, 64–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–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–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].

For this investigation, we selected ten representative isomers in the (H_{2}O)_{1-6} series: each single isomer of (H_{2}O)_{1-3,} the two isomers of (H_{2}O)_{4,} two isomers of (H_{2}O)_{5} out of a total of at least four [34]_{,} and three isomers of (H_{2}O)_{6} out of a total of at least twelve [34] (cf. Fig 1). The isomers were calculated at the UHF/6-31G* [(H_{2}O)_{1-6}] and /6-31G** [(H_{2}O)_{1-5}] levels. When two or more isomers of a (H_{2}O)_{n} cluster are considered, they are depicted (cf. Fig 1) and listed (cf. Table 1) in their increasing order of energies [e.g. hexamer a (prism), hexamer b (cage), and hexamer c since ${E}_{\text{Hexamer b}}^{\text{Cage}}-{E}_{\text{Hexamer a}}^{\text{Prism}}$ = 4.2 kJ/mol and ${E}_{\text{Hexamer c}}^{}-{E}_{\text{Hexamer a}}^{\text{Prism}}$ = 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-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.

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 ${\mathbf{R}}_{{\text{H}}^{+}}^{0}=(b\ge 0,\phantom{\rule{0.25em}{0ex}}0,\phantom{\rule{0.5em}{0ex}}+30\phantom{\rule{0.25em}{0ex}}\text{a}.\text{u}.)$ and momentum ${\mathbf{P}}_{{\text{H}}^{+}}^{0}=(0,\phantom{\rule{0.25em}{0ex}}0,-{p}_{{\text{H}}^{+}}^{z})$, where *b* ≥ 0 is the projectile impact parameter measured from the (H_{2}O)_{1-6} centers of mass, and ${p}_{{\text{H}}^{+}}^{z}$ corresponds to *E*_{Lab} = 100 keV (cf. Fig 2, panel I). Having set ${\mathbf{R}}_{{\text{H}}^{+}}^{0}$ and ${\mathbf{P}}_{{\text{H}}^{+}}^{0}$, 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 ${\mathbf{R}}_{{\text{H}}^{+}}^{i}$ and ${\mathbf{P}}_{{\text{H}}^{+}}^{i}$ are therefore obtained by rotating ${\mathbf{R}}_{{\text{H}}^{+}}^{0}$ and ${\mathbf{P}}_{{\text{H}}^{+}}^{0}$ 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}, 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 ${\mathbf{R}}_{{\text{H}}^{+}}^{0}$ and ${\mathbf{P}}_{{\text{H}}^{+}}^{0}$ projections on the x-y plane and *β* is measured directly from the +z axis (cf. Fig 2). These rotations define projectile-target relative orientations Ω_{i} = (*α*_{i}, *β*_{i}, *γ*_{i}). The definite initial conditions of the H^{+} projectile for the simulations, ${\mathbf{R}}_{{\text{H}}^{+}}^{i}$ and ${\mathbf{P}}_{{\text{H}}^{+}}^{i}$, 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: Ω_{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 $\underset{0}{\overset{2\pi}{\int}}{\displaystyle \underset{0}{\overset{\pi}{\int}}{\displaystyle \underset{0}{\overset{2\pi}{\int}}{D}_{{M}^{\prime}M}^{J}}}}\phantom{\rule{0.25em}{0ex}}\left(\alpha ,\phantom{\rule{0.25em}{0ex}}\beta ,\phantom{\rule{0.25em}{0ex}}\gamma \right)\phantom{\rule{0.25em}{0ex}}\mathrm{sin}\beta d\alpha \phantom{\rule{0.25em}{0ex}}d\beta d\gamma =0$ for 2 ≤ *J* ≤ 5). For a given orientation Ω_{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 and run for a total time of 30.00 a.u. (0.7257 fs); this simulation time ensures that the final projectile-target separation is at least equal to the initial one (30.00 a.u.). Simulations for the prediction of fragmentations utilize the grids denoted as [*b*]_{2} in 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.

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]

$${\sigma}_{1-ET}=\frac{1}{4\pi}{\displaystyle \underset{0}{\overset{\infty}{\int}}{\displaystyle \underset{0}{\overset{2\pi}{\int}}{\displaystyle \underset{0}{\overset{\pi}{\int}}{\displaystyle \underset{0}{\overset{2\pi}{\int}}b\phantom{\rule{0.5em}{0ex}}{P}_{1\text{-ET}}\phantom{\rule{0.25em}{0ex}}\left(\alpha ,\phantom{\rule{0.25em}{0ex}}\beta ,\phantom{\rule{0.25em}{0ex}}\gamma ,\phantom{\rule{0.25em}{0ex}}b\right)db\phantom{\rule{0.5em}{0ex}}\mathrm{sin}\beta}}}}\phantom{\rule{0.25em}{0ex}}d\alpha \phantom{\rule{0.25em}{0ex}}d\beta \phantom{\rule{0.25em}{0ex}}d\gamma $$

(7)

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 potentials, *P*_{1−ET} (*α*, *β*, *γ*, *b*) → *P*_{1−ET} (*b*), and Eq (7) simplifies to the familiar expression: ${\sigma}_{1-\text{ET}}=2\pi {\displaystyle \underset{0}{\overset{\infty}{\int}}b}\phantom{\rule{0.5em}{0ex}}{P}_{1-\text{ET}}\left(b\right)\phantom{\rule{0.25em}{0ex}}db$ [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]

$$\begin{array}{l}{P}_{0-\text{ET}}\left(\alpha ,\phantom{\rule{0.25em}{0ex}}\beta ,\phantom{\rule{0.25em}{0ex}}\gamma ,\phantom{\rule{0.25em}{0ex}}b\right)=\left(1-{N}_{\alpha}\right)\left(1-{N}_{\beta}\right);\\ {P}_{1-\text{ET}}\left(\alpha ,\phantom{\rule{0.25em}{0ex}}\beta ,\phantom{\rule{0.25em}{0ex}}\gamma ,\phantom{\rule{0.25em}{0ex}}b\right)={N}_{\alpha}\left(1-{N}_{\beta}\right)+{N}_{\beta}\left(1-{N}_{\alpha}\right);\\ {P}_{2-\text{ET}}\left(\alpha ,\phantom{\rule{0.25em}{0ex}}\beta ,\phantom{\rule{0.25em}{0ex}}\gamma ,\phantom{\rule{0.25em}{0ex}}b\right)={N}_{\alpha}{N}_{\beta}\end{array}$$

(8)

where *N*_{α} and *N*_{β} are the number of *α*- and *β*-spin electrons in the outgoing projectile calculated from their respective electron densities ${\rho}_{\alpha}^{out.\phantom{\rule{0.25em}{0ex}}proj.}\left(\mathbf{r}\right)$ and ${\rho}_{\beta}^{out.\phantom{\rule{0.25em}{0ex}}proj.}\left(\mathbf{r}\right)$:

$${N}_{\alpha}={\displaystyle \int {\rho}_{\alpha}^{out.\phantom{\rule{0.25em}{0ex}}proj.}\left(\mathbf{r}\right)}\phantom{\rule{0.25em}{0ex}}d\mathbf{r};\phantom{\rule{0.5em}{0ex}}{N}_{\beta}={\displaystyle \int {\rho}_{\beta}^{out.\phantom{\rule{0.25em}{0ex}}proj.}\left(\mathbf{r}\right)}\phantom{\rule{0.25em}{0ex}}d\mathbf{r}$$

(9)

Eqs 8 and 9 are evaluated at the final simulation time, when the outgoing projectile and the clusters are well separated by at least 30.00 a.u. of length; therefore, *N*_{α} and *N*_{β} are the number of electrons unequivocally assigned to the distant outgoing projectile. Moreover, at those separations, *N*_{α} and *N*_{β} from Eq 9 and those from any ordinary electron population analyses (Mulliken, Löwdin, etc.) [53] become identical; this assures that *N*_{α} and *N*_{β} are free of any arbitrary criterion for electrons’ distributions over the projectile and clusters. Other calculated properties are the orientation-averaged 1-ET probabilities ${\overline{P}}_{1-\text{ET}}\left(b\right)$, ${\overline{P}}_{1-\text{ET}}\left(b\right)={\left(8{\pi}^{2}\right)}^{-1}{\displaystyle {\int}_{0}^{2\pi}{\displaystyle {\int}_{0}^{\pi}{\displaystyle {\int}_{0}^{2\pi}{P}_{1-\text{ET}}}}}\phantom{\rule{0.25em}{0ex}}\left(\alpha ,\phantom{\rule{0.25em}{0ex}}\beta ,\phantom{\rule{0.25em}{0ex}}\gamma ,\phantom{\rule{0.25em}{0ex}}b\right)\phantom{\rule{0.25em}{0ex}}\mathrm{sin}\beta $, and their *b* – weighted counterparts, $b{\overline{P}}_{1-\text{ET}}\left(b\right)$, where ${\sigma}_{1\text{-ET}}=2\pi {\displaystyle {\int}_{0}^{\infty}b}\phantom{\rule{0.5em}{0ex}}{\overline{P}}_{1-\text{ET}}\phantom{\rule{0.25em}{0ex}}\left(b\right)\phantom{\rule{0.25em}{0ex}}db$. ${\overline{P}}_{1-\text{ET}}\left(b\right)$ and $b{\overline{P}}_{1-\text{ET}}\left(b\right)$ reveal more mechanistic details than *σ*_{1-ET}.

The first property calculated in this investigation is the cluster-to-proton total 1-ET ICS, *σ*_{1−ET}, for the H^{+} + (H_{2}O)_{n} systems, *n* = 1–6, at *E*_{Lab} = 100 keV at the SLEND/6-31G* (*n* = 1–6) and SLEND/6-31G** (*n* = 1–5) levels. Both SLEND *σ*_{1−ET} for the monomeric system (*n* = 1) are listed in Table 2 along with their available counterparts from four experiments denoted as Exp. A to D [70–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, ${\overline{\sigma}}_{1-ET}^{Exp.}$, and its average relative error, *e*_{Exp.}, are 1.27 Å^{2} and ± 10.62%, respectively. The theoretical ${\sigma}_{1-ET}^{Theo.}$ ‘s and their average relative deviations ${\overline{\mathrm{\Delta}}}_{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 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 (${\overline{\mathrm{\Delta}}}_{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 (${\overline{\mathrm{\Delta}}}_{Theo.}$ = 63.0%). The reason and remediation of the SLEND/6-31** *σ*_{1−ET} overestimation will be discussed in detail in the Further Analysis Sub-section. 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\phantom{\rule{0.25em}{0ex}}\left(n\right)=\left(4/3\right)\phantom{\rule{0.25em}{0ex}}\pi \phantom{\rule{0.5em}{0ex}}{R}_{n}^{3}$ and $A\left(n\right)=4\pi \phantom{\rule{0.5em}{0ex}}{R}_{n}^{2}$, 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 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., ${\sigma}_{1-ET}^{6-31{G}^{*}}$ = 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 Tables22 and and33 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 *σ*_{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–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 Figs44 and and55 show the orientation-averaged 1-ET probabilities ${\overline{P}}_{1-ET}^{}(b)$ and their *b* – weighted counterparts $b{\overline{P}}_{1-ET}^{}(b)$ vs. *b* for the considered (H_{2}O)_{1-6} isomers. With both basis sets, the ${\overline{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 ${\overline{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 $b{\overline{P}}_{1-ET}^{}(b)$ vs. *b* plots exhibit similar patterns to those of the ${\overline{P}}_{1-ET}^{}(b)$ but modulated by the *b* value. As the integrands of the *σ*_{1−ET}, the $b{\overline{P}}_{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.).

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 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*) is a superposition of various UHF states corresponding to various products’ channels [69, 75], e.g. ${\text{H}}_{proj.}^{+}+{\text{H}}_{2}\text{O}\to {\text{H}}_{proj.}^{{q}_{1}=+1}+{\text{H}}_{2}^{{q}_{2}=0}+{\text{O}}_{}^{{q}_{3}=0}$ or $\to {\text{H}}_{proj.}^{{q}_{1}=0}+{\text{H}}_{2}^{{q}_{2}=+1}+{\text{O}}_{}^{{q}_{3}=0}$ or $\to {\text{H}}_{proj.}^{{q}_{1}=-1}+{\text{H}}_{2}^{{q}_{2}=+1}+{\text{O}}_{}^{{q}_{3}=+1}$, etc., where ${\text{H}}_{proj.}^{{q}_{1}=0,+1,-1}$ 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 ${\text{H}}_{proj.}^{{q}_{1}}$) are not necessarily integer numbers corresponding to canonical chemical species (e.g. ${\text{H}}_{proj}^{{q}_{1}=+1}$, ${\text{H}}_{proj}^{{q}_{1}=0}$ and ${\text{H}}_{proj}^{{q}_{1}=-1}$) 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. ${\text{H}}_{proj.}^{+1}+{\left({\text{H}}_{2}\text{O}\right)}_{n}\to {\text{H}}_{proj.}^{{q}_{1}}+{\left({\text{H}}_{2}\text{O}\right)}_{n}^{1-{q}_{1}}$, 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 ${\text{H}}_{proj.}^{{q}_{1}}$ and clusters fragments A^{q2} 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–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 the number of total simulations by using a finer grid. Such improvements entail further computational cost and will be attempted later.

The predicted fragmentations are:

**H**^{+}
**+ H**_{2}**O simulations:** 9 out of 252 simulations exhibited the H_{2}O target fragmenting into: H + OH (2 simulations), H + H + O (6 simulations) and H_{2} + O (1 simulation, cf. Fig 6).

**H**^{+}
**+ (H**_{2}**O)**_{2}
**simulations:** 9 out of 540 simulations exhibited the (H_{2}O)_{2} target fragmenting into: H_{2}O + HO + H (2 simulations), H_{3}O + O + H (1 simulation), H_{2}O + 2H + O (5 simulations) and H_{3}O + OH (1 simulation, cf. Fig 7).

**H**^{+}
**+ (H**_{2}**O)**_{3}
**simulations:** 3 out of 540 simulations exhibited the (H_{2}O)_{3} target fragmenting into: H_{3}O + OH + H_{2}O (2 simulations) and H + OH + 2H_{2}O (1 simulation).

**H**^{+}
**+ (H**_{2}**O)**_{4}
**(tetramer a) simulations:** 10 out of 748 simulations exhibited the (H_{2}O)_{4} target fragmenting into: (H_{2}O)_{2} + 2H_{2}O (4 simulations), H_{3}O + OH + 2H_{2}O (4 simulations), an H + OH + 3H_{2}O (2 simulations)

In conclusion, present SLEND/6-31G* simulations predict the DNA-damaging radicals H, OH, O and H_{3}O, and the innocuous species H_{2}O and (H_{2}O)_{2} as water radiolysis products.

To illustrate the predicted fragmentations, we present a few animation stills from some representative simulations. Fig 6 shows five sequential stills of the animation of ${\text{H}}_{proj.}^{+}+{\text{H}}_{2}\text{O}\to {\text{H}}_{proj.}^{+}+{\text{H}}_{2}+\text{O}$ and a sixth panel plotting the Mulliken populations of the H_{proj.}, H_{2} and O moieties vs. time. This last panel permits the observation of the time evolution of the electrons. Mulliken populations are basis-set dependent and, more importantly, somewhat arbitrary in the way they distribute electrons over neighboring fragments. However, in previous SLEND 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). The incoming projectile ${\text{H}}_{proj.}^{+}\phantom{\rule{0.25em}{0ex}}$ passes in between the H atoms of H_{2}O, hits the O atom and bounces back. As a result of this collision, the H_{2} and O moieties of H_{2}O break apart; the ejected H_{2} moiety undergoes a series of strong oscillations ranging from the near dissociation of H_{2} into H atoms to these atoms’ recombination back to H_{2}. Finally, Fig 7 shows the nine predicted fragments in H^{+} + (H_{2}O)_{2}.

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 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 bound-to-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 $|{\mathrm{\Psi}}_{e}^{SLEND}\left({t}_{f}\right)\u27e9$, Eq 6, was projected onto the ground |0 and excited states |… *h* → *p* … 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: ${\sigma}_{DI}^{SLEND}$ = 13.9Å and ${\sigma}_{DI}^{EXPT.}$ = 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** ${\sigma}_{CT}^{SLEND}$ via Eq 8 can be corrected by subtracting the ${\sigma}_{DI}^{SLEND}$ part from it: ${\sigma}_{CT}^{SLEND}$ = 2.06 Å2–0.92 Å^{2} = 1.14 Å^{2}; this places SLEND/6-31G** ${\sigma}_{CT}^{SLEND}$ within the range of the four experimental values and closest to that from Exp. D [73]: ${\sigma}_{CT}^{Expt.D}$ = 1.13 Å Δ_{Theo.} = 0.8%. A similar correction might occur with the SLEND/6-31G* ${\sigma}_{CT}^{SLEND}$ 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* … from |0 forming a full CI expansion should be generated so that $|{\mathrm{\Psi}}_{e}^{SLEND}\left({t}_{f}\right)\u27e9$ in Eq 6 can be projected on all those states. This more demanding capability is not currently available in PACE but is under development.

To model microscopic processes in PCT [1–7], the SLEND method was applied to the H^{+} + (H_{2}O)_{n} systems at *E*_{Lab} = 100 keV with the 6-31G* (*n* = 1–6) and 6-31G** (*n* = 1–5) basis sets. Ten (H_{2}O)_{1–6} clusters were selected for this study: eight exhibit mono- and di-cyclic quasi-planar/multiplanar structures [33–35] and two others, the prism and cage (H_{2}O)_{6} isomers, exhibit multi-cyclic, three-dimensional, drop-like structures [34, 35]. These “smallest-drop” clusters were purposely included in this study in an attempt to reproduce early manifestations of bulk-water properties in PCT. Short-time SLEND/6-31G* (*n* = 1–6) and /6-31G** (*n* = 1–5) simulations render cluster-to-projectile total 1-ET ICS, *σ*_{1−ET}, and 1-ET probabilities. In absolute quantitative terms, SLEND/6-31G**σ*_{1−ET} compares satisfactorily with alternative experimental [70–73] and theoretical[22, 23] results only available for *n* = 1, and exhibits almost the same accuracy of the best alternative theoretical result from BGM [23] calculations. SLEND/6-31G** overestimates *σ*_{1−ET} and a detail account about the cause and remediation of this effect was presented. In relative quantitative terms, both SLEND/6-31* and /6-31G** *σ*_{1−ET} precisely fit into physically justified scaling formulae *σ*_{1−ET} (*n*) = *cn*^{2/3} as a function of the cluster size *n*. Long-time SLEND/6-31G* (*n* = 1–4) simulations predict the formation of the DNA-damaging radicals H, OH, O and H_{3}O. While “smallest-drop” isomers were included, no incipient manifestations of bulk-water PCT properties are observed. Therefore, to capture bulk-water effects, simulations with larger water clusters are currently underway.

The authors thank Prof. Takao Oi (Sophia University, Japan) for providing the optimized geometry of the prism isomer of (H_{2}O)_{6} [34]. Present calculations were performed at the Texas Tech University High Performance Computer Center and the Texas Advanced Computing Center at the University of Texas at Austin.

This work was supported by 1) Science without Borders Fellowship Program of the National Council for Scientific and Technological Development of Brazil, to EST, http://www.cienciasemfronteiras.gov.br/web/csf-eng/faq: 2) Cancer Prevention and Research Institute of Texas to JAM, grant: RP140478 http://www.cprit.state.tx.us/. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Data Availability

All relevant data are within the paper.

1. Levin WP, Kooy H, Loeffler JS, DeLaney TF. Proton beam therapy. Br J Cancer. 2005;93:849–54. doi: 10.1038/sj.bjc.6602754
[PMC free article] [PubMed]

2. Solov’yov AV, Surdutovich E, Scifoni E, Mishustin I, Greiner W. Physics of ion beam cancer therapy: A multiscale approach. Phys Rev E. 2009;79(1):011909. [PubMed]

3. Surdutovich E, Yakubovich A. V., Solov'yov AV. Multiscale approach to radiation damage induced by ion beams: complex DNA damage and effects of thermal spikes. European Physical Journal D. 2010;60:101–8.

4. Surdutovich E, Solov'yov AV. Multiscale approach to the physics of radiation damage with ions. European Physical Journal D. 2014;68:353.

5. Girdhani S, Sachs R, Hltaky L. Biological Effects of Proton Radiation: What We Know and Don't Know. Radiat Res. 2013;179:257–72. doi: 10.1667/RR2839.1
[PubMed]

6. McLaurin PM, Privett A, Stopera C, Grimes TV, Perera A, Morales JA. In Honor of N. Yngve Öhrn: Surveying Proton Cancer Therapy Reactions with Öhrn’s Electron Nuclear Dynamics Method. Aqueous Clusters Radiolysis and DNA-Bases Damage by Proton Collisions. Mol Phys. 2015;113:297–313.

7. Jäkel O. Medical Physics Aspects of Particle Therapy. Radiat Prot Dosim. 2009;137:156–66. [PubMed]

8. Mayani DD. Proton therapy for cancer treatment. Journal of Oncology Pharmacy Pratice
2010;17(186–190). [PubMed]

9. Pichl L, Kimura M, Li Y, Buenker RJ. Branching ratios for secondary processes of water ions induced by proton beams in radiation therapy of cancer. IEEE T Nucl Sci. 2004;51(4):1407–11.

10. Chang H, Oehri W, Elsner P, Thiele JJ. The Role of H2O2 as a Mediator of UVB-induced Apoptosis in Keratinocytes. Free Radical Res. 2003;37:655–63. [PubMed]

11. Shikazono N, Noguchi M, Fujii K, Urushibara A, Yokoya A. The Yield, Processing, and Biological Consequences of Clustered DNA Damage Induced by Ionizing Radiation. J Radiat Res. 2009;50:27–36. [PubMed]

12. Friedland W, Jacob P, Bernhardt P, Paretzke HG, Dingfelder M. Simulation of DNA Damage after Proton Irradiation. Radiat Res. 2003;159(3):401–10. [PubMed]

13. Watanabe R, Saito K. Monte carlo simulation of strand-break induction on plasmid DNA in aqueous solution by monoenergetic electrons. Radiat Environ Bioph. 2002;41:207–15. [PubMed]

14. Agostinelli S, Allison J, Amako K, Apostolakis J, Araujo H, Arce P, et al.
Geant4—a simulation toolkit. Nucl Instrum Meth A. 2003;506(3):250–303.

15. Illescas C, Errea LF, Mendez L, Pons B, Rabadan I. Ion-water collisions at intermediate energies. XXVII International Conference on Photonic, Electronic and Atomic Collisions (ICPEAC 2011)
Journal of Physics: Conference Series. 2012;388:102007.

16. Nikjoo H, Uehara S, Emfietzoglou D, Cucinotta FA. Track-structure codes in radiation research,. Radiat Meas. 2006;41:1052–74.

17. Incerti S, Ivanchenko A, Karamitros M, Mantero A, Moretto P, Tran HN, et al.
Comparison of GEANT4 very low energy cross section models with experimental data in water. Med Phys. 2010;37(9):4692–708. doi: 10.1118/1.3476457
[PubMed]

18. Fassò A., Ferrari A, Sala PR. Radiation transport calculations and simulations. Radiat Prot Dosim. 2009;137:118–33 [PubMed]

19. Garcia-Molina R, Abril I, Denton CD, Heredia-Avalos S, Kyriakou I, Emfietzoglou D. Calculated depth-dose distributions for H+ and He+ beams in liquid water,. Nucl Instrum Meth B. 2009;267:2647–52.

20. Turner JE, Hamm RN, Wright HA, Ritchie RH, Magee JL, Chatterjee A, et al.
Studies to link the basic radiation physics and chemistry of liquid water,. Radiat Phys Chem. 1988;32:503–10.

21. Champion C, Weck PF, Lekadir H, Galassi ME, Fojón OA, Abufager P, et al.
Proton-induced single electron capture on DNA/RNA bases. Phys Med Bio. 2012;57:3039–49. [PubMed]

22. Champion C, Quinto MA, Monti JM, Galassi ME, Weck PF, Fojón OA, et al.
Water versus DNA: new insights into proton track-structure modelling in radiobiology and radiotherapy. Phys Med Bio. 2015;60:7805–28. [PubMed]

23. Murakami M, Kirchner T, Horbatsch M, Lüdde HJ. Single and multiple electron removal processes in proton–water-molecule collisions. Phys Rev A. 2012;85(5):052704.

24. Cabrera-Trujillo R, Sabin JR, Deumens E, Öhrn Y. Orientation Effects in Energy Deposition by Protons in Water. In: Sabin JR, Brandas EJ, editors. Adv Quantum Chem. 482005. p. 47–57.

25. Quinet O, Deumens E, Öhrn Y. Proton collisions with the water dimer at keV energies. Int J Quantum Chem. 2009;109(2):259–65.

26. Stopera C, Grimes TV, McLaurin PM, Privett A, Morales JA. Some Recent Developments in the Simplest-Level Electron Nuclear Dynamics Method: Theory, Code Implementation, and Applications to Chemical Dynamics. In: Sabin JR, Brandas EJ, editors. Adv Quantum Chem. Vol. 66, Chapter 3: Elsevier; 2013. p. 113–49.

27. Privett A, Morales JA. Electron Nuclear Dynamics of Proton Collisions with DNA/RNA Bases at ELab = 80 keV: A Contribution to Proton Cancer Therapy Research. Chem Phys Lett. 2014; 603:82–8.

28. Hedström M, Morales JA, Deumens E, Öhrn Y. Electron nuclear dynamics of H+ + H2O collisions. Chem Phys Lett. 1997;279(3–4):241–6.

29. Deumens E, Diz A, Longo R, Öhrn Y. Time-dependent theoretical treatments of the dynamics of electrons and nuclei in molecular systems. Rev Mod Phys. 1994;66(3):917–83.

30. Perera SA, McLaurin PM, Grimes TV, Morales JA. Time-dependent density-functional theory method in the electron nuclear dynamics framework. Chem Phys Lett. 2010;496(1–3):188–95.

31. Thouless DJ. Stability conditions and nuclear rotations in the Hartree-Fock theory. Nucl Phys. 1960;21:225–32.

32. Tabet J, Eden S, Feil S, Abdoul-Carime H, Farizon B, Farizon M, et al.
Absolute total and partial cross sections for ionization of nucleobases by proton impact in the Bragg peak velocity range. Phys Rev A. 2010;82(2):022703.

33. Xantheas SS, Dunning TH. Ab initio studies of cyclic water clusters (H2O)n, n = 1–6. i. optimal structures and vibrational spectra. J Chem Phys. 1993;99:8774–92.

34. Oi T. Molecular Orbital Estimation of Reduced Partition Function Ratios of Small Water Clusters. J Nucl Sci Technol. 2002;39:419–25.

35. Saykally RJ, Wales DJ. Pinning Down the Water Hexamer. Science. 2012;336:814–5. doi: 10.1126/science.1222007
[PubMed]

36. Jungwirth P. How Many Waters Are Necessary To Dissolve a Rock Salt Molecule?
J Phys Chem A. 2000;104(1):145–8.

37. Yoshikawa A, Morales JA. The onset of dissociation in the aqueous LiOH clusters: a solvation study with the effective fragment potential model and quantum mechanics methods. J Mol Struct-THEOCHEM. 2004;681(1–3):27–40.

38. Surdutovich E, Solov'yov AV. Shock wave initiated by an ion passing through liquid water. Physical Review E. 2010;82:051915. [PubMed]

39. Hagelberg F. Electron Dynamics in Molecular Interactions: Principles and Applications
1ed.
Singapure: World Scientific; 2014.

40. Kramer P, Saraceno M. Geometry of the time-dependent variational principle in quantum mechanics
New York, New York: Springer; 1981 1981.

41. Delos JB. Thoery of electronic transitions in slow collisions. Review of Modern Physics. 1981;53(2):287–358.

42. Child MS. Molecular Collision Theory. 1 ed.
New York: Academic Press, Inc.; 1974.

43. Longo R, Diz A, Deumens E, Öhrn Y. Influence of Electronic Nuclear Coupling on Dynamics. Chem Phys Lett. 1994;220(3–5):305–11.

44. Coutinho-Neto M, Deumens E, Öhrn Y. Abstraction and exchange mechanisms for the D-2+NH3+ reaction at hyperthermal collision energies. J Chem Phys. 2002;116(7):2794–802.

45. Longo R, Deumens E, Ohrn Y. H+ +H, He, and H-2 Scattering Using a New Time-Dependent Method for Electron-Nuclear Dynamics. J Chem Phys. 1993;99(6):4554–65.

46. Cabrera-Trujillo R, Sabin JR, Ohrn Y, Deumens E. Direct differential-cross-section calculations for ion-atom and atom-atom collisions in the keV range. Phys Rev A. 2000;61(3):032719.

47. Cabrera-Trujillo R, Ohrn Y, Deumens E, Sabin JR. Stopping cross section in the low- to intermediate-energy range: Study of proton and hydrogen atom collisions with atomic N, O, and F. Phys Rev A. 2000;62(5):art. no.-052714.

48. Cabrera-Trujillo R, Deumens E, Ohrn Y, Sabin JR. Impact parameter dependence of electronic and nuclear energy loss of swift ions: H+ -> He and H+ -> H. Nuclear Instruments & Methods in Physics Research Section B-Beam Interactions with Materials and Atoms. 2000;168(4):484–92.

49. Flocke N, Lotrich V. Efficient electronic integrals and their generalized derivatives for object oriented implementations of electronic structure calculations. J Comput Chem. 2008;29(16):2722–36. doi: 10.1002/jcc.21018
[PubMed]

50. Lotrich V, Flocke N, Ponton M, Yau AD, Perera A, Deumens E, et al.
Parallel implementation of electronic structure energy, gradient, and Hessian calculations. J Chem Phys. 2008;128(19):194104
doi: 10.1063/1.2920482
[PubMed]

51. Valiev M, Bylaska EJ, Govind N, Kowalski K, Straatsma TP, Dam HJJV, et al.
NWChem: A comprehensive and scalable open-source solution for large scale molecular simulations. Comp Phys Comm. 2010;181(9):1477–89.

52. Schmidt MW, Baldridge KK, Boatz JA, Elbert ST, Gordon MS, Jensen JH, et al.
General atomic and molecular electronic structure system. J Comput Chem. 1993;14(11):1347–63.

53. Szabo A, Ostlund NS. Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory. 1st, rev ed.
Mineola, NY: Dover Publications Inc; 1989.

54. Liu X, Lu W-C, Wang CZ, Ho KM. Energetic and fragmentation stability of water clusters (H20)n, n = 2–30. Chem Phys Lett. 2011;508:270–5.

55. Hincapié G, Acelas N, Castaño M, David J, Restrepo A. Structural Studies of the Water Hexamer. J Phys Chem A. 2010;114:7809–14. doi: 10.1021/jp103683m
[PubMed]

56. Wang Y, Babin V, Bowman JM, Paesani F. The Water Hexamer: Cage, Prism, or Both. Full Dimensional Quantum Simulations Say Both. J Am Chem Soc. 2012;134
11116–9. doi: 10.1021/ja304528m
[PubMed]

57. Shields RM, Temelso B, Archer KA, Morrell TE, Shields GC. Accurate Predictions of Water Cluster Formation, (H2O)n = 2–10. J Phys Chem A. 2010;114
11725–37. doi: 10.1021/jp104865w
[PubMed]

58. Bates DM, Tschumper GS. CCSD(T) Complete Basis Set Limit Relative Energies for Low-Lying Water Hexamer Structures. J Phys Chem A
2009;113:3555–9. doi: 10.1021/jp8105919
[PubMed]

59. Dahlke EE, Olson RM, Leverentz HR, Truhlar DG. Assessment of hte Accuracy of Density Functionals for Prediction of Relative Energies and Geometries of Low-Lying Isomers of Water Hexamers. J Phys Chem A. 2008;112:3976–84. doi: 10.1021/jp077376k
[PubMed]

60. Tainter CJ, Skinner JL. The water hexamer: Three-body interactions, structures, energetics, and OH-stretch spectroscopy at finite temperature. J Chem Phys. 2012;137:104304
doi: 10.1063/1.4746157
[PubMed]

61. Lee HM, Suh SB, Lee JY, Tarakeshwar P, Kim KS. Structures, energies, vibrational spectra, and electronic properties of water monomer to decamer. J Chem Phys. 2000;112:9759.

62. Losada M, Leutwyler S. Water hexamer clusters: Structures, energies, and predicted mid-infrared spectra. J Chem Phys. 2002;117:2003.

63. Gillan MJ, Alfè D, Michaelides A. Perspective: Hoow good is DFT for water?
J Chem Phys. 2016;144:130901
doi: 10.1063/1.4944633
[PubMed]

64. Keutsch FN, Cruzan JD, Saykally RJ. The Water Trimer. Chem Rev. 2003;103:2533–78. doi: 10.1021/cr980125a
[PubMed]

65. Keutsch FN, Saykally RJ. Water clusters: Untangling the mysteries of the liquid, one molecule at a time,. P Natl Acad Sci USA. 2001;98
10533–40. [PubMed]

66. Pérez C, Muckle MT, Zaleski DP, Seifert NA, Temelso B, Shields GC, et al.
Structures of Cage, Prism, and Book Isomers of Water Hexamer from Broadband Rotational Spectroscopy. Science. 2012;336:897–901. doi: 10.1126/science.1220574
[PubMed]

67. Zare RN. Angular Momentum, Understanding Spatial Aspects in Chemisty and Physics. New York: Wilesy-Interscience Publication; 1988.

68. Mamone S, Pileio G, Levitt MH. Orientational Sampling Schemes Based on Four Dimensional Polytopes. Symmetry. 2010;2(3):1423–49.

69. McLaurin PM. Chemistry PhD Dissertation: New Applications of the Electron Nuclear Dynamics Theory to Scattering Processes and Chemical Reactions: Tool Development, Method Validation, and Computer Simulation. PhD Dissertation, Texas Tech University, Department of Chemistry and Biochemistry. 2011.

70. Rudd ME, Goffe TV, DuBois RD, Toburen LH. Cross sections for ionization of water vapor by 7–4000-keV protons. Phys Rev A. 1985;31(1):492–4. [PubMed]

71. Gobet F, Farizon B, Farizon M, Gaillard MJ, Carré M, Lezius M, et al.
Total, Partial, and Electron-Capture Cross Sections for Ionization of Water Vapor by 20–150 keV Protons. Phys Rev Lett. 2001;86(17):3751–4. doi: 10.1103/PhysRevLett.86.3751
[PubMed]

72. Gobet F, Eden S, Coupier B, Tabet J, Farizon B, Farizon M, et al.
Ionization of water by (20–150)-keV protons: Separation of direct-ionization and electron-capture processes. Phys Rev A. 2004;70(6):062716.

73. Luna H, de Barros ALF, Wyer JA, Scully SWJ, Lecointre J, Garcia PMY, et al.
Publisher’s Note: Water-molecule dissociation by proton and hydrogen impact [Phys. Rev. A 75, 042711 (2007)]. Phys Rev A. 2007;75(5): 059902

74. Baccarelli I, Gianturco FA, Grandi A, Lucchese RR, Sanna N. Electron-driven molecular processes induced in biological systems by eletromagnetic and other ionizing sources. In: Sabin JR, Brandas EJ, editors. Adv Quantum Chem. 52: Elsevier; 2007. p. 189–230.

75. Morales J, Diz A, Deumens E, Öhrn Y. Electron nuclear dynamics of H++H2 collisions at Elab = 30 eV. The Journal of Chemical Physics. 1995;103(23):9968–80.

76. Gobet F, Eden S, Coupier B, Tabet J, Farizon B, Farizon M, et al.
Electron-loss and target ionization cross sections for water vapor by 20–150 keV neutral atomic hydrogen impact. Chem Phys Lett. 2004;421:68–71.

77. Luna H, Barros ALF, Wyer JA, Scully SWJ, Lecointre J, Garcia PMY, et al.
Water-molecule dissociation by proton and hydrogen impact. Physical Review A. 2007;75:042711.

78. Werner U, Beckord K, Becker J, Lutz HO. 3D Imaging of the Collision-Induced Coulomb Fragmentation of Water Molecules. Phys Rev Lett. 1995;74:1962
doi: 10.1103/PhysRevLett.74.1962
[PubMed]

79. Jacquemin D, Morales JA, Deumens E, Ohrn Y. Electron nuclear dynamics of proton collisions with methane at 30 eV. J Chem Phys. 1997;107(16):6146–55.

80. Maiti B, Sadeghi R, Austin A, Morales JA. Coherent-states dynamics of the H^{+} + HF reaction at E_{Lab} = 30 eV: A complete electron nuclear dynamics investigation. Chem Phys. 2007;340(1–3):105–19.

81. Shah MB, Elliott DS, Gilbody HB. Ionisation of atomic hydrogen by 9–75 keV protons. J Phys B: At, Mol Opt Phys. 1987;20:2481.

Articles from PLoS ONE are provided here courtesy of **Public Library of Science**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |