PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of plosonePLoS OneView this ArticleSubmit to PLoSGet E-mail AlertsContact UsPublic Library of Science (PLoS)
 
PLoS One. 2017; 12(4): e0174456.
Published online 2017 April 4. doi:  10.1371/journal.pone.0174456
PMCID: PMC5380356

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

Ilia Solov'yov, Editor

Abstract

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.

Introduction and background

Proton cancer therapy (PCT) is an approved type of radiotherapy that utilizes high-energy H+ projectiles to fight cancer [17]. The ultimate effect of this radiation is to damage the DNA of cancerous cells [17]. If left unrepaired, this damage produces various anomalies in cancerous cells that eventually lead to their death (apoptosis) [17]. 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 [17]. 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 [17].

In all radiotherapies, the radiation predominantly interacts with cellular H2O 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+ + H2O → H+ + H·+ OH·)[24, 6, 9], (b) secondary ions (e.g. H+ + H2O → 2H+ + OH-)[24, 6, 9], (c) reactive molecules (e.g. H+ + 2H2O → H+ + H2 + H2O2) [10], and (d) localized heat in the medium [24, 6, 9]. These species can react with H2O 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 [25, 11].

PCT comprises various processes spanning different space (l = 10−10–100 m) and time (t = 10−21–105 s) scales [24]. For instance, water radiolysis reactions, DNA damage at the genome level [12] and tumor remissions lie at the microscopic (10−10l≤10−8 m), mesoscopic (10−8l≤10−3 m) and macroscopic (10−3l≤100 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 [1214]. 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, 1220], and the latter calculate proper radiation doses for treatments [18].

Although PCT is clinically approved, various PCT details at the microscopic scale remain uncertain [17]. Knowledge of those details is essential for a rational design of PCT seeking to maximize its therapeutic power and minimize its side effects [17, 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 [2123] and time-dependent dynamics methods [6, 2427] have been applied to computationally feasible prototypes of PCT reactions (e.g. H+ + (H2O)1–4 to model water radiolysis reactions [6, 2126, 28] and H+ + DNA/RNA bases to model DNA proton damage [6, 21, 22, 26, 27]).

Among quantum-mechanics methods for PCT [6, 2427], 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+ + (H2O)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, 2126, 28], we present herein an exploratory SLEND study of PCT water radiolysis reactions with the H+ + (H2O)1-6 prototypes at ELab = 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+ + H2O [2124, 28], whose “cluster” is the farthest from being bulk water. While these H+ + H2O studies were indeed useful for investigating radiolysis processes [2124, 28], they could not completely capture the processes occurring in bulk water. For that reason, previous SLEND studies explored H+ + (H2O)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+ + (H2O)1-6 prototypes that include ten isomers in a larger series of clusters (H2O)1-6 (cf. Fig 1) and involve a larger number of proton-cluster orientations (60) and simulations (25,020).

Fig 1
UHF/6-31G** [(H2O)2-5] and UHF/6-31G*[(H2O)6] optimized geometries of the water cluster isomers selected for this study.

The selected (H2O)1-6 series contains the initial terms of the long progression from molecular H2O to bulk water. Specifically, the first terms in this progression, (H2O)1-5, have mono- and di-cyclic quasi-planar/multiplanar structures [3335] exhibiting no “bulky” shapes (cf. Fig 1), but two isomers in the last term—the prism and cage isomers of (H2O)6—have multi-cyclic, three-dimensional structures exhibiting drop-like shapes [34, 35] (cf. Fig 1). In fact, these two (H2O)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+ + (H2O)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 [1520] and radiation dosages [1214, 1620] 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 [24] 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 |ΨTotalEND=|ΨNEND|ΨeEND, which consists of nuclear |ΨNEND and electronic |ΨeEND wavefunctions, and treats |ΨTotalEND 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 |ΨNSLEND for a system having NN nuclei is the product of 3NN frozen, narrow Gaussian wave packets:

X|ΨNSLEND(t)=X|R(t),P(t)=A=13NNexp{[XARA(t)2ΔRA]2+iPA(t)[XARA(t)]}
(1)

with average positions RA(t), average momenta PA(t) and widths ΔRA. To lower computational cost, the zero-width limit, ΔRA → 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 |ΨeSLEND is a spin-unrestricted, single-determinantal wavefunction in the Thouless representation [31]. Specifically, taking Ne occupied {ψh}, and KNe virtual {ψp} molecular spin-orbitals (MSOs), |ΨeSLEND is [31]:

|ΨeSLEND=|z(t);R(t),P(t)=det{χh[xh;z(t),R(t),P(t)]}=exp[h=1,p=Ne+1Ne,Kzph(t)bpbh]|0;
(2)

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

χh[x;z(t),R(t),P(t)]=ψh[x;R(t),P(t)]+p=Ne+1Kzph(t)ψp[x;R(t),P(t)];(1hNe)
(3)

with complex-valued molecular coefficients {zph(t)}. The MSOs are obtained at initial time at the UHF level. The MSOs are constructed with travelling atomic basis set functions ϕAi(riRAPA)—i.e., Slater-type orbitals in terms of contracted Gaussian-type orbitals on the moving nuclear centers RA(t) and augmented with electron translation factors (ETFs) [41] to include explicit nuclear momenta PA(t) effects. MSOs and DSOs are spin-unrestricted; therefore, the unrestricted determinant |z(t);R(t),P(t)[right angle bracket] 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] LSLEND=ΨTotalSLEND|(i/2)(/t/t)H^Total|ΨTotalSLEND/ΨTotalSLEND|ΨTotalSLEND1 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] ASLEND, δASLEND=δt1t2LSLEND(t)dt=0. The resulting SLEND dynamical equations are [26, 29]:

[iC0iCRiCP0iC*iCR*iCP*iCRiCRTCRRI+CRPiCPiCPTI+CPRCRP][dzdtdz*dtdRdtdPdt]=[ETotalz*ETotalzETotalRETotalP]
(4)

where Etotal is the total (nuclear and electronic) energy and the dynamic metric matrices are [26, 29]

(CXY)ik,jl=2Im2lnSXik'Yjl|R=RP'=P;(CXik)ph=2lnSzph*Xik|R=RP'=P;Cph,qg=2lnSzph*zqg|R=RP'=P
(5)

where Xik and Yjl denote either RA = i,k or PA = j,l and S = [left angle bracket]z(t), R′(t), P′(t)|z(t), R(t), P(t)[right angle bracket]. CR and CRR 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 CP = CPR = 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[4548] 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 {RAi}, momenta {PAi} and Thouless electronic state |z(i) = 0, R(i)[right angle bracket] = |0[right angle bracket], 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)[right angle bracket] becomes a superposition of ground |0[right angle bracket] and excited |… hp[right angle bracket] UHF states (cf. Eq 2)

|ΨeSLEND(t)=|z(t);R(t),P(t)=h=1,p=Ne+1K,Ne+1zph(t)|hp+
(6)

Computational details

Software

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.

Water cluster structures at the initial states

Present SLEND simulations start with the super-molecular systems H+ + (H2O)1-6 optimized at the UHF level and having projectile-to-target [H+-to-(H2O)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+ + (H2O)1-6] and 6-31G** [53] [for H+ + (H2O)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).

Fig 2
H+ + (H2O)1-6 initial conditions.

Numerous theoretical [33, 34, 5463] and experimental [35, 6466] 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 (H2O)2-3 clusters present one isomer each [34, 55]; however, the (H2O)n, n ≥ 4, clusters present a variable number of isomers that rapidly increase with the cluster size n [34, 55]. (H2O)1-5 present mono- and di-cyclic quasi-planar/multiplanar structures [3335], whereas (H2O)n, n ≥ 6, present multi-cyclic, three-dimensional structures in addition to quasi-planar/multiplanar ones [34, 35] (cf. Fig 1). The three-dimensional (H2O)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 (H2O)1-5 isomers, but discrepancies arise in the energy orders of the (H2O)n, n ≥ 6, isomers [35, 56]. Recent spectroscopy experiments have identified the cage isomer as the lowest-energy (H2O)6 structure [35, 66]. However, ab initio calculations at various levels of accuracy have disagreed on whether the prism [34, 5759], the cage [60, 61], both of them [56], or the chair [62] isomer(s) is(are) the lowest-energy (H2O)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 (H2O)6 structure [58].

For this investigation, we selected ten representative isomers in the (H2O)1-6 series: each single isomer of (H2O)1-3, the two isomers of (H2O)4, two isomers of (H2O)5 out of a total of at least four [34], and three isomers of (H2O)6 out of a total of at least twelve [34] (cf. Fig 1). The isomers were calculated at the UHF/6-31G* [(H2O)1-6] and /6-31G** [(H2O)1-5] levels. When two or more isomers of a (H2O)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 EHexamer bCageEHexamer aPrism = 4.2 kJ/mol and EHexamer cEHexamer aPrism = 15.4 kJ/mol with UHF/6-31G*). Notice that the UHF prism isomer is the lowest-energy (H2O)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 (H2O)1-6 are not intended to contribute to the resolution of the discrepancies about the (H2O)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.

Table 1
Grids for the projectile impact parameter b for the SLEND 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+ + (H2O)1-6 are assembled for the initial conditions (cf. Fig 2). The (H2O)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 (H2O)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 RH+0=(b0,0,+30a.u.) and momentum PH+0=(0,0,pH+z), where b ≥ 0 is the projectile impact parameter measured from the (H2O)1-6 centers of mass, and pH+z corresponds to ELab = 100 keV (cf. Fig 2, panel I). Having set RH+0 and PH+0, various projectile-target relative orientations can be generated by rotating a (H2O)1-6 target according to ordinary Euler angles prescriptions [67]. However, such a procedure involves the electronic re-optimization of the (H2O)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 (H2O)1-6 target fixed while rotating the H+ projectile around. The definite initial conditions of the H+ projectile RH+i and PH+i are therefore obtained by rotating RH+0 and PH+0 through the extrinsic Euler angles [67] in the order: 1st, 00γ < 3600, 2nd, 00β ≤ 1800, and 3rd, 00α < 3600, 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 RH+0 and PH+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, RH+i and PH+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 02π0π02πDMMJ(α,β,γ)sinβdαdβdγ=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+ + (H2O)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+ + (H2O)1–6 → H + different cluster products[69]

σ1ET=14π002π0π02πbP1-ET(α,β,γ,b)dbsinβdαdβdγ
(7)

where P1−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, P1−ET (α, β, γ, b) → P1−ET (b), and Eq (7) simplifies to the familiar expression: σ1ET=2π0bP1ET(b)db [42]. In the present systems, an outgoing projectile can in practice capture up to two electrons: H+ + A → H1−n + An, 0 ≤ n ≤ 2, because the probability of forming unstable H1−n with n ≥ 3 is negligible. Under these conditions, Pn−ET (α, β, γ, b), 0 ≤ n ≤ 2, are [69]

P0ET(α,β,γ,b)=(1Nα)(1Nβ);P1ET(α,β,γ,b)=Nα(1Nβ)+Nβ(1Nα);P2ET(α,β,γ,b)=NαNβ
(8)

where Nα and Nβ are the number of α- and β-spin electrons in the outgoing projectile calculated from their respective electron densities ραout.proj.(r) and ρβout.proj.(r):

Nα=ραout.proj.(r)dr;Nβ=ρβout.proj.(r)dr
(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 P¯1ET(b), P¯1ET(b)=(8π2)102π0π02πP1ET(α,β,γ,b)sinβ, and their b – weighted counterparts, bP¯1ET(b), where σ1-ET=2π0bP¯1ET(b)db. P¯1ET(b) and bP¯1ET(b) reveal more mechanistic details than σ1-ET.

Results and discussion

Electron-transfer properties

The first property calculated in this investigation is the cluster-to-proton total 1-ET ICS, σ1−ET, for the H+ + (H2O)n systems, n = 1–6, at ELab = 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 [7073], 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, σ¯1ETExp., and its average relative error, eExp., are 1.27 Å2 and ± 10.62%, respectively. The theoretical σ1ETTheo. ‘s and their average relative deviations Δ¯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 (Δ¯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 (Δ¯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.

Table 2
Water-to-proton 1-electron-transfer total integral cross sections (ICSs) σ1−ET for H+ + H2O at ELab = 100 keV from experiments, SLEND theory and alternative theories: basis generator method (BGM) and continuum distorted wave-eikonal initial ...

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+ + (H2O)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) = cn2/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 (H2O)n clusters should be approximately proportional to their size n, V(n) [proportional, variant] 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)πRn3 and A(n)=4πRn2, respectively, where Rn is the radius of the (H2O)n sphere; therefore, A(n) [proportional, variant] V (n)2/3 [proportional, variant] n2/3. In turn, the σ1−ET are proportional to the effective external area A(n) of the (H2O)n exposed to the incident H+ for ET processes [42]; therefore, σ1−ET [proportional, variant] A(n) [proportional, variant] n2/3 [implies] σ1−ET (n) = cn2/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) = cn2/3 with correlation factors R2 = 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., σ1ET631G* = 5.96 and 6.51 Å for n = 7 and 8, respectively; these estimated values should be interpreted as average σ1−ET over the various (H2O)7 and (H2O)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 H2O 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 ELab = 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) = cn2/3 formulae differing in their coefficient c would exclusively connect values from different sets of clusters—e.g. a single σ1−ET (n) = cn2/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 (H2O)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 (H2O)1-6 isomers. Thus, unlike the case of structural and solubility properties [3537], the “magic number” of six waters in the prism and cage (H2O)6 isomers do not bring about any hint of water radiolysis processes in bulk water. Likely, an extension of the current (H2O)1-6 series with the (H2O)7-20 isomers may bring about some type of bulk-water manifestations.

Fig 3
SLEND/6-31G* and /6-31G** target-to-proton total 1-ET ICSs σ1−ET for H+ + (H2O)1-6 vs. the water cluster size n.
Table 3
SLEND cluster-to-proton 1-electron-transfer integral cross sections σ1−ET for H+ + (H2O)n, n = 2–6, at ELab = 100 keV.

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 P¯1ET(b) and their b – weighted counterparts bP¯1ET(b) vs. b for the considered (H2O)1-6 isomers. With both basis sets, the P¯1ET(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¯1ET(b) vs. b plots show a variable number of maximum peaks (from one to three) depending on the considered water cluster. The bP¯1ET(b) vs. b plots exhibit similar patterns to those of the P¯1ET(b) but modulated by the b value. As the integrands of the σ1−ET, the bP¯1ET(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.).

Fig 4
Orientation-averaged target-to-proton 1-ET probabilities P¯1ET(b) at the SLEND/6-31G* (left panel) and /6-31G** (right panel) levels vs. the impact parameter b for the investigated (H2O)1-6.
Fig 5
Orientation-averaged impact-parameter-weighted target-to-proton 1-ET probabilities bP¯1ET(b) at the SLEND/6-31G* (left panel) and /6-31G** (right panel) levels vs. b for the investigated (H2O)1-6.

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 the 1-ET processes. To study those reactions, the SLEND/6-31G* simulations to calculate σ1−ET of H+ + (H2O)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)[right angle bracket] is a superposition of various UHF states corresponding to various products’ channels [69, 75], e.g. Hproj.++H2OHproj.q1=+1+H2q2=0+Oq3=0 or Hproj.q1=0+H2q2=+1+Oq3=0 or Hproj.q1=1+H2q2=+1+Oq3=+1, etc., where Hproj.q1=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 qi (e.g. q1 = 1−NαNβ for the final Hproj.q1) are not necessarily integer numbers corresponding to canonical chemical species (e.g. Hprojq1=+1, Hprojq1=0 and Hprojq1=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. Hproj.+1+(H2O)nHproj.q1+(H2O)n1q1, with qi continuously varying in the range −1≤qi≤+1). However, as in previous SLEND studies of H+ + (H2O)1–4 at ELab = 1 keV [6], the present simulations leading to clusters fragmentations always bring about outgoing projectiles Hproj.q1 and clusters fragments Aq2 with q1 = +1 and q2 = 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 [7678]. 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.

Fig 6
Example of a SLEND/6-31G* simulation of H+ + H2O leading H2O to fragment into H2 and O.

The predicted fragmentations are:

H+ + H2O simulations: 9 out of 252 simulations exhibited the H2O target fragmenting into: H + OH (2 simulations), H + H + O (6 simulations) and H2 + O (1 simulation, cf. Fig 6).

H+ + (H2O)2 simulations: 9 out of 540 simulations exhibited the (H2O)2 target fragmenting into: H2O + HO + H (2 simulations), H3O + O + H (1 simulation), H2O + 2H + O (5 simulations) and H3O + OH (1 simulation, cf. Fig 7).

Fig 7
Target fragments from individual SLEND/6-31G* simulations of H+ + (H2O)2..

H+ + (H2O)3 simulations: 3 out of 540 simulations exhibited the (H2O)3 target fragmenting into: H3O + OH + H2O (2 simulations) and H + OH + 2H2O (1 simulation).

H+ + (H2O)4 (tetramer a) simulations: 10 out of 748 simulations exhibited the (H2O)4 target fragmenting into: (H2O)2 + 2H2O (4 simulations), H3O + OH + 2H2O (4 simulations), an H + OH + 3H2O (2 simulations)

In conclusion, present SLEND/6-31G* simulations predict the DNA-damaging radicals H, OH, O and H3O, and the innocuous species H2O and (H2O)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 Hproj.++H2OHproj.++H2+O and a sixth panel plotting the Mulliken populations of the Hproj., H2 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 Hproj.+ passes in between the H atoms of H2O, hits the O atom and bounces back. As a result of this collision, the H2 and O moieties of H2O break apart; the ejected H2 moiety undergoes a series of strong oscillations ranging from the near dissociation of H2 into H atoms to these atoms’ recombination back to H2. Finally, Fig 7 shows the nine predicted fragments in H+ + (H2O)2.

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 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 (H2O)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 P1−ET (cf. Eqs 7 and 8). However, as derived in Ref. [69] and supposed in previous SLEND studies [29, 75], the Pn−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[center dot] + A[center dot]+ [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[center dot]+ [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 Pn−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 Pn−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 Pn−ET in Eq 8 and resulting σnET 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 σnET 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 ≤ ELab ≤ 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 P1−ET in Eq 8 for σ1−ET, the final-time electronic wavefunction |ΨeSLEND(tf), Eq 6, was projected onto the ground |0[right angle bracket] and excited states |… hp[right angle bracket] 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 ≤ ELab ≤ 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 ELab = 60 keV: σDISLEND = 13.9Å and σDIEXPT. = 13.8 Å [implies] 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+ + (H2O)n simulations and ruled it out as a source of the SLEND/6-31G**σ1−ET overestimation. For H+ + H at ELab = 100 keV, SLEND σDI is 0.92 Å2. If P1−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+ + H2O, the overestimated SLEND/6-31G** σCTSLEND via Eq 8 can be corrected by subtracting the σDISLEND part from it: σCTSLEND = 2.06 Å2–0.92 Å2 = 1.14 Å2; this places SLEND/6-31G** σCTSLEND within the range of the four experimental values and closest to that from Exp. D [73]: σCTExpt.D = 1.13 Å [implies] ΔTheo. = 0.8%. A similar correction might occur with the SLEND/6-31G* σCTSLEND 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+ + H2O is more complicated than that of H+ + H because the former has more than one electron. For H+ + H2O, numerous excited states |… hp[right angle bracket] from |0[right angle bracket] forming a full CI expansion should be generated so that |ΨeSLEND(tf) in Eq 6 can be projected on all those states. This more demanding capability is not currently available in PACE but is under development.

Conclusions

To model microscopic processes in PCT [17], the SLEND method was applied to the H+ + (H2O)n systems at ELab = 100 keV with the 6-31G* (n = 1–6) and 6-31G** (n = 1–5) basis sets. Ten (H2O)1–6 clusters were selected for this study: eight exhibit mono- and di-cyclic quasi-planar/multiplanar structures [3335] and two others, the prism and cage (H2O)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 [7073] 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) = cn2/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 H3O. 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.

Acknowledgments

The authors thank Prof. Takao Oi (Sophia University, Japan) for providing the optimized geometry of the prism isomer of (H2O)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.

Funding Statement

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

Data Availability

All relevant data are within the paper.

References

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 ELab = 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