Advances in multiscale modeling for novel and emerging technologies

Abstract

Computational multiscale modeling encompasses a wide range of end-products and a great number of technological applications. This paper provides an overview of the computational multiscale modeling approach based on utilization of MBN Explorer and MBN Studio software packages, the universal and powerful tools for computational modeling in different areas of challenging research arising in connection with the development of novel and emerging technologies. Three illustrative case studies of multiscale modeling are reviewed in relation to: (i) the development of novel sources of monochromatic high-energy radiation based on the crystalline undulators, (ii) controlled fabrication of nanostructures using the focused electron-beam induced deposition, and (iii) ion-beam cancer therapy. These examples illustrate the key algorithms and unique methodologies implemented in the software.

Graphic abstract

draftps

1 Introduction

Computational multiscale modeling encompasses a wide range of end-products and a significant number of applications in (i) the avionics and automobile industry for designing nanostructured materials, functionalized surface coating, as well as stronger and lighter materials for aircraft and cars; (ii) mechanical engineering for virtual design of superhard nanostructured materials; (iii) medical applications for novel materials technologies for implants and tissue regeneration; (iv) electronic and chemical industry for constructing highly efficient batteries and catalyzers; (v) the pharmaceutical industry for drug design, etc. In most of these applications, it is necessary to identify and design specific properties of the system determined by its molecular structure on the nanoscale and ensure the transfer of these properties to the macroscopic scale to make them functional and usable. Such a transition implies multiscale modeling approaches that rely on the combined use of quantum mechanics methods together with classical molecular dynamics (MD), or linkage of MD and Monte Carlo (MC) simulations, or the application of efficient computational algorithms allowing to perform simulations across the scales.

This paper provides an overview of the computational multiscale modeling approach [1] based on utilization of the powerful and universal software packages MBN Explorer [2] and MBN Studio [3]. The paper begins with a brief overview of different theoretical methods for modeling Meso-Bio-Nano (MBN) systems and the limitations of these methods. It is followed by a brief description of the MBN Explorer and MBN Studio software packages and their application areas. Next, several case studies of multiscale modeling by means of MBN Explorer and MBN Studio are reviewed, illustrating the main algorithms and unique methodologies implemented in the software.

There are many concrete examples of novel and emerging technologies benefiting from computational multiscale modeling. As an illustration, a few of them are highlighted in the paper: (i) development of novel sources of monochromatic high-energy radiation based on the crystalline undulators, (ii) controlled fabrication of nanostructures using the focused electron-beam-induced deposition (FEBID), and (iii) ion-beam cancer therapy. These examples illustrate the new possibilities that computational multiscale modeling provides novel and newly emerging technologies.

2 Computational approaches in meso-bio-nano-science

MBN Science is an interdisciplinary field of research studying structure formation and dynamics of animate and inanimate matter on the nano- and mesoscales. This field bundles up several traditional topics in theoretical physics and chemistry at the interface with life sciences and materials research under a common theme.

The range of open challenging scientific problems in the field of MBN Science is very broad [1, 4]. They may include: structure and dynamics of clusters, nanoparticles, biomolecules, and many other nanoscopic and mesoscopic systems; clustering, self-organization, growth and structure-formation processes, and their multiscale nature; assemblies of clusters/nanoparticles and bio-macromolecules; hybrid bio-nanosystems; nanostructured materials; surface phenomena; nanoscale phase and morphological transitions; thermal, optical and magnetic properties; collective and many-body phenomena; electron transport and molecular electronics; collision, fusion, fission and fragmentation processes; channeling effects; radiation effects; radiobiological effects, and many more.

MBN Science has been boosted immensely over the past two decades by the fast development of computer powers and related computational techniques that became broadly available. This development resulted in a significant increase in the efficiency of available computer codes for scientific research. Such codes are usually focused on particular systems, their particular sizes and phenomena involved (see Fig. 1) and thus are limited in their ability to model physical, chemical or biological phenomena that go across the scales. Therefore, much effort of research communities has been devoted to the computational approaches and modeling techniques to overcome this drawback and open new horizons in theoretical and computational research.

Fig. 1

An illustration of entities in the different computational techniques

2.1 Atomic and nanoscales

The characteristic quantum processes in MBN systems often involve electronic excitation, relaxation, fragmentation, or transport dynamics. Often, these effects in molecular and nanosystems occur when they are embedded into larger-scale molecular environments. The quantitative description of the mentioned processes requires inclusion of both short-time/(sub-)nanoscale quantum aspects and long-time/nano- and mesoscale environment effects. Dynamical descriptions of molecular systems on the quantum atomic/subnano- and nano-/mesoscales are presently performed with disconnected theoretical and simulation tools. There is no serious hope to explore in the near future both ranges by simply extending one tool to the other domain. Bridging this gap implies interfacing quantum atomic/molecular and nano-/mesoscopic approaches. These efforts will result from merging into a common model a time-dependent quantum mechanical (QM) approach and classical molecular dynamics (MD). This will enable exploring new types of dynamics in a broad range of molecular systems, not accessible for analysis by other theoretical and computational means.

Attempts have been made during recent years in various realizations of quantum mechanics/molecular mechanics (QM/MM) methods. State-of-the-art achievements of these methods have been reviewed in several publications [5,6,7], and the field will experience a significant development in the future. Many well-established quantum and classical MD codes, being developed for a long time entirely independently, now reach the point when they require coupling to go further across disciplines. As this goal could not be achieved in a single universal approach, its realization often requires the solution for each particular pair of quantum and classical codes with their particular application areas and the research communities standing behind.

2.2 Nano- and meso-scales

Description of structure and dynamics of molecular and nanosystems based on classical MD principles provides significant computational advantages compared to the quantum descriptions due to the simplicity of the formalism involved. Numerical realizations of this approach nowadays enable to treat the structures consisting of tens of millions of atoms and access dynamical timescales of molecular processes up to tens of microseconds, although for smaller systems, typically on the scale of hundreds of thousands of atoms. Achievements made on the construction of software packages built on these principles have been reviewed in [8,9,10]. The existing codes are based on the utilization of different classical force fields suitable for describing certain types of molecular systems, e.g., biomolecular systems, carbon or metallic nanosystems, etc. However, many of the MD packages do not have the possibility (at least in the sufficiently elaborated form) to combine different force fields for the description of hybrid molecular systems and molecular environments due to the high level of optimization tuned for the specific type of force fields.

In recent years, it has become clear that a detailed understanding of numerous quantum molecular processes happening in larger-scale molecular environments becomes possible due to new advances in theoretical and experimental tools developed in atomic and molecular physics and nanotechnology. An impressive example is the discovery of resonant mechanisms at the origin of radiation damage caused by low-energy electrons [11, 12] via the dissociative electron attachment (DEA) mechanism. The identification of such a quantum effect in this context had an enormous impact: it paved the way to new research on possible low-energy radiation damage effects [13]. Another example where the detailed knowledge of atomic and molecular processes in larger-scale molecular environments is crucial concerns novel radiosensitizers, e.g., metallic nanoparticles, which can enhance radiation effects on biological targets [14]. An important feature of the aforementioned processes is that they occur on relatively short times (tens of femtoseconds) and in relatively small spatial domains (up to a few nanometers).

By exploiting the locality of the QM processes, one can isolate the active part of the system and describe it at the pure QM level. However, in this case, possible large-scale effects coming from the other parts of the system might be missing. For instance, this can be the case for a chromophore molecule in solvating shells or impurities in a solid. The size of the QM system can be enlarged, but the convergence of such a strategy, in terms of accuracy and/or simulation time, can become slow. Instead, one can explicitly include the environment and treat it classically at the MM level. The development of such hierarchical techniques over past decades mainly operates in the field of biological chemistry or surface science, where large sizes and/or environmental effects are crucial. These methods allow one to tackle systems of thousands of atoms (a few tens at the QM level).

All-atom MD simulations involving very large systems may require significant computer resources that are not easily affordable. Similarly, simulations of processes on long timescales (beyond $\sim$1 ms) are prohibitively expensive because they require carrying out too many integration time steps. In these cases, one can sometimes tackle the problem by using reduced representations of the system where, instead of explicitly representing every atom of the system, groups of atoms are represented by “pseudo-atoms”. Such representations are sometimes called the coarse-grained models [15]. The parameterization of these models must be done empirically by matching the model’s behavior to appropriate experimental data or all-atom simulations. Coarse-grained models have been used successfully to examine a wide range of questions in structural biology, liquid crystal organization, and polymer glasses [15].

2.3 Monte Carlo approach and finite element method

Kinetic Monte Carlo (KMC) method is designed to model the time evolution of many-particle systems stepwise in time. Instead of solving dynamical equations of motion, the KMC approach assumes that the system undergoes a structural transformation at each step of evolution with a certain probability. The new configuration of the system is then used as the starting point for the next evolution step. The transformation of the system is governed by several kinetic rates chosen according to the model considered. Due to its probabilistic nature, this methodology permits studying dynamical processes involving complex molecular systems on timescales significantly exceeding the characteristic timescales of conventional MD simulations [16,17,18]. The KMC method is ideal in situations when certain minor details of dynamic processes become inessential, and the major transition of the system to new states can be described by a few kinetic rates, determined through the corresponding physical parameters.

Quantitative structural analysis of macroscopically large systems and the description of some processes occurring on the macroscopic scale (such as heat transfer, fluid flow, mass transport, etc.) can be effectively carried out by means of finite element methods (FEM) [19]. It is generally impossible to obtain analytical mathematical solutions for systems involving complicated geometries, loadings, and material properties. FEM enables formulating the problem in terms of a set of algebraic equations rather than requiring the solution of differential equations. These numerical methods yield approximate values of the unknowns at discrete numbers of points in the continuum. Hence, this process of modeling a body by dividing it into an equivalent system of smaller bodies or units (finite elements) interconnected at points common to two or more elements (nodal points or nodes) and/or boundary lines and/or surfaces is called discretization. In FEM, instead of solving the problem for the entire body in one operation, one formulates the equations for each finite element and combines them to obtain the solution of the whole body.

3 MBN Explorer and MBN Studio

Multiscale modeling of MBN systems is one of the hot topics of modern theoretical and computational research. To fully understand and exploit all the richness and complexity of the MBN world, especially its dynamics, one needs to consult many disciplines ranging from physics and chemistry to materials and life sciences, exploiting technologies from software engineering and high-performance computing. This general trend brought up the idea and then the development of MBN Explorer [2] and MBN Studio [3]. These software packages have been designed as powerful and universal instruments of computational research in the field of MBN Science, which are capable to explore, simulate, record, and visualize both structure and dynamics of MBN systems, reproduce its known features, and predict the new ones.

MesoBioNano Explorer (MBN Explorer) [2] is a software package for the advanced multiscale simulations of the structure and dynamics of various MBN systems. It has many unique features and a wide range of applications in physics, chemistry, biology, and materials science. It is suitable for classical, non-relativistic and relativistic molecular dynamics (MD), Euler dynamics, reactive and irradiation-driven molecular dynamics (RMD and IDMD) simulations, as well as for stochastic dynamics or kinetic Monte Carlo (KMC) simulations of various randomly moving MBN systems or processes. These algorithms are applicable to a broad range of systems, such as nano- and biological systems, nanostructured materials, composite/hybrid materials, gases, plasmas, liquids, solids, and their interfaces, with sizes ranging from atomic to mesoscopic.

3.1 Main features of MBN Explorer

MBN Explorer enables calculations of energies of a large variety of MBN systems and optimization of their structures. The software package supports different types of molecular dynamics for MBN systems. The program operates with an extensive library of interatomic potentials, thus allowing to model many different molecular systems. It is also possible to simulate the dynamics of MBN systems in the presence of external fields—electric, magnetic, gravitational, and electromagnetic waves.

Apart from the standard algorithms, MBN Explorer is equipped with unique algorithmic implementations that significantly enhance the computational modeling capacities in various research and technological areas. The complete list of algorithms implemented in MBN Explorer can be found in [1, 20, 21].

In particular, MBN Explorer has unique algorithmic implementations for multiscale modeling. By means of stochastic dynamics (KMC, random walk dynamics) algorithms, MBN Explorer allows simulating the dynamics of MBN systems on timescales significantly exceeding the limits for the conventional atomistic MD simulations. Such a multiscale dynamics approach is ideal for systems in which details of their atomistic dynamics become excessive, and the overall behavior of a system can be reproduced through kinetic rates for the dominating modes of motion and probabilities of the key processes occurring in the system. This important feature of MBN Explorer significantly expands its application areas and goes beyond the limits of other MD codes usually unable to deal with multiscale modeling.

Fig. 2

(a) Magnetic undulator for the European XFEL and (b) a Si$_{1-x}$Ge$_x$ superlattice Crystalline Undulator (CU) build atop the silicon substrate with the face normal to the $\langle 100 \rangle$ crystallographic direction [1]. The (110) planes in the superlattice are bent periodically. The picture of a CU is courtesy of J.L. Hansen, A. Nylandsted and U. Uggerhøj (University of Aarhus)

In the case of ultrarelativistic charged particle propagation through different media, the implemented algorithms enable simulations of particle dynamics on macroscopically large distances with atomistic accuracy. These algorithms enable to obtain the necessary atomistic insights into macroscopically large systems and processes occurring therein. For instance, one can simulate the operation of novel intensive sources of high-energy monochromatic $\gamma$-rays based on irradiation of oriented crystals by beams of ultra-relativistic electrons and positrons.

3.2 Main features of MBN Studio

In order to facilitate the practical work with MBN Explorer, a special multi-task software toolkit, called MBN Studio, has been developed [3]. It simplifies modeling of MBN systems, setting up and starting calculations with MBN Explorer, monitoring their progress and examining the calculation results. The software can be utilized for any type of calculations supported by MBN Explorer.

MBN Studio enables the Project set-up (standard as well as application-specific). Application-specific projects are usually designed for utilization in specific application areas, e.g., related to novel or emerging technologies. Such application-specific projects often involve special algorithms. A special modeling plug-in allows one to construct and prepare application-specific projects for simulation quickly and efficiently.

MBN Studio has an advanced MBN system modeler, a built-in tool for the computational design of MBN systems. By means of this plug-in, one can easily assemble molecular systems of different geometries and compositions for their further simulations with MBN Explorer.

MBN Studio supports various standard input/output data formats and links to numerous online databases and libraries with coordinates and geometries for atomic clusters, nanoparticles, biomolecules, crystals and other molecular systems, which can be utilized in simulations with MBN Explorer.

MBN Studio is equipped with the output data handling, visualization and analytic tools that allow calculation and analysis of specific characteristics determined by the output of MD simulations. Examples include calculations of diffusion coefficients, heat capacities, melting temperatures for solids, atomic radial distribution functions and many others.

MBN Studio also enables video rendering of the dynamics of MBN systems simulated with MBN Explorer.

Fig. 3

Schematic representation of a CU [23, 31, 32]. Closed circles mark the atoms of crystallographic planes which are periodically bent with the amplitude a and period $\lambda _{\text {u}}$. Thin dotted line illustrates the trajectory of the particle (open circles), which propagates along the center line (the undulator motion) and simultaneously undergoes so-called channeling oscillations. The periodic motion leads to the emission of the undulator-type radiation, and, under certain conditions, may result in the stimulated radiation

To summarize, MBN Explorer and MBN Studio are powerful tools for computational modeling in numerous challenging research areas arising in connection with the development of the aforementioned technologies. There are several such areas in which simulations performed by means of MBN Explorer and MBN Studio contributed immensely to their development. For instance, one of the areas concerns the development of novel light sources based on charged particles channeling in crystalline undulators [22,23,24]. Another example concerns controllable fabrication of nanostructures with nanometer resolution using tightly focused electron beams [25,26,27]. The third example deals with simulations of the nanoscopic molecular processes playing the key role in the ion-beam cancer therapy [13, 28, 29]. MBN Explorer combined with the visualization interface of MBN Studio in many cases can substitute expensive laboratory experiments by computational modeling, making the software play a role of a “computational nano- and microscope”. The following sections illustrate the capabilities of MBN Explorer and MBN Studio for advanced multiscale modeling in the three aforementioned technological areas.

4 Computational modeling of novel crystal-based gamma-ray light sources

MBN Explorer and MBN Studio provide unique possibilities for computational multiscale modeling of the physical processes and phenomena [23, 24, 30,31,32,33] that may facilitate the design and practical realization of novel gamma-ray crystal-based light sources (CLS) [22]. Such light sources can be constructed through the exposure of oriented crystals (linear, bent, periodically bent) to beams of ultra-relativistic charged particles (e.g., electrons or positrons). The construction of novel CLSs is a challenging task that constitutes a highly interdisciplinary field entangling a broad range of correlated activities. CLSs provide a low-cost alternative to conventional X-ray LSs based on free electron lasers (FEL) and have an enormous number of applications in the basic sciences, technology, and medicine [22].

The development of LSs for wavelengths $\lambda$ well below 1 Å (corresponding photon energies $E_{\mathrm{ph}} > rsim 10$ keV) is a challenging goal of modern physics. Sub-angstrom wavelength powerful spontaneous and, especially, coherent radiation will have many applications in nuclear and solid-state physics and life sciences. At present, several FEL facilities (FERMI, FLASH, LCLS, SACLA) provide femtosecond short laser-like photon pulses to user experiments. Their wavelengths range from the EUV and soft X-rays (FERMI, FLASH) to hard X-rays (LCLS, SACLA) [34, 35]. However, no laser system has yet been commissioned for wavelengths $\lambda$ significantly shorter than 1 Å due to the limitations of permanent magnet and accelerator technologies. Modern synchrotron facilities [36, 37] provide radiation of shorter wavelengths but of much less intensity which falls off very rapidly as $\lambda$ decreases.

Therefore, to create a powerful LS for $\lambda$ well below 1 Å, i.e., in the hard X- and gamma-ray band, one should consider new approaches and technologies. Novel gamma-ray CLS can generate radiation in the photon energy range where the technologies based on the charged particles’ motion in the fields of permanent magnets become inefficient or incapable. The limitations of conventional LSs are overcome by exploiting very strong crystalline fields that can be as high as $\sim$10$^{10}$ V/cm, which is equivalent to a magnetic field of 3000 T, whilst modern superconducting magnets provide 1–10 T [38]. The orientation of a crystal along the beam significantly enhances the strength of the particles’ interaction with the crystal due to strongly correlated scattering from lattice atoms. This enables the guided motion of particles through crystals of different geometries and the enhancement of radiation.

Practical realization of CLSs often relies on the channeling effect. The basic phenomenon of channeling is in a large propagation distance of a projectile particle that moves along a crystallographic plane or axis and experiences collective action of the electrostatic fields of the lattice atoms [39]. In the planar regime, positrons channel between two adjacent planes, whereas electrons propagate in the vicinity of a plane, thus experiencing more frequent collisions. As a result, a typical distance covered by a particle before it leaves the channeling mode due to uncorrelated collisions (the so-called dechanneling length) for positrons is much larger than for electrons. To ensure enhancement of the emitted radiation due to the dechanneling effect, the crystal length must be chosen of the order of the dechanneling length [30,31,32].

The projectile’s motion and the radiation emission in bent and periodically bent crystals (BCs and PBCs) are similar to those in magnet-based synchrotrons and undulators. In the latter, the particles and photons move in vacuum. In contrast, in crystals, the particles propagate in a medium, leading to several limitations for the crystal length, bending curvature, and beam energy. However, the crystalline fields are so strong that they steer ultra-relativistic particles more effectively than the most advanced superconducting magnets. The strong crystalline fields enable to reduce the bending radius in BC down to the cm range and bending period $\lambda _{\text {u}}$ in PBCs (see Fig. 3) to the hundred or even ten microns range. The size and the corresponding cost of such devices are orders of magnitude smaller than that for the analogous devices based on magnets [34]. Figure 2 shows the comparison of the magnetic undulator for the European XFEL with the CU manufactured in the University of Aarhus and used in recent experiments [40].

Fig. 4

a A schematic representation of the channeling process [33]. The z-axis is aligned with the incident beam direction ($\mathbf{v}_0$ denotes the initial velocity) and is parallel to the crystallographic direction, along which the channeling is simulated. The y-axis is perpendicular to the crystal plane. At the entrance, the x and y coordinates of the particle are randomly chosen to lie in the central part of the (xy)-plane (the highlighted rectangle). b Illustration of the trajectory of an axially channeling negatively charged particle

One of the important practical realizations of the CLSs is based on the concept of a crystalline undulator (CU) [23, 30, 31]. A CU device contains a PBC and a beam of ultra-relativistic positrons or electrons undergoing planar channeling, see Fig. 3. In such a system, in addition to the channeling radiation (ChR) [41], the undulator radiation appears due to the periodic motion of the particles, which follow the bending of the planes. A light source based on a CU can generate photons in the energy range from tens of keV up to the GeV region [42]. (The corresponding wavelengths range starts at 0.1 Å and goes down to $10^{-6}$ Å.) The intensity and characteristic frequencies of the CU radiation (CUR) can be varied by changing the beam energy, the parameters of bending and the type of crystal. Under certain conditions, a CU can become a source of the hard X- and gamma-ray laser light within the range $\lambda = 10^{-2} - 10^{-1}$ Å [23], which cannot be reached in existing and planned FELs based on magnets.

The mechanism of the photon emission by means of CU is illustrated by Fig. 3, which presents a cross section of a single crystal. The z axis is aligned with the mid-plane of two neighboring non-deformed crystallographic planes (not shown) spaced by the distance d. The closed circles denote the nuclei of the planes, which are periodically bent with the amplitude a and period $\lambda _{\text {u}}$. The harmonic (sinusoidal) shape of periodic bending, $y(z) = a \sin {(2\pi z/\lambda _\mathrm{u})}$, is of a particular interest since it results in a specific pattern of the spectral-angular distribution of the radiation emitted by a beam of ultra-relativistic charged particles (the open circles in the figure) propagating in the crystal following the periodic bending.

The operational principle of a CU does not depend on the projectile type. Provided certain conditions are met, a particle undergoes channeling in a PBC [23, 30,31,32]. Its trajectory contains two elements. First, there are channeling oscillations due to the action of the interplanar force. Their frequency $\Omega _{\mathrm{ch}}$ depends on the projectile energy $\varepsilon$, on the maximal value of the interplanar force and the interplanar distance d. Second, there are undulator oscillations due to periodicity of the bending whose frequency is $\Omega _{\mathrm{u}} \approx 2\pi c/\lambda _{\mathrm{u}}$. The spontaneous emission is associated with both of these oscillations. The typical frequency of ChR is $\omega _{\mathrm{ch}} \approx 2 \gamma ^2 \Omega _{\mathrm{ch}}$, where $\gamma = (1 - v^2/c^2)^{-1/2}$ is the relativistic Lorenz factor with v being the speed of a projectile and c the speed of light. The frequency of CUR is $\omega _{\text {u}} \approx 2 \gamma ^2 \Omega _{\text {u}}$. If $\Omega _{\text {u}} \ll \Omega _{\mathrm{ch}}$, then the ChR and CUR frequencies are well separated. In this case, the characteristics of CUR are virtually independent on the channeling oscillations [23, 30, 31], and the operational principle of a CU is the same as of a magnet-based one (e.g., [43,44,45]) in which the monochromaticity of radiation is due to constructive interference of the photons emitted from similar parts of trajectory.

The feasibility for the construction of a PBC-based light source, a CU, was verified theoretically relatively recently [30,31,32, 46, 47]. In the cited papers, as well as in the subsequent publications (see the latest review [23]), the essential conditions and limitations which must be met were formulated. These papers gave rise to coordinated theoretical, computational, technological, and experimental studies of several related phenomena.

MBN Explorer enables simulations of propagation of various (positively and negatively charged, light and heavy) particles in various media, such as hetero-crystalline structures (including superlattices), bent and periodically bent crystals, amorphous solids, liquids, nanotubes, fullerites, biological environment, and many more. The applicability of the code to different structures can be adjusted either by choosing a proper interaction potential or, if necessary, by including a new potential [33].

By these means, the channeling phenomenon [23, 39], which takes place when a charged particle enters a crystal at small angles with respect to a crystallographic direction, can be modeled [24]. Depending on the orientation, one can distinguish planar and axial channeling when a projectile enters the crystal at small angles with respect to a crystallographic plane or axis, respectively. The particle becomes confined and forced to move through the crystal, preferably along the crystallographic direction, experiencing collective action of the electrostatic field of the lattice ions (see Fig. 4). Since the field is repulsive for positively charged particles, they are steered into the interatomic region, while negatively charged projectiles move in close vicinity of ion strings or planes. Figure 4b illustrates the axial channeling of a negatively charged projectile.

4.1 Relativistic equations of motion

In order to simulate motion of relativistic particles, MBN Explorer considers [33] relativistic equations of motion:

$$\begin{aligned} \left\{ \begin{array}{ll} {\dot{\mathbf{v}}} = {\frac{1}{m \gamma } \, \left( \mathbf{F} - \mathbf{v} \frac{\mathbf{F } \cdot \mathbf{v }}{c^2} \right) }, &{} \\ {\dot{\mathbf{r}}} = \mathbf{v } &{} \end{array} \right. \end{aligned}$$

(1)

where $\gamma$ is the relativistic Lorenz factor. The force $\mathbf{F} = -\mathbf{\nabla } U(\mathbf{r})$ acting on the projectile is due to the interaction with the crystal constituents. This interaction is characterized by the electrostatic potential energy $U(\mathbf{r})$ which is considered as a sum of atomic potentials $U_\mathrm{at}$:

$$\begin{aligned} U(\mathbf{r}) = \sum _j U_\mathrm{at}(|\mathrm{r} - \mathrm{R}_j|), \end{aligned}$$

(2)

where $\mathrm{R}_j$ is the position vector of the jth atom. In MBN Explorer, the interaction between the charged projectiles and the constituent atoms is described by means of the widely used Molière approximation as well as the more recent approximation by Pacios [23]. Both interaction potentials are written as a product of the Coulomb potential and a screening function $\chi (r_{ij})$, see Refs. [1, 23, 24] for details.

Applied to the propagation in a medium, equations (1) describe the classical motion of a particle in the electrostatic field of the medium atoms. As a first step in simulating the projectile’s motion along a particular crystallographic direction, the crystalline structure inside the simulation box of the size $L_x \times L_y \times L_z$ is generated. The z-axis is aligned with the chosen crystallographic direction, see Fig. 4a. The integration of equations of motion (1) starts at $t=0$ when the particle enters the crystal at $z=0$. The initial coordinates $x_0$ and $y_0$ are chosen randomly within the central part of the (xy)-plane, shown by the red rectangle in Fig. 4a. The initial velocity $\mathbf{v}_0 = (v_{0_x}, v_{0_y}, v_{0_z})$ is predominantly oriented along the z-axis. The transverse components $v_{0_x}, v_{0_y}$ can be chosen with account for the beam emittance and energy distribution of particles in the beam [33].

4.2 Dynamic simulation box

To simulate the propagation of particles through an infinite medium, MBN Explorer utilizes the dynamic boundary conditions. In this case, all particles in the system are separated into two groups: the moving particles and the fixed ones.

Fig. 5

Illustration of the dynamic simulation box algorithm [33]. Once an X-marked projectile approaches the simulation box boundary (left panel), a new box of the same size is generated (right panel) with the particle placed approximately in its center. Positions of atoms (small shadowed circles) located in the intersection of the old and new boxes do not change. In the rest part of the new box, the atomic positions are generated anew

For a single moving particle interacting with fixed atoms of the medium, the dynamic boundary algorithm is illustrated in Fig. 5. Once the particle approaches the border of the simulation box, a new box of the same size is generated with the particle placed at the geometrical center of the box. To avoid spurious change in the force acting on the projectile, the positions of the atoms located in the intersection of the two boxes are conserved. The remaining part of the new box is filled with atoms of the medium following the specifications in the input file.

The dynamic simulation box moves following the propagation of the particle. Thus, it allows the “on-the-fly” modeling of the crystalline structure in the course of integrating the relativistic equations of motion, Eq. (1). Atoms of the medium are generated in the vicinity of the projectile as many times as it is necessary to propagate the projectile through the crystal of the thickness L. Periodic boundary conditions will not be effective in this case as the trajectory of the simulated particle is not periodic at all but instead rather unique at every instant. Similarly, one can consider the motion of many projectiles propagating through the medium built of fixed constituents. In this case, the computational scheme illustrated in Fig. 5 is applied to each projectile.

The dynamic simulation box algorithm permits to study at the atomistic level of detail the channeling phenomenon in mesoscopically large crystals, being micron-to-cm in length (see Fig. 6 as an illustrative example). Systems of such size cannot be handled by means of the all-atom MD approach with the standard periodic boundary conditions. Further details of this implementation are given in Refs. [23, 33].

Fig. 6

Channeling of 6.7 GeV positrons (a) and electrons (b) in a 105-$\mu$m-thick silicon crystal. The plots show typical trajectories of the particles initially collimated along Si(110) crystallographic planes. The trajectories are simulated by means of the dedicated channeling module of MBN Explorer [33] employing the Molière potential. Horizontal dashed lines indicate the planes separated by the distance $d = 1.92$ Å

MBN Explorer enables to simulate classical trajectories for different types of projectiles and their energies, different crystals (straight, bent and periodically bent), and crystallographic directions. The latest version of the software enables accounting for the ionization energy losses and the radiation damping, as well as for the atom knockout processes.

Figure 6 illustrates representative trajectories of 6.7 GeV positrons (panel A) and electrons (panel B) propagating in a 105-$\mu$m-thick silicon crystal. The particles were initially collimated along Si(110) crystallographic planes. Horizontal dashed lines indicate the planes separated by the distance $d = 1.92$ Å, that is the (110) interplanar distance in the silicon crystal at $T = 300$ K. The trajectories shown in Fig. 6 have been selected from a bunch of simulated particle trajectories. In any particular case study, between $5 \times 10^3$ and $5 \times 10^4$ independent trajectories are typically considered. Each simulated trajectory describes the propagation of a single projectile through the medium. In the simulations, (i) transverse coordinates and velocities of the projectile at the crystal entrance and (ii) positions of the lattice atoms due to the thermal fluctuations are randomly sampled [33]. Therefore, each simulated trajectory corresponds to a unique crystalline environment and all simulated trajectories are statistically independent.

Figure 6 illustrates a difference in channeling motion experienced by electrons and positrons. Positrons (Fig. 6a) move along the (110) planes bouncing between two neighboring planes, and their oscillations exhibit nearly harmonic character. In contrast, electrons (Fig. 6b) propagate oscillating in the vicinity of a plane. As a result, their trajectories are strongly non-harmonic, and they tend to leave the channeling mode of motion earlier as compared to positrons.

Simulated trajectories can be used further for statistical characterization of the radiation emitted by a projectile. For each set of simulated trajectories of the total number $N_\mathrm{{traj}}$, MBN Explorer provides an option to calculate the spectral distribution of the energy $\mathrm{d} E$ emitted within the cone $\theta \le \theta _{0} \ll 1$ along the direction of the incident particle beam:

$$\begin{aligned} \left\langle \frac{\mathrm{d} E}{\mathrm{d} (\hbar \omega )} \right\rangle= & {} \frac{1}{N_\mathrm{{traj}}} \sum _{n=1}^{N_\mathrm{{traj}}} \frac{\mathrm{d} E_n }{ \mathrm{d} (\hbar \omega )} \ , \nonumber \\ \frac{\mathrm{d} E_n }{ \mathrm{d} (\hbar \omega )}= & {} \int \limits _{0}^{2\pi } \mathrm{d} \phi \int \limits _{0}^{\theta _{0}} \theta \mathrm{d}\theta \, \frac{\mathrm{d}^2 E_n }{ \mathrm{d} (\hbar \omega )\, \mathrm{d}\Omega }. \end{aligned}$$

(3)

Here, $\mathrm{d}^2 E_n / \mathrm{d} (\hbar \omega )\, \mathrm{d}\Omega$ stands for the spectral-angular distribution emitted by a projectile moving along the jth trajectory. In MBN Explorer, this quantity is calculated within the framework of the quasi-classical approximation [48] that combines the classical description of the particle’s motion with the quantum corrections due to the radiative recoil and using the algorithm described in Refs. [23, 33]. The sum in Eq. (3) is carried out over all simulated trajectories, i.e., its takes into account the contribution of the channeling segments of the trajectories as well as of those corresponding to the non-channeling regime.

Fig. 7

Radiation spectra from 6.7 GeV electrons (solid line) and positrons (dashed line) channeling through a 105-$\upmu$m-thick Si(110) crystal [33]. The lines are drawn over ca. 200 photon energy points in which the spectra were calculated. The symbols mark a small fraction of the points and are drawn to illustrate typical statistical errors (due to a finite number of the simulated trajectories) in different parts of the spectrum

Figure 7 shows the dependence $\mathrm{d}E/\mathrm{d}(\hbar \omega )$ calculated for 6.7 GeV electrons and positrons aligned along Si(110) crystallographic plane at the crystal entrance [33]. Statistical uncertainties due to the finite number of the analyzed trajectories ($\approx 500$ trajectories in each case) are indicated by the error bars which correspond to the probability $\alpha = 0.999$.

The spectral distribution of emitted radiation, Eq. (3), is calculated by means of the quasi-classical formalism from the simulated trajectories. This means that all the radiation mechanisms associated with particle dynamics that are present in the system are accounted for. The bremsstrahlung radiation is much weaker (by one-two orders of magnitude) than the channeling radiation providing particles move in the channeling regime [23]. The quasi-classical approximation [48] used in the simulations reduces to the Bethe–Heitler theory (used to describe the bremsstrahlung in the energy range considered) if particles leave channels and experience incoherent interactions with atoms of the target crystal. The enhancement of the radiation in the channeling regime arises due to the coherent interaction of ultrarelativistic particles with the oriented crystalline planes or axes. These are all very well-known effects discussed, for example, in Refs. [23, 24]. In the case of periodically bent crystals (as schematically shown in Fig. 3), due to the periodicity of the particles’ trajectories following the bending of the crystalline planes, the photon emission spectrum attains the prominent features of the undulator radiation, which are automatically captured by the quasi-classical formalism used in the simulations. These phenomena are described in detail in the book [23] and in the recent review papers [22, 24].

In connection with the problems of constructing the high-quality undulator material and carrying out a quantitative analysis of the structures manufactured by different methods, MBN Explorer enables all-atom MD simulations to model periodically bent structures under various external conditions. This approach, combined with modern numerical algorithms and advanced computational facilities, will bring the predictive power of the software up to the accuracy level comparable or higher than that achievable experimentally. It can turn computational modeling into a very useful tool, which in many cases could substitute expensive laboratory experiments by computational modeling and thus reduce the experimental and technological costs.

Verification of the channeling module of MBN Explorer against experimental data as well as predictions of other theoretical models has been carried out in a large number of publications [24, 33, 42, 49,50,51,52,53], where channeling of electrons and positions in straight, bent and periodically bent crystals was studied, and the radiation spectra were calculated and successfully compared to experimental data.

5 Computational modeling of the focused electron-beam-induced deposition process

The controllable fabrication of nanostructures with nanoscale resolution remains a considerable scientific and technological challenge [54]. To address such a challenge, novel techniques have been developed [26], which exploit irradiation of nanosystems with collimated electron and ion beams. One of such techniques, focused electron beam-induced deposition (FEBID) [25,26,27], is based on the irradiation of precursor molecules [55] by high-energy electrons, whilst they are being deposited on a substrate. Electron-induced decomposition releases the metallic component of the precursor, which forms a deposit on the surface with a size similar to that of the incident electron beam (typically, a few nanometers) [56].

FEBID enables reliable direct-write fabrication of complex, free-standing 3D structures [56, 57]. Still, as the intended resolution falls below 10 nm, even FEBID struggles to yield the desired size, shape and chemical composition [56, 58], which primarily originates from the lack of molecular-level understanding of the irradiation-driven chemistry (IDC) underlying nanostructure formation and growth [27, 56]. Computational multiscale modeling provides a new methodology for understanding IDC and, consequently, advancing controllable fabrication of nanostructures.

Fig. 8

Electron-induced fragmentation rates for W(CO)$_6$ precursor molecules irradiated with PE beams of $E_0 = 1$ keV (a) and $E_0 = 30$ keV (b). The green transparent surface depicts the PE beam area [59]

FEBID operates through successive cycles of organometallic precursor molecules replenishment on a substrate and irradiation by a focused electron beam, which induces the release of metal-free ligands and the growth of metal-enriched nanodeposits. It involves a complex interplay of phenomena that require dedicated computational approaches: (a) deposition, diffusion and desorption of precursor molecules on the substrate; (b) multiple scattering of the primary electrons (PE) through the substrate, with a fraction of them being reflected (backscattered electrons, BSE) and the generation of additional secondary electrons (SE) by ionization; (c) electron-induced dissociation of the deposited molecules; and (d) the follow-up chemistry along with potential thermomechanical effects. While processes (b) and (c) typically happen on the femtosecond-to-picosecond timescale, (a) and (d) may require up to microseconds or even longer.

Until recently, most computer simulations of FEBID and the nanostructure growth have been performed using a Monte Carlo (MC) approach and diffusion-reaction theory [26, 60, 61], which allow simulations of the average characteristics of the process concerning local growth rates and the nanostructure composition. However, these approaches do not provide any molecular-level details regarding structure (crystalline, amorphous, mixed) and the IDC involved. At the atomic/molecular level, ab initio methods permit the precise simulation of electronic transitions or chemical bond reorganization [62, 63], although their applicability is typically limited to the femtosecond–picosecond timescales and to relatively small molecular sizes. In between these approaches, classical MD [1] and particularly reactive MD [64] have proved to be very useful in the atomistic-scale analysis of molecular fragmentation and chemical reactions up to nanoseconds and microseconds [64, 65].

A breakthrough in the atomistic simulation of FEBID was achieved recently by means of irradiation-driven molecular dynamics (IDMD), a novel and general methodology for computer simulations of irradiation-driven transformations of complex molecular systems [66]. This approach overcomes the limitations of previously used computational methods and describes FEBID-based nanostructures at the atomistic level by accounting for chemical transformation of surface-adsorbed molecular systems under focused electron beam irradiation [1, 59, 66].

Within the IDMD framework, various quantum processes occurring in an irradiated system (e.g., ionization, bond dissociation via electron attachment, or charge transfer) are treated as random, fast, and local transformations incorporated into the classical MD framework in a stochastic manner with the probabilities elaborated on the basis of quantum mechanics [66]. Major transformations of irradiated molecular systems are simulated by means of MD with reactive CHARMM (rCHARMM) force field [64, 67] using the MBN Explorer [2] and MBN Studio [3] software packages.

In the pioneering study [66], IDMD was successfully applied for the simulation of FEBID of W(CO)$_6$ precursors on a SiO$_2$ surface and enabled to predict the morphology, molecular composition, and growth rate of tungsten-based nanostructures emerging on the surface during the FEBID process. The follow-up study [59] introduced a novel multiscale computational methodology that couples event-by-event MC simulations of electron transport [68, 69] with IDMD for simulating the IDC processes with atomistic resolution. The developed multiscale modeling approach enables simulation of radiation transport and effects in complex systems where all the FEBID-related processes (deposition, irradiation, replenishment) are considered. The spatial and energy distributions of secondary and backscattered electrons emitted from a SiO$_2$ substrate were used to simulate electron-induced formation and growth of metal nanostructures obtained after deposition of W(CO)$_6$ precursors on SiO$_2$.

Within the IDMD framework, the space-dependent rate for bond cleavage in molecules on the substrate surface is given by:

$$\begin{aligned} P(x,y)= & {} \sigma _{\mathrm{frag}}(E_0) \, J_{\mathrm{PE}}(x,y,E_0) \nonumber \\&+ \sum _i \sigma _{\mathrm{frag}}(E_i) \, J_{\mathrm{SE/BSE}}(x,y,E_i) \ , \end{aligned}$$

(4)

where $E_0$ is the initial energy of the electron beam, $E_i < E_0$ a discrete set of values for the electron energies lower than $E_0$; $J_{\mathrm{PE/SE/BSE}}(x, y, E_i)$ are space- and energy-dependent fluxes of PE/SE/BSE (electrons per unit area and unit time), and $\sigma _{\mathrm{frag}}(E_i)$ is the energy-dependent molecular fragmentation cross section. The PE beam flux at the irradiated circular spot of radius R is:

$$\begin{aligned} J_0 = \frac{I_0}{e \, S_0} \ , \end{aligned}$$

(5)

where $I_0$ corresponds to the PE beam current, $S_0 = \pi R^2$ to its area and e is the elementary charge. The electron distributions were simulated using the MC radiation transport code SEED (Secondary Electron Energy Deposition) [68, 69]. Molecular fragmentation and further chemical reactions were simulated by means of MBN Explorer [2], while MBN Studio [3] was employed for constructing the molecular system, performing the precursor molecule replenishment phases, as well as for analyzing the IDMD simulation results.

Fig. 9

Evolution of the number of atoms in the largest simulated islands for PE beams of energies 1, 10 and 30 keV, for different currents as indicated [59]

Figure 8 illustrates the space-dependent fragmentation rates induced by uniform 1-keV (panel A) and 30-keV (panel B) beams of unit PE flux $J_0 = 1$ nm$^{-2}$fs$^{-1}$ within a circular area of radius $R = 5$ nm. Although the number of BSE/SE electrons for 30 keV is small, their large cross section (in relation to PE) produces a significant fragmentation probability, but less than that due to PE at the beam area. However, for 1 keV, the fragmentation probability due to BSE/SE ( $\sim$ 80–90$\%$ exclusively due to SE) is very large and significantly extends beyond the PE beam area. These results demonstrate the very different scenarios to be expected for beams of different energies, which will importantly influence the deposit properties and the prominent role of low-energy SE on molecular fragmentation.

Each irradiation phase lasts for a time known as dwell time, whose typical duration in experiments ($\ge \mu$s) is still computationally demanding for MD. To address this challenge, the irradiation phase was simulated for 10 ns and simulated PE fluxes $J_0$ (and hence PE beam currents $I_0$) were then scaled to match the same number of PE per unit area and per dwell time as in experiments [66]. As for replenishment, its characteristic times are also typically very long ($\sim$ms). In simulations, the CO molecules desorbed to the gas phase are removed during the replenishment stages and new W(CO)$_6$ molecules are deposited.

As the irradiation-replenishment cycles proceed, atomic clusters and islands of different sizes and compositions appear on the substrate as the result of IDC, and the process of nucleation of metal-enriched islands and their coalescence starts [66]. Figure 9 shows the number of atoms (either W, C or O) in the largest island for three simulation conditions close to reported in experiments [70]: 30 keV at $I_0 = 0.28$ nA, 10 keV at $I_0 = 2.3$ nA and 1 keV at $I_0 = 3.7$ nA. Smaller clusters tend to merge with time, giving rise to larger structures. The jumps in the island size observed with some frequency are due to the merging of independent clusters that grow on the substrate.

Experimental measurements performed to date have been limited to particular values of energy and current due to the characteristics of the electron source [70]. In contrast, the IDMD simulation method permits the exploration of a much broader range of electron beam parameters. Full symbols in Fig. 10a depict the simulated metal contents of the deposits as a function of experimentally equivalent current $I_{\mathrm{exp}}$. Error bars show the standard deviations obtained from three independent simulations for each case. Experimental results [70] are shown by open symbols. Numbers next to symbols represent the PE beam energies in keV. It is clearly seen that the results from simulations are within the range of experimental uncertainties, which indicates the predictive capabilities of the simulations.

Fig. 10

Compositions and morphologies of the deposits created by FEBID [59]. (a) Dependence of the deposit metal content on the beam energy $E_0$ and current $I_{\text {exp}}$ from experiments (open symbols) [70] and simulations (full symbols) [59]. Numbers next to symbols represent the beam energy in keV for each case. Panels b and c show the top views of the deposits produced by 10keV@2.3nA and 1keV@3.7nA beams, respectively. The green area marks the PE beam spot, while blue, white, and red spheres represent, respectively, W, C and O atoms; the SiO$_2$ substrate is represented by a yellow surface

This analysis provides a detailed “map” of the attainable metal content in the deposits as a function of the beam parameters, which is a valuable outcome for optimizing FEBID with W(CO)$_6$ on SiO$_2$. Dashed lines in Fig. 10a correspond to the limiting values of PE beam energy and current studied. These results clearly show that, within the analyzed energy domain, a decrease in the beam energy and an increase in the current promote the faster growth of the deposit, as well as the augment in its metal content. Simulation results provide the grounds for clearly understanding such trends: an increment in the current means a larger number of PE per unit time, while a reduction in the energy produces an increase in the SE yield. These lead to both the greater size of the deposit and its larger metal content due to the increased probability for bond cleavage (see Fig. 8).

Figure 10b, c shows top views of the simulated deposits for 10keV@2.3nA and 1keV@3.7nA, after 7 and 5 irradiation cycles, respectively. (The number of atoms in the largest island is approx. 12000 in both cases.) The green circle marks the area covered by the PE beam (having a radius of 5 nm). These figures show that different energy-current regimes lead to distinct deposit microstructures and edge broadenings. While the more energetic 10 keV energy beam produces a deposit almost entirely localized within the PE beam area, the 1 keV beam produces a more sparse deposit (at least during the early stage of the FEBID process) that significantly extends beyond the PE beam area, producing an undesired edge broadening of the structure.

The results presented in this section demonstrate how the novel MC-IDMD approach provides the necessary molecular insights into the key processes behind FEBID, which can be used for its further optimization and development. The simulations [59, 66] have demonstrated a great predictive power, yielding fabricated nanostructure compositions and morphologies in good agreement with available experimental data [70]. The IDMD methodology provides a wide range of possibilities for atomistic-level study of FEBID and many other processes in which the irradiation of molecular systems and irradiation-driven chemistry play a key role.

6 Computational modeling of ion-induced DNA damage in relation to ion-beam cancer therapy

MBN Explorer and MBN Studio can be utilized to evaluate radiobiological damage created by heavy ions propagating in different media, including biological. Such an analysis is behind the important biomedical technology known as the ion-beam cancer therapy (IBCT) [13, 28, 29, 71]. IBCT allows delivery of high doses into tumors, maximizing cancer cell destruction and simultaneously minimizing the radiation damage of surrounding healthy tissue. The full potential of such therapy can only be realized if the fundamental mechanisms leading to lethal cell damage under ion irradiation are well understood. The key question is whether it is possible to quantitatively predict macroscopic biological effects caused by ion radiation on the basis of physical and chemical effects related to the ion-medium interactions on a nanometer scale.

Recent review papers [13, 72] and the book [29] presented an overview of the main ideas of the multiscale approach to the physics of radiation damage with ions (MSA). This approach has the goal of developing knowledge about biodamage at the nanoscale and molecular level and finding the relation between the characteristics of incident particles and the resultant biological damage. The MSA is unique in distinguishing essential phenomena relevant to radiation biodamage at a given time, space and energy scale, and assessing the damage [13, 29, 72]. Temporal and spatial scales are schematically shown in Fig. 11.

Fig. 11

Features, processes, and disciplines, associated with radiation therapy, shown in a space-time diagram, which indicates temporal and spatial scales of the phenomena [13]. The history from ionization/excitation to biological effects is shown in the main figure, and features of ion propagation are shown in the inset

Radiation damage due to ionizing radiation is initiated by the ions incident on tissue. Initially, they have energy ranging from a few to hundreds of MeV per nucleon. In the process of propagation through tissue, the ions lose energy due to ionization, excitation, nuclear fragmentation, etc. Most of the energy loss of the ion is transferred to the tissue. Naturally, radiation damage is associated with this transferred energy, and the dose (i.e., deposited energy density) is a common indicator for the assessment of the damage [13, 28, 73]. The profile of the linear energy transfer (LET), that is the energy absorbed by the medium per unit length of the projectile’s trajectory, along the ion’s path is characterized by a plateau followed by a sharp Bragg peak. The position of this peak depends on the initial energy of the ion and marks the location of the most severe damage in the target DNA molecule. In the process of radiation therapy, a tumor is being “scanned” with the Bragg peak both laterally and longitudinally.

Fig. 12

The scenario of biological damage with ions [13]. Ion propagation ends with a Bragg peak, shown in the top right corner. A segment of the track at the Bragg peak is shown in more detail. Secondary electrons and radicals propagate away from the ion’s path damaging biomolecules (central circle). They transfer the energy to the medium resulting in the rapid temperature and pressure increase inside the cylinder. The ion-induced shock wave (shown in the expanding cylinder) due to the pressure increase damages biomolecules by stress (left circle), but it also effectively propagates reactive species to larger distances (right circle). A living cell responds to all shown DNA damage by creating foci (visible in the stained cells), in which enzymes attempt to repair the induced lesions. If these efforts are unsuccessful, the cell dies (the lower right corner)

However, the deposition of large doses in the vicinity of the Bragg peak does not explain how the radiation damage occurs, since projectiles themselves only interact with a few biomolecules along their trajectory and this direct damage is only a small fraction of the overall damage. It is commonly understood that the secondary electrons and free radicals produced in the processes of ionization and excitation of the medium with ions are largely responsible for the vast portion of the biodamage.

Fig. 13

Survival probabilities as a function of deposited dose for different human (a) and Chinese hamster (b) cell lines. The survival probabilities calculated within the MSA [13] at the indicated values of LET are shown with lines. Experimental data for cells measured at a specific dose, either at normal or hypoxic conditions, are shown by symbols. See Refs. [74, 75] for further details

Secondary electrons are produced during a short time of $10^{-18}$–$10^{-17}$ s following the ion’s passage. The energy spectrum of secondary electrons has been extensively discussed in the literature, see, e.g., the review [13]. The main result is that most secondary electrons for ions at the Bragg peak region have energy below 50 eV. This has several important consequences. First, the ranges of propagation of these electrons in tissue are rather small, below 10 nm [76]. Second, the angular distribution of their velocities, as they are ejected from their original host and as they scatter further, is largely uniform [77]; this allows one to consider their transport using a random walk approach [78,79,80,81].

The next timescale $10^{-16} - 10^{-15}$ s corresponds to the propagation of secondary electrons in tissue. In liquid water, the mean free paths of elastically scattered and ionizing 50-eV electrons are about 0.43 and 3.5 nm, respectively [77]. This means that they ionize a molecule after about seven elastic collisions, while the probability of second ionization is small [82]. Thus, the secondary electrons are losing most of their energy within first 20 collisions, and this happens within 1–2 nm of the ion’s path [83]. After that they continue propagating, elastically scattering with the molecules of the medium until they get bound or solvated electrons are formed. It is important to notice that these low energy electrons remain important agents for biodamage since they can attach to biomolecules like DNA causing dissociation [84]. The solvated electrons may play an important role in the damage scenario as well [85, 86].

The energy lost by secondary electrons in the processes of ionization and excitation of the medium is transferred to its heating (i.e., vibrational excitation of molecules) due to the electron–phonon interaction. As a result, the medium within a 1–2-nm region (for ions not heavier than iron) surrounding the ion’s path is heated up rapidly [83, 87]. The pressure inside this narrow region increases by several orders of magnitude (e.g., by a factor of 10$^3$ for a carbon ion at the Bragg peak [87]) compared to the pressure in the medium outside that region. This pressure builds up by about $10^{-14}$–$10^{-13}$ s, and it is a source of a cylindrical shock wave [88], which propagates through the medium for about $10^{-13}$–$10^{-11}$ s. Its relevance to the biodamage is as follows. If the shock wave is strong enough (the strength depends on the distance from the ion’s path and the LET), it may inflict damage directly by breaking covalent bonds in a DNA molecule [67, 83, 89,90,91,92,93]. Besides, the radial collective motion of the medium induced by the shock wave is instrumental in propagating the highly reactive molecular species, such as hydroxyl radicals and solvated electrons, to large radial distances (up to tens of nanometers), thus increasing the area of an ion’s impact [93,94,95].

The assessment of the primary damage to DNA molecules and other parts of cells due to the above mentioned effects is done within the MSA. This damage consists of various lesions on DNA and other biomolecules. Some of these lesions may be repaired by the living system, but some may not and the latter may lead to cell inactivation.

The scenario described above is illustrated in Fig. 12 [13]. The detailed comparison of its outcomes with experimental observations on the survival probability of irradiated cells is performed in Refs. [74, 75]. Figure 13 highlights this comparison for several exemplar case studies.

6.1 Simulation of the thermomechanical damage of the DNA by the ion-induced shock wave

The direct thermomechanical damage of the DNA molecule as a result of interaction with the ion-induced shock wave has been explored using MBN Explorer [67, 83, 89, 90, 93]. In the earlier investigations [83, 89, 90], the DNA damage by ion-induced shock waves was studied by means of classical MD simulations using non-reactive molecular mechanics force fields. In those simulations, the potential energy stored in a particular DNA bond was monitored in time as the bond length varied around its equilibrium distance [83, 90]. When the potential energy of the bond exceeded a given threshold value, the bond was considered broken.

Fig. 14

a The cylindrical shock wave front in water (on the right; ion’s path is the axis of this cylinder, perpendicular to the figure plane) interacts with a segment of a DNA molecule on the surface of a nucleosome (on the left). The yellow dot indicates the place where damage occurs. The medium is very dense following the wave front and is rarefied in the wake. b The dependence of the logarithm of the normalized number of the covalent bond energy records for the selected DNA backbone region per 0.01 eV energy interval on the bond energy for four values of LET: 900, 1730, 4745, and 7195 eV/nm, corresponding to the Bragg peak values for carbon, neon, argon, and iron ions, respectively. Straight lines correspond to the fits of these distributions. The figures are adapted from Ref. [83]

In the pioneering study [83], the MD simulations were focused on the interaction of the cylindrical shock wave originating from ion’s path with a fragment of a DNA molecule situated on the surface of a nucleosome, see Fig. 14a. Nucleosomes, histone-protein octamers wrapped about with a DNA double helix, are the primary structural units of chromatin, which is a principal component of the cell nucleus in eukaryotic cells. The simulations [83] were done for four values of LET, namely 900, 1730, 4745, and 7195 eV/nm, corresponding to the Bragg peak values for carbon, neon, argon, and iron ions, respectively. Carbon ions are clinically used for cancer treatment, whereas heavier ions up to iron are present in galactic cosmic rays, being potentially damaging for humans during space missions [96, 97].

The simulations were performed using the standard CHARMM force field [98], which implies the harmonic approximation for describing the interaction potentials for covalent bonds and thus does not allow to observe bond breaking events directly. Therefore, to study whether the covalent bonds in the DNA backbone can be broken during the shock wave action, the energy temporarily deposited to these bonds was calculated. The analysis of MD simulations [83] performed for the four indicated values of LET gives the distributions of the bond energy records. These records can be represented by a histogram that assigns to every interval of energy ($\varepsilon , \varepsilon + \delta \varepsilon$) the number of records corresponding to the bond energies from this interval. For each value of LET, the bond energy distribution was constructed. These distributions (normalized to the total number of records $N_r$ for each value of LET) are shown in Fig. 14b, where $\ln {(1/N_r dN/d\varepsilon )}$ is plotted versus the corresponding energy interval. Next, the number of energy records of selected covalent bonds of the DNA backbone exceeding a given threshold was counted. This was done by direct counting of bond energy records for which $E > E_0$, where $E_0$ is a variable threshold. The records counted for $E_0 = 2.5$ eV are shown in Fig. 14b to the right of the dashed vertical line.

A more quantitative description of the ion-induced shock wave phenomenon has become possible [93] by means of reactive MD simulations that permit explicit simulation of covalent bond rupture and formation [64]. A recent study [67] presented a detailed computational protocol for modeling the shock wave induced DNA damage by means of the rCHARMM force field [64].

The target DNA molecule studied in [67, 93] contained 30 complementary DNA base pairs. The molecule was placed in a water box extending 17 nm from the DNA in the x- and y-directions and 8 nm in the z-direction. The total system size was 1,010,994 atoms.

Fig. 15

a Illustration of the propagation of the ion-induced shock wave in the vicinity of the DNA segment containing 30 base pairs [67]. Water molecules inside the 1-nm radius “hot” cylinder surrounding the ion track (shown by the yellow arrow) are highlighted. b The density of water in the radial direction from the ion’s path is shown at different instants from 0 to 12 ps after irradiation

It is widely established that one of the key events of radiation-induced DNA damage concerns the formation of single- and double-strand breaks (SSBs and DSBs) of the sugar-phosphate backbone. Therefore, the rCHARMM force field was used to describe interatomic interactions in the C$_3^{\prime }$–O, C$_4^{\prime }$–C$_5^{\prime }$, C$_5^{\prime }$–O and P–O bonds in the DNA backbone, which connect the sugar ring of one nucleotide and the phosphate group of an adjacent nucleotide. Bond dissociation energies and cutoff distances for bond breakage/formation were determined from density functional theory (DFT) calculations [67]. Covalent interactions in other parts of the DNA molecule were modeled using the standard CHARMM force field.

As described above, low-energy electrons with the average kinetic energy of about 40 eV, which are predominantly produced for ions at the Bragg peak region, lose most of their energy by ionizing and exciting molecules of the medium within approximately 1 nanometer from the ion’s path [13, 94]. Therefore, in the MD simulations the energy lost by the propagating ion has been deposited into the kinetic energy of water molecules located inside a “hot” cylinder of 1 nm radius around the ion’s path. The equilibrium velocities of all atoms inside that spatial region are increased by a factor $\alpha$ such that the kinetic energy of these atoms reads as [83, 89, 99]:

$$\begin{aligned} \sum _i^N \frac{1}{2} m_i \left( \alpha v_i \right) ^2 = \frac{3 N k_{B} T}{2} + S_e \, l \ . \end{aligned}$$

(6)

Here $S_e$ is the LET of the projectile ion, l is the length of the simulation box in the z-direction (parallel to the ion’s path), and N is the total number of atoms within the “hot” cylinder. The first term on the right-hand side of Eq. (6) is the kinetic energy of the 1-nm radius cylinder at the equilibrium temperature, $T = 300$ K, whereas the second term describes the energy loss by the ion as it propagates through the medium. Note that non-uniform radial distribution of deposited energy, which accounts for the transport of energetic $\delta$-electrons, was addressed in an earlier study [100]. The cited study demonstrated that the more accurate radial dose distribution (as compared to the uniform energy distribution within the cylinder of 1 nm radius) does not affect the shock wave dynamics for ions in the Bragg peak region in any significant way.

The shock wave propagates in the molecular system radially away from the ion track, see Fig. 15a. The range of the shock wave propagation in the aqueous environment can be determined by monitoring the radial density of the water molecules in time. The ion track was oriented in this case directly through the geometrical center of the DNA strand. The results of this analysis are shown in Fig. 15b. The wavefront moves toward the edge of the simulation box as time passes, and the wave profile becomes lower and broader, showing that the shock wave relaxes as time passes. This indicates that the impact of the shock wave weakens over time. Figure 15b shows that the maximal density of the wave does not change significantly in the range of 2–6 nm from the ion track; thus, a DNA strand placed in this range is expected to receive the strongest impact from the shock wave.

Fig. 16

The number of strand breaks in the 30-base pair-long DNA segment induced by the thermomechanical stress by the argon ion with an LET of 2890 eV/nm as a function of simulation time [67]. Two insets in the top show the breakage of the DNA strands at simulation time instances of 4 ps and 15 ps

The impact of a shock wave on DNA can be characterized through the probability of strand breaks formation, which can be used to quantify the amount of biodamage induced by the shock wave mechanism [13]. Figure 16 shows the number of strand breaks in the DNA segment caused by a shock wave induced by an argon ion with $S_e = 2890$ eV/nm, as a function of simulation time. The number of breaks rises quickly within the first 7 ps of the simulation, where the shock wave front hits both DNA strands of the target molecule. After this moment, the DNA damage rate becomes lower lasting until approximately 15 ps; hereafter, the number of breaks remains steady. Even if the number of breaks generally rises with time, some fluctuations of this quantity can be seen locally. This effect might be attributed to the broken bonds, which can be rejoined if the atoms involved get close to each other after the initial bond breakage. This analysis shows that multiple or even complex stand breaks might be induced by the generated shock waves.

The methodology reviewed in this section can be applied in further computational studies considering irradiation of DNA with different ions and different orientations between the ion’s path and the DNA molecule. Such analysis is important for understanding the radiation damage with ions on a quantitative level, focusing on particular physical, chemical, and biological effects that bring about lethal damage to cells exposed to ion beams [13, 74, 75].

7 Conclusions

MBN Explorer and MBN Studio are powerful tools for computational modeling in different areas of challenging research arising in connection with the development of the novel and emerging technologies. Illustrative case studies of multiscale modeling have been presented in this paper in relation to three emerging technologies, namely (i) development of novel sources of monochromatic high-energy radiation based on the crystalline undulators, (ii) controlled fabrication of nanostructures using the focused electron-beam-induced deposition, and (iii) ion-beam cancer therapy. These examples illustrate that the unique algorithms and methodologies implemented in MBN Explorer, combined with the visualization interface and additional tools of MBN Studio, in many cases can substitute expensive laboratory experiments by computational multiscale modeling, making the software play a role of a “computational nano- and microscope”.