**4. Computation of an Interval Stress Tensor for Materials with a Covalent Chemical Bond**

Let us consider an applied problem of computational materials science, within the framework of which the stresses arising during the deformation of an ideal crystal are calculated [12]. Different angles are possible between the orientation of the crystal lattice and the direction of stretching with a fixed stretch value. The stress tensor thus becomes interval. This problem is solved using molecular dynamics simulation. The motion of atoms is described by the classical equations of dynamics:

$$\begin{cases} \mathbf{r}\_{i}^{\prime} = \mathbf{v}\_{i\prime} \\ \mathbf{v}\_{i}^{\prime} = \frac{1}{m\_{i}} \mathbf{F}\_{i\prime} \end{cases}$$

where **r***<sup>i</sup>* is the radius vector of the atom with the number *i*, **v***<sup>i</sup>* is its velocity, *mi* is its mass, and **F***<sup>i</sup>* is the force acting on the atom, in this case **F***<sup>i</sup>* = −∇*iE*, where *E* is the total energy of the system, and ∇*<sup>i</sup>* is the gradient along the position of the atom with the number *i*.

In this problem, materials with a covalent interatomic bond are considered. The total energy of interaction between atoms of such materials is well described using the Tersoff potential [25]:

$$\begin{cases} \begin{aligned} &E = \frac{1}{2} \sum\_{i} \sum\_{j \neq i} V\_{ij}, \; V\_{ij} = f\_{\mathbb{C}}(r\_{ij}) \left( f\_{R} \left( r\_{ij} \right) + b\_{ij} f\_{A} \left( r\_{ij} \right) \right), \\ &f\_{\mathbb{C}}(r) = \begin{cases} &1, r < R - D, \\ \frac{1}{2} \left( 1 - \sin \left( \frac{\pi}{2} \frac{r - R}{D} \right) \right), R - D \le r < R + D, \\ &0, R + D \le r, \\ &f\_{\mathbb{R}}(r) = A \exp(-\lambda\_{1} r), \; f\_{A}(r) = -B \exp(-\lambda\_{2} r), \\ &b\_{ij} = \left( 1 + \beta^{n} \xi\_{ij}^{n} \right)^{-\frac{1}{2n}}, \\ \xi\_{ij} = \sum\_{k \neq i, j} f\_{\mathbb{C}}(r\_{ik}) \xi \left( \theta\_{ijk} \right) \exp\left( \lambda\_{3}^{m} \left( r\_{ij} - r\_{ik} \right)^{m} \right), \\ &g(\theta) = \gamma \left( 1 + \frac{c^{2}}{d^{2}} - \frac{c^{2}}{d^{2} + \left( \cos(\theta) - \cos(\theta\_{0}) \right)^{2}} \right), \end{aligned} \end{cases}$$

where *E* is total system energy, *Vij* is the contribution to the interaction energy of atoms with numbers *i* and *j*, *rij* is the distance between atoms with numbers *i* and *j*, *fC*(*r*) is a cut-off function, *fR*(*r*) and *fA*(*r*) are the repulsion and attraction potentials, respectively, and *R*, *D*, *A*, *B*, *n*, *m*, *λ*1, *λ*2, *λ*3, *β*, *γ*, *c*, *d* and cos(*θ*0) are potential parameters that are selected in order to reproduce the properties of the simulated material. Methods of parametric identification of the Tersoff potential parameters are considered in papers [26,27].

The initial positions of atoms and their number are determined by the structure of the crystal lattice and the restriction on the minimum size of the simulated space is specified by the structure of the potential.

Consider crystalline silicon as a typical material. The simulated sample is represented by eight unit cells of a diamond crystal lattice making up a cube of 2 × 2 × 2 in size, with periodic boundary conditions; each unit cell contains eight unique atoms (Figure 14), so a total of 64 atoms are involved in the simulation. Initial speeds are considered to be zero. The initial conditions for a dynamical system can be represented as follows:

$$\begin{cases} \mathbf{r}\_i(0) = \left( (\mathbf{x}, y, z)^T + (dx, dy, dz)^T \right) a\_\star(\mathbf{x}, y, z) \in \text{Base}, dx, dy, dz \in \{0, 1\}, \mathbf{v}\_i(0) = (0, 0, 0)^T, \\\quad \text{Base} = \left\{ \begin{array}{l} (0, 0, 0), (0, 0.5, 0.5), (0.5, 0, 0.5), (0.5, 0.5, 0), (0.75, 0.75, 0.75), \\ (0.75, 0.25, 0.25), (0.25, 0.75, 0.25), (0.25, 0.25, 0.75) \end{array} \right\}\_{\prime} \end{cases}$$

where *<sup>a</sup>* <sup>=</sup> 5.429 <sup>×</sup> <sup>10</sup>−<sup>10</sup> m is the linear size of a unit cell of a silicon crystal.

**Figure 14.** Unit cell of a silicon crystal.

To take into account deformation, 4 additional variables are introduced: one of them reflects the elongation value, and three more are responsible for the angle between the orientation of the lattice and the direction of stretching. In this case, the elongation is set to be fixed, and the variables responsible for the rotation are taken as interval. Rotation is generated evenly using quaternions [28].

The final system looks like this:

⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ **r***i* = **v***i*, **v***<sup>i</sup>* <sup>=</sup> <sup>−</sup> <sup>1</sup> 2*mi* ∑ *i* ∑ *j*=*i* ∇*i fC rij fR rij* + *bij fA rij*, *rij*(*s*, *μ*1, *μ*2, *μ*3) = -**S**(*s*)**R**(*μ*1, *μ*2, *μ*3)**r***ij*-, **S**(*s*) = diag(1 + *s*, 1, 1), **R**(*μ*1, *μ*2, *μ*3) = ⎛ ⎝ 1 − 2 *q*2 <sup>2</sup> + *<sup>q</sup>*<sup>2</sup> 3 2(*q*1*q*<sup>2</sup> − *q*3*q*0) 2(*q*1*q*<sup>3</sup> + *q*2*q*0) 2(*q*1*q*<sup>2</sup> + *q*3*q*0) 1 − 2 *q*2 <sup>1</sup> + *<sup>q</sup>*<sup>2</sup> 3 2(*q*2*q*<sup>3</sup> − *q*1*q*0) 2(*q*1*q*<sup>3</sup> − *q*2*q*0) 2(*q*2*q*<sup>3</sup> + *q*1*q*0) 1 − 2 *q*2 <sup>1</sup> + *<sup>q</sup>*<sup>2</sup> 2 ⎞ ⎠, (*q*0, *q*1, *q*2, *q*3) = 1 − *μ*<sup>1</sup> sin(2*πμ*2), 1 − *μ*<sup>1</sup> cos(2*πμ*2), √*μ*<sup>1</sup> sin(2*πμ*3), √*μ*<sup>1</sup> sin(2*πμ*3) , *t* ∈ B 0, 10−12C , (13)

where *<sup>m</sup>* <sup>=</sup> 4.65 <sup>×</sup> <sup>10</sup>−<sup>26</sup> kg is the mass of atoms, *<sup>s</sup>* <sup>=</sup> 0.1 is the relative elongation of the sample, *μ*<sup>1</sup> ∈ [0.1, 0.9], *μ*<sup>2</sup> ∈ [0.1, 0.9], and *μ*<sup>3</sup> ∈ [0.1, 0.9] are stretching direction parameters, *<sup>R</sup>* <sup>=</sup> 2.85 <sup>×</sup> <sup>10</sup>−<sup>10</sup> m, *<sup>D</sup>* <sup>=</sup> 0.15 <sup>×</sup> <sup>10</sup>−<sup>10</sup> m, *<sup>A</sup>* <sup>=</sup> 6.12 <sup>×</sup> <sup>10</sup>−<sup>16</sup> J, *<sup>B</sup>* <sup>=</sup> 1.81 <sup>×</sup> <sup>10</sup>−<sup>17</sup> J, *<sup>c</sup>* <sup>=</sup> 9.69, *<sup>d</sup>* <sup>=</sup> 2.35, *<sup>n</sup>* <sup>=</sup> 4.16, *<sup>β</sup>* <sup>=</sup> 0.132, *<sup>λ</sup>*<sup>1</sup> <sup>=</sup> 3.36 <sup>×</sup> <sup>10</sup><sup>10</sup> <sup>m</sup>−1, *<sup>λ</sup>*<sup>2</sup> <sup>=</sup> 1.27 <sup>×</sup> <sup>10</sup><sup>10</sup> <sup>m</sup>−1, *<sup>λ</sup>*<sup>3</sup> <sup>=</sup> 1.19 <sup>×</sup> <sup>10</sup><sup>10</sup> <sup>m</sup>−1, *<sup>γ</sup>* <sup>=</sup> 5.71, and cos(*θ*0) = <sup>−</sup>0.408 are the parameters of the potential.

Integration of the resulting ordinary differential system (13) was carried out using the Verlet velocity method with a constant integration step of 10−<sup>15</sup> s. As a result, the interval stress tensor was obtained (values are given in Pascals):

$$
\begin{pmatrix}
[-1.58 \times 10^{10} \text{J} - 1.43 \times 10^{10}] & [-1.35 \times 10^9, 1.35 \times 10^9] & [-1.35 \times 10^9, 1.35 \times 10^9] \\
[-1.35 \times 10^9, 1.35 \times 10^9] & [-4.79 \times 10^9, -1.46 \times 10^9] & [-1.42 \times 10^9, 1.42 \times 10^9] \\
[-1.35 \times 10^9, 1.35 \times 10^9] & [-1.42 \times 10^9, 1.42 \times 10^9] & [-4.79 \times 10^9, -1.46 \times 10^9]
\end{pmatrix}
$$

For *<sup>ε</sup>* <sup>=</sup> <sup>10</sup>−2, the obtained result was *<sup>I</sup>* <sup>=</sup> 1079132.3 and *error* <sup>=</sup> <sup>2</sup> <sup>×</sup> <sup>10</sup>−1.

Note that the possibilities of the proposed approach are not limited to a specific type of interatomic interaction potential in a material. The method can be applied to the simulation of the stress–strain state of materials with various types of chemical bonds, including the modeling of composite materials.

## **5. Discussion**

In the previous sections, the proposed approach was tested on representative interval problems of nonlinear dynamics and computational materials science. It is found that, thanks to sparse grids, it is possible to integrate ODE systems with a large number of interval uncertainties in a reasonable time. To estimate the computational costs, a criterion was used that is equal to the equivalent number of sampling points from the original uncertainty region.

Table 1 shows estimates of the computational costs for the ODE system given by Equation (9) describing a conservative oscillator with two interval initial conditions for two values *ε*. For *ε* = 10−3, the approach proposed in the paper works 1.5 to 5 times faster than the classical adaptive interpolation algorithm.

Table 2 shows the computational costs when integrating the ODE system given by Equation (10), describing the Lotka-Volterra model with two interval initial conditions and one interval parameter. Here, for *ε* = 10−3, the use of adaptive sparse grids gives an acceleration of 1.9 <sup>−</sup> 2.8 times compared to the classical algorithm, and for *<sup>ε</sup>* <sup>=</sup> <sup>10</sup>−5, the acceleration is from 1.25 time to 15 time.

For ODE systems given by Equations (11) and (12), it was possible to obtain a solution in a reasonable time only using the approach proposed in the paper since the number of interval uncertainties is quite large. To solve system given by Equation (11), the equivalent number of sampling points was about 80 thousand, and in the case of using the classical algorithm with a degree of polynomial 4, the value would be about 10 million. Thus, using sparse grids in this problem gives an acceleration of at least 125 times.

The tables show that adaptive sparse grids work faster than regular sparse grids, and even faster than full grids. This fact is in line with the theoretical foundations. The classical adaptive interpolation algorithm in example (9) with two interval uncertainties showed itself slightly better only when *ε* = 10−<sup>5</sup> and *p* = 4. This is primarily due to the chosen value of p. It is known that the greater the degree of the interpolation polynomial, the faster the error decreases with increasing mesh density. Therefore, it seems promising to use sparse grids on a nonlinear basis. We should also note the possibilities of applying adaptive grids to ODE systems not only with interval uncertainties but also with stochastic uncertainties, including applied nonlinear systems.
