Chieh-Min
Hsieh
a,
Alexander
Dellwisch
a and
Tim
Neudecker
*abc
aUniversity of Bremen, Institute for Physical and Theoretical Chemistry, Leobener Str. 6, D-28359 Bremen, Germany. E-mail: neudecker@uni-bremen.de
bBremen Center for Computational Materials Science, Am Fallturm 1, D-28359 Bremen, Germany
cMAPEX Center for Materials and Processes, Bibliothekstr. 1, D-28359 Bremen, Germany
First published on 11th June 2025
Estimating the projected range of high-energy particles is important for ion implantation and designing shielding strategies for space devices. In this work, we propose a molecular dynamics (MD) workflow to calculate the projected range of protons and demonstrate its capabilities by calculating the projected range of protons in graphite and poly(methyl methacrylate) (PMMA). The results show excellent agreement with reference data. Besides, we investigate irradiation-induced bond breaking by simulating the proton bombardment of a perylene-3,4,9,10-tetracarboxylic dianhydride (PTCDA) molecule and analyze the strain energy accumulated in the system using quantum chemical tools. The findings indicate a correlation between strain energy and the kinetic energy of the primary knock-on atom.
One common approach to estimate the projected range is the binary collision approximation (BCA),9 which models the recoil event as a series of collisions between two particles. This method is effective and computationally efficient for simulating particle collisions. The widely used Stopping and Range of Ions in Matter (SRIM) code,10 based on the BCA and Monte Carlo11 approach, is a prominent example in this category. However, a limitation of this method is that it models the trajectories of particles only in amorphous materials. Later, molecular dynamics (MD) for calculating the projected range gained attention. A notable first attempt was the MDRANGE code,12 developed by Nordlund in the 1990s. Although MD is computationally more expensive than BCA, it models the crystal structures atomistically, enabling more accurate range prediction for crystalline materials. As a result, MD has successfully estimated the range of several crystalline materials, showing its potential in this field.13,14
Although the aforementioned methods have been effective, they were developed decades ago and have significant limitations, as they primarily target inorganic solids. For instance, polymers—important materials for space applications—may not be accurately analyzed using these methods due to their more complex chemical bonding patterns. In this paper, we revisit calculating the projected range of protons in matter using the state-of-the-art Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) MD code.15 Specifically, we employed MD to calculate the projected range of high-velocity protons in two chemically very diverse materials: graphite and poly(methyl methacrylate) (PMMA). The former has a stacked layer structure, while the latter is a polymer. For these simulations, we used the reactive force field (ReaxFF) potential,16 which includes elements commonly found in organic materials (e.g., H, C, O, and N) and is considered a well-established method for modeling bond rupture events in such systems. For both graphite and PMMA, which contain typical organic bonds (CO, C–H, C–C, and aromatic structures), our use of ReaxFF yielded equilibrium densities in reasonable agreement with experimental values, supporting the validity of our chosen potential. Last but not least, we explored the degradation scenarios, such as bond breaking, and applied the quantum chemical Judgement of Energy DIstribution (JEDI) strain analysis17,18 to investigate the strain energy induced by the impacting proton in a single perylene-3,4,9,10-tetracarboxylic dianhydride (PTCDA) molecule.
![]() | (1) |
When a projectile interacts with atoms, the force acting between them is reduced by a friction force, as
![]() | (2) |
The other key parameter is the Ziegler–Biersack–Littmark (ZBL)19 repulsive interatomic potential, which is defined as
![]() | (3) |
![]() | (4) |
Finally, adaptive time steps are essential to make the time step inversely proportional to the projectile's velocity.20 This means that a shorter time step is applied when the projectile's velocity is high, while a longer time step is used when the projectile slows down. This approach ensures both the accuracy and efficiency of the calculation.
In our approach, the proton is initially placed at a random position at the boundary of the simulation cell with a specified initial velocity in the z-direction, as demonstrated in Fig. 1. An MD simulation is then initiated. During the simulation, the system evolves without a thermostat or barostat, and the temperature rises during the simulation, reflecting an adiabatic process where the proton deposits energy into the material over a short timescale. In our previous tests, applying NVT or NPT conditions in this step led to a significant underestimation of the proton's projected range. Simultaneously, the proton's kinetic energy is monitored until the simulation terminates under one of two conditions: (1) the proton reaches the boundary of the simulation cell again. In this situation, the proton's velocity is recorded and used to initialize a new MD simulation. (2) the proton's kinetic energy drops below 1 eV (indicating the proton is stopped). In this case, the simulation is done, and the projected range is determined as the sum of the penetration distances in the z-direction across all individual MD runs. The workflow of the simulation process is illustrated in Fig. 2.
The MD simulations of the projected range were carried out using the LAMMPS code. Initially, the system was equilibrated using the NPT ensemble with a Nosé–Hoover thermostat22,23 and a Nosé–Hoover barostat with Martyna–Tuckerman–Klein (MTK) corrections24 at 50 K and 1 atm. A proton was then positioned at the boundary of the simulation cell on the xy-plane with its z-coordinate set to 0, while its x- and y-coordinate were assigned randomly. The proton was given an initial velocity corresponding to a kinetic energy of 1, 5, 10, 20, or 50 keV, directed along the negative z-direction. The simulation then proceeded by following the MD workflow (Fig. 2) to determine the projected range. For each initial proton kinetic energy value, 20 individual trajectories were calculated to obtain a statistically meaningful distribution.
Periodic boundary conditions were applied in all calculations. The ReaxFF potential (ffield.reax.FC)25,26 distributed with the LAMMPS code was used to model the systems, with the ReaxFF algorithm implemented in LAMMPS, as described in the corresponding reference.27 The ZBL19 potential with inner and outer cutoffs of 3.0 and 4.0 Å, respectively, was applied to account for the repulsive interactions between the projectile and the atoms within the target material. In addition, stopping power from the PSTAR database was employed to account for the inelastic energy loss of the proton during penetration. It should be noted that stopping powers below 1 keV are not available in the PSTAR database. Therefore, the stopping power provided by SRIM was used for this energy range. Finally, the time step was adjusted using the dt/reset command to limit the proton's maximum travel distance per simulation step to 0.005 Å, ensuring the time step was appropriately assigned.
To validate our simulation, we utilized the PSTAR database and the SRIM code as references. The PSTAR database is provided by the National Institute of Standards and Technology (NIST). It contains a stopping power table of protons, which is derived from a combination of experimental data and theoretical predictions based on the Bethe formula.28 This database also provides the projected range R, which is calculated by integrating the reciprocal stopping power over the energy range, as shown in
![]() | (5) |
SRIM calculations were performed using the 2013 version of the SRIM code.10 The ion was hydrogen with a mass of 1.008 amu. The graphite target consisted of carbon atoms with a density of 1.70 g cm−3, while the PMMA target was composed of hydrogen, carbon, and oxygen atoms in an atomic ratio of 8:
5
:
2, with a density of 1.18 g cm−3.
To increase the density, an NPT run using Nosé–Hoover thermostat and Nosé–Hoover barostat with MTK corrections was performed at 300 K and 1 atm for 0.5 ns. This was followed by a 21-step relaxation procedure34 for enhanced equilibration, involving various compression and decompression steps within a pressure range of 1 to 50000 atm and a temperature range of 300 to 1000 K over a total simulation time of 1.56 ns. This 21-step equilibration process has been used for the preparation of various systems in the past.35–37 During equilibration, the simulation cell was resized, resulting in a density of 1.11 g cm−3, which is close to the experimental value of 1.17 g cm−3.38 The equilibrated system obtained from the 21-step equilibration procedure was then used as the initial configuration for calculating the projected range.
![]() | ||
Fig. 3 A PTCDA molecule with the highlighted rectangular area indicating the region impacted by the proton. |
1 keV | 5 keV | 10 keV | 20 keV | 50 keV | |
---|---|---|---|---|---|
Graphite | |||||
Avg. MD range | 182 | 793 | 1314 | 2472 | 4994 |
SD | 57 | 175 | 297 | 424 | 409 |
SRIM range | 220 | 933 | 1630 | 2745 | 5385 |
PSTAR range | 179 | 837 | 1502 | 2602 | 5237 |
Dev. from SRIM range | −17% | −15% | −19% | −10% | −7% |
Dev. from PSTAR range | 2% | −5% | −13% | −5% | −5% |
PMMA | |||||
Avg. range | 225 | 876 | 1616 | 2958 | 5727 |
SD | 70 | 264 | 263 | 220 | 415 |
SRIM range | 281 | 1188 | 2065 | 3450 | 6668 |
PSTAR range | 224 | 1026 | 1809 | 3074 | 6092 |
Dev. from SRIM range | −20% | −26% | −22% | −14% | −14% |
Dev. from PSTAR range | 0% | −15% | −11% | −4% | −6% |
![]() | ||
Fig. 4 Comparison of the projected range of protons in graphite and PMMA calculated in this study with reference data from PSTAR and SRIM. Error bars represent one standard deviation. |
To investigate this discrepancy, we compared the stopping power of the two materials, defined as the proton's energy loss per unit penetration distance (eqn (1)), as provided by SRIM and PSTAR. As illustrated in Fig. 5, the stopping power for PMMA obtained from PSTAR is slightly higher than that from SRIM. This difference is expected to result in a smaller calculated projected range when the stopping power from PSTAR is used. It is important to note that our calculations combined the stopping power data from both PSTAR and SRIM. Specifically, for the kinetic energy values from 1 keV to 50 keV, the stopping power was taken from PSTAR, while for kinetic energy values below 1 keV, the stopping power was obtained from SRIM, as PSTAR does not provide data for this range. This leads to a larger calculated range, especially for calculations with a low-energy proton. As seen in the trend of the blue curve in Fig. 4—when compared to the red and green curves—we observed that the data points for graphite at 1 and 5 keV, as well as PMMA at 1 keV, are higher than expected. This anomaly can be attributed to the lower stopping power from SRIM, which leads to a larger calculated projected range.
![]() | ||
Fig. 5 Comparison of the stopping power of protons in graphite and PMMA over the kinetic energy range of 1 keV to 100 keV. |
It is worth mentioning that, during the simulation, the proton can change its direction, indicating the proton collided with an atom in the material and rebounded in a different direction. This behavior can be attributed to the repulsive force generated by the incorporated ZBL potential, which is crucial for decelerating the proton. In our preliminary calculations, we observed that without the ZBL potential, the projected range of the proton was highly overestimated.
Moreover, we conducted five additional simulations for graphite with proton energies of 1 keV and 5 keV, respectively. Our results indicate that the energy loss due to electronic stopping is trajectory-dependent. For 1 keV protons, the fraction of kinetic energy lost to electronic stopping ranges from 54% to 70%, whereas for 5 keV protons, this increases to 83% to 92%. This trend is expected, as higher-energy protons travel farther and interact with more electrons, resulting in greater cumulative energy loss via electronic stopping. In contrast, lower-energy protons undergo more frequent nuclear collisions, which contribute more significantly to nuclear stopping—that is, energy loss due to repulsive interactions described by the ZBL potential.
In conclusion, the overall satisfactory results lend credibility to our simulation approach, and this opens up the possibility to simulate the projected range of protons in a wide range of organic materials in the future.
Shi et al.43 applied classical MD to simulate defect generation in graphene bombarded by a proton. Specifically, they employed the Tersoff44 potential combined with the ZBL potential. They compared their MD results with ab initio MD (AIMD) simulations and observed consistency. In our study, we adopted a similar approach but used the ReaxFF potential instead of the Tersoff potential due to its greater compatibility with diverse atom types. Since both Tersoff and ReaxFF are bond-order potentials, they are well-suited for modeling bond breaking in MD simulations.
We conducted an MD simulation of a proton impacting a PTCDA molecule, which, similar to graphene, has a two-dimensional layered structure with a π-conjugated system. Besides, it contains hydrogen and oxygen atoms, which are common in organic compounds. The initial kinetic energy of the proton was set to 1 and 10 keV. Given the D2h symmetry of PTCDA, only one-fourth of the impact area needed to be analyzed. To thoroughly evaluate the whole impact region, a two-dimensional scan with a resolution of 0.1 Å in the xy-plane was performed, as illustrated in Fig. 3.
The results of bond breaking analysis suggest that the impact position must be close to the primary knock-on atom (PKA) to break bonds, as illustrated in Fig. 6, where the red circles (bond breaking occurs) align with the nuclear configuration of a PTCDA molecule. This observation is consistent with the findings of Shi et al.43 Besides, the number of red circles demonstrates the stability of the corresponding atom. For instance, the carbon atoms in the aromatic ring of the perylene core are more stable than the carbon atom in the carboxylic anhydride group. Moreover, as the proton's kinetic energy increases, the number of red circles decreases, implying a reduced probability of bond breaking. This trend arises because higher-energy protons have a shorter interaction time to transfer energy to the PKA.43
We also analyzed the maximum kinetic energy of the PKA during our simulations for different impact positions. As shown in Fig. 7, a comparison of the kinetic energy of the PKA for bond-breaking and non-bond-breaking cases reveals a significant energy difference, confirming that the energy transferred to the PKA is decisive in damaging chemical bonds.
![]() | ||
Fig. 7 MD simulation results showing the maximum kinetic energy of the PKA in the highlighted area of Fig. 6, expressed in kcal mol−1. |
Fig. 8 presents the MD snapshots of the PTCDA molecule just before the C–C bond breaking, with a color scale highlighting the strain energy distribution. The analysis reveals that the strain energy is highly localized near the impact position. The strain energy stored in the C–C bond before bond breaking amounts to 96 kcal mol−1. Considering the resonance structure of the aromatic ring, the C–C bond strength should lie between a C–C single bond (ethane, 90 kcal mol−1) and a CC double bond (ethylene, 174 kcal mol−1).45 Although the strain energy calculated by JEDI falls within this range, the value appears to be slightly lower than expected. In cases where bond rupture did not occur, the strain energy was redistributed throughout the entire molecule, resulting in molecular vibrations. Videos demonstrating the dynamic strain energy distribution are available in the ESI.†
Furthermore, we compared the maximum strain energy of the distorted PTCDA molecule and the maximum kinetic energy of the PKA during the MD simulation, as demonstrated in Table 2. The maximum strain energy was determined for each impact position by analyzing the geometry of the PTCDA molecule in each MD trajectory using JEDI. The strain energies for four of the MD trajectories could not be calculated because the JEDI strain analysis is based on the harmonic approximation, which becomes invalid when bond breaking occurs. For the remaining cases, the kinetic energy value correlates with the strain energy, indicating that the kinetic energy of the PKA is converted into strain energy within the distorted structure. However, some relatively higher deviations are observed, which could be attributed to differences in bond length estimated by the ReaxFF potential compared to DFT calculations underlying the JEDI analyses. In addition, the following points need to be kept in mind: (1) the harmonic approximation used in JEDI generally leads to an underestimation of strain energy in compressive cases and an overestimation in stretching cases.18 (2) In dynamic JEDI, only potential energy is quantified; the contribution of kinetic energy of the strained system is neglected, and the relative contributions of each component remain unknown.17 (3) Strain induced by thermal vibrations may also contribute to the total strain of the molecular system. Combined, these effects explain the deviations between the kinetic energy of the proton and the strain energy of the system calculated with JEDI.
6–14 | 7–14 | 8–14 | 9–14 | |
---|---|---|---|---|
max.ke.PKA | 108 | 265 | 246 | 94 |
max.strain | 105 | 249 | 239 | 112 |
Dev. | −3% | −6% | −3% | 19% |
6–13 | 7–13 | 8–13 | 9–13 | |
---|---|---|---|---|
max.ke.PKA | 236 | 1834 | 1394 | 191 |
max.strain | 181 | — | — | 164 |
Dev. | −23% | — | — | −14% |
6–12 | 7–12 | 8–12 | 9–12 | |
---|---|---|---|---|
max.ke.PKA | 198 | 974 | 816 | 164 |
max.strain | 171 | — | — | 151 |
Dev. | −14% | — | — | −8% |
6–11 | 7–11 | 8–11 | 9–11 | |
---|---|---|---|---|
max.ke.PKA | 79 | 159 | 150 | 70 |
max.strain | 103 | 145 | 146 | 106 |
Dev. | 30% | −9% | −3% | 51% |
Our JEDI strain analysis indicates that the strain induced by the proton is highly localized, and concentrated near the impact position. The strain energy is then redistributed throughout the molecule if the bond rupture does not occur, leading to increased atomic oscillations and, ultimately, heating of the system. Additionally, by comparing the kinetic energy of PKAs and the strain of the distorted PTCDA molecule, we conclude that most of the PKA's kinetic energy is converted into strain energy.
In our simulations, the proton is treated as a neutral particle, which is a simplified approximation. While this assumption has minimal influence on the calculated projected range, it may affect the material in other ways. For example, it could impact proton diffusion and the generation of radicals within the structure.46,47 To address this limitation, the accuracy provided by quantum mechanics (QM) is desirable. We see the potential of the electron force field (EFF) method as a solution to this issue.48,49 This method approximates quantum effects through the quantum wave packet, thereby providing a more accurate description of the electronic structure. Another promising approach is using machine learning force field (MLFF) potentials,50 which could enhance accuracy while maintaining the computational efficiency of classical MD simulations.
Footnote |
† Electronic supplementary information (ESI) available: Scripts for calculating the projected range of protons using LAMMPS, scripts for modeling proton bombardment on a PTCDA molecule, and videos showing dynamic strain distribution in the PTCDA molecule under proton bombardment are available. See DOI: https://doi.org/10.1039/d5cp00223k |
This journal is © the Owner Societies 2025 |