Introduction

Conducting theoretical studies of the formation of clusters of point defects in semiconductor structures is an important task on the way to improving technologies for obtaining new materials. One of the reasons for the formation of various defects in the semiconductor structure, including point, extended, their clusters and complexes is the radiative forcing. As is known from a number of experimental studies [1,2], as a result of irradiation, such structures acquire various properties that can be used in the manufacture of new generation devices for nanoelectronics. In this work, the multiscale approach is used to calculate defects that form in semiconductor structures as a result of radiation exposure [3,4]. It is based on:

-selection of the main scale levels, including atomic-crystalline and molecular-dynamic; -development and application for solving the set tasks of mathematical models and algorithms specific for each level;

-combining models and algorithms in the general computing process. Due to the large labor-intensiveness and versatility of the methods used for numerical simulation of semiconductor materials, a software package was created for multiscale modeling of their structural features, which allowed studying the formation of clusters of point and extended defects in a computational experiment. The paper presents the results of calculations obtained using the software for paralleling the computations placed on the hybrid high-performance computing complex of the FIC "Informatics and Control" of the Russian Academy of Sciences.

At present, a theoretical study of the processes occurring in complex defect structures is a very urgent task. One of the common approaches to conducting such research, which provides a compromise between the speed of calculations and the accuracy of the results, is the method of molecular-dynamic calculations. However, the problem of molecular dynamic modeling of complex defective structures remains resource-intensive due to the need to consider systems of high dimensionality and to take into account the interaction of a huge number of atoms with each other. In this regard, it is necessary to develop algorithms that allow accelerating such calculations without sacrificing the size of the simulated atomic structures.

The task of developing a new computational algorithm for molecular dynamics calculations was solved in the framework of the study of clusters of point defects in single-crystal silicon. Crystal defects call any violation of the translational symmetry of a crystal -the ideal periodicity of the crystal lattice. There are several types of defects: point, linear, flat, bulk. Point defects are a local violation of the crystal structure, the dimensions of which in all three dimensions are comparable with one or several interatomic distances. These are defects associated with impurities, with displacement or with the replacement of a small group of atoms. Such defects usually occur during heating during crystal growth, during radiation exposure, and also as a result of the addition of impurities. One type of point defects is a vacancy -a free, unoccupied atom, a lattice site.

The occurrence of radiation defects is an inevitable side effect of such a method of modifying materials as ion implantation. Radiation defects occur due to exposure to the material of neutrons or gamma rays. Such an effect is characteristic of materials found in a nuclear reactor. In space, the impact on materials of electrons and protons, as well as heavy ions with low energy, is characteristic.

Theoretical and experimental studies show that defects in silicon can form complex extended structures. For numerical simulation of such structures and the study of their stability, it is necessary to take into account in the calculations a large number of interacting atoms. This, in turn, leads to an increase in the computational complexity of the problem and an increase in the time required to solve it. The use of new approaches to the parallelization of molecular dynamic calculations on high-performance computer clusters allows us to solve such problems today.

In this paper, a new efficient approach is proposed for parallel computations when solving the problem of the molecular dynamic description of the structure of silicon with interacting vacancy defects.

Multiscale Modeling of Clusters of Point Defects in Silicon

Multiscale approach was applied to calculations of ordered cluster configurations of vacancies and the interstitial atoms in Si [3.4]. Two large-scale levels (apart from zero) -atomic and crystal and molecular and dynamic were selected. Within set-theoretic representations it can be set by means of multiscale composition in which are involved: 𝐾 0,1,2 (0,14;1,1;1,2;2,1;2,2) = 𝐾 0,1,2

(𝑆𝑖)(1)

In Figure 1 shows the structure of a multiscale composition for calculating ordered cluster configurations of vacancies and interstitial atoms in Si. Specimens of basic compositions and the sequence of their use in the computational process are indicated. At the first level, data on the chemical composition and atomic crystal structure of Si (diamond structure), obtained using the basic composition 𝑪 1 1 ("CRYSTAL CHEMICAL FORMULA"), were used. Further, they were used as input data in the base composition 𝑪 1 2 ("QUANTUM-MECHANICAL CELL") during the first-principle calculations in the framework of the electron density functional theory using the VASP software package [5]. In the first-principle modeling of the structure of ideal silicon, a periodic cell consisting of 64 atoms of dimension (2 × 2 × 2) was used.

Refined the atomic-crystalline and electronic structure of silicon with defects, was calculated 𝐸 𝑐𝑜ℎ . The computational resources of the Interdepartmental Supercomputer Center of the Russian Academy of Sciences and Moscow State University of M. V. Lomonosov were used for the calculations.

At the second scale level, questions were studied of the time variation of the structure of silicon with defects and with defective clusters. A composition of computational models consisting of basic compositions 𝑪 2 1 ("NUCLEAR CLUSTER -STATIC") and 𝑪 2 2 ("NUCLEAR CLUSTER -DYNAMICS") was used. Moreover, when forming the input data 𝑪 2 1 we used the results of first-principle calculations [5], obtained using the basic composition 𝑪 1 2 . They were used as reference. The cohesive energy of the system was the global parameter transmitted from the first to the second scale 𝐸 𝑐𝑜ℎ .

Model Description

For the numerical simulation of silicon single crystal, an elementary silicon cell consisting of 8 atoms was used (Figure 2). Further, using parallel transfers, the unit cell multiplied to silicon structures containing 616, 1160, and 4504 atoms. A structure was simulated with a complex of point defects, two vacancies located in neighboring lattice sites (Figure 3), with a selected frequency of defect recurrencethrough a cell. The Solution of the Problem Molecular dynamics modeling was carried out using the multiparameter potential of Tersoff [6], which proved well in solving problems of modeling compounds with covalent bonds. The calculation of the cohesive energy of the system of atoms was carried out as follows:

𝐸 = ∑ 𝐸 𝑖 = 1 2 𝑖 ∑ 𝑉 𝑖𝑗 𝑖(7)𝑏 𝑖𝑗 = 𝜒 𝑖𝑗 (1 + 𝛾 𝑖 𝑛𝑖 𝜁 𝑖𝑗 𝑛𝑖 ) − 𝑛𝑖 2(11)

Potential parameters were selected as a result of solving the optimization problem [7]:

𝐹(𝜉) =

where 𝜉 = (𝜉 1 … 𝜉 𝑚 ) -parameters of Tersoff potential 𝐸 𝑐𝑜ℎ 𝑓𝑝𝑐 , 𝑎 𝑓𝑝𝑐 , 𝐵 𝑓𝑝𝑐 , 𝐶′ 𝑓𝑝𝑐 , 𝐶 44 𝑓𝑝𝑐 , 𝜁 𝑓𝑝𝑐 -reference values obtained using ab initio calculations and from publication [8]. Selected values of the potential parameters are given in Table 1.

where 𝑛𝜖[1, 𝑁]; 𝑚 𝑛mass of the n-th atom, 𝑛𝜖[1, 𝑁];

𝐹 𝑛force acting on the particle n.

To integrate the Cauchy problem, we used the Verlet velocity method [9]:

𝑟 𝑛 𝑘+1 = 𝑟 𝑛 𝑘 + 𝜏 𝑘 𝑣 𝑛 𝑘 − 𝜏 𝑘 2 2 𝜕𝑈 𝑛 𝑘 𝜕𝑟 𝑛 𝑘 (15)𝑣 𝑛 𝑘+1 = 𝑣 𝑛 𝑘 + 𝜏 𝑘 2 ( 𝜕𝑈 𝑛 𝑘+1 𝜕𝑟 𝑛 𝑘+1 + 𝜕𝑈 𝑛 𝑘 𝜕𝑟 𝑛 𝑘 ) (16)

This method is a compromise between accuracy and speed of implementation. The method is stable and accurate, as well as self-starting due to the fact that its velocity is taken into account to obtain the next position of the particle. 512 GB) Features of the implementation are associated with the use of the potential of Tersoff, which is difficult to parallelize due to its complex structure, and is also resource-intensive from the point of view of computation.

To implement the parallelized algorithm, the OpenMP standard was used. To do this, we determined the maximum available number of threads for the program instance being started, after which the input data were divided into the number of blocks equal to the number of available threads. Each input block was launched in its own stream. During the molecular dynamics simulation, at the end of each time step, the blocks exchanged data in order to synchronize the parameters of atoms that are common to different blocks (Figure 5). This approach has reduced the computational complexity by reducing the number of atoms processed in each stream. In turn, this made it possible to significantly accelerate the process of modeling, thanks to which it became possible to take into account a larger number of atoms in the calculations.

For a visual comparison of the speed of the parallel and sequential algorithms, the modeling process was limited to 10 time steps. The results of the calculations are presented in Table 2. According to the obtained results, the application of the developed algorithm on a personal computer allows you to significantly speed up the calculations, and the demonstrated result of acceleration on a supercomputer is more significant than on a quad-core processor. The results obtained using the developed software package were also compared with the results obtained using the LAMMPS package. A separate program was written to obtain results within the LAMMPS package. In it, the elementary cell of silicon was modeled from 16 atoms with the help of the Tersoff potential at different temperatures. The obtained data were compared with a similar experiment set in the developed software package. The bond energy (E_pair) and system temperature (Temp) were compared. From the comparison performed, it can be seen that the binding energy is close to the energy values obtained using the LAMMPS software, and the temperature is comparable to the results obtained in LAMMPS. The implemented algorithm made it possible to calculate the characteristics of monocrystalline silicon, in particular, the value of the binding energy for monocrystalline silicon was calculated (𝐸 𝑐𝑜ℎ (𝜉) = −4.6314е𝑉). The obtained value turned out to be close to the value obtained earlier with the help of ab initio calculations (𝐸 𝑐𝑜ℎ 𝑓𝑝𝑐 = −4.6305е𝑉) [10]. Thus, the developed software using the selected parameters of the Tersof potential allows one to fairly accurately describe the geometric and energy properties of monocrystalline silicon. This approach can be further applied to the study of more complex structures of vacancy clusters in a silicon single crystal.

Conclusion

In this paper, software was developed that enables molecular dynamics modeling, effectively parallelizing computational flows. As part of the task of studying the stability of defective clusters in monocrystalline silicon, the developed software was used to calculate the divacancy binding energy in the structure of monocrystalline silicon. The proposed algorithm allows us to significantly accelerate molecular dynamics calculations, making it possible to take into account a larger number of interacting atoms in the simulation. In turn, this will allow us to study the properties of more complex defective structures in silicon. This approach can be applied further to simulate interacting atomic systems described by different potentials.

The calculations were performed by Hybrid high-performance computing cluster of FRC CS RAS) [19,20] and Shared Facility Center "Data Center of FEB RAS" (Khabarovsk) [21].

Application

The designation of variables used in formulas.

𝑬 -total energy of the system, (eV); 𝑽 𝒊𝒋 -potential energy of the interaction of two particles i and j, (eV); 𝒇 𝑪 (𝒓 𝒊𝒋 ) -cutoff-function 𝒓 𝒊𝒋 -distance between two particles i and j, (Å); 𝑹 𝒄𝒖𝒕 -cutoff distance (Å); 𝒇 𝑹repulsive potential between two atoms, (eV); 𝒇 𝑨potential of attraction between two atoms, (eV); В 𝒊𝒋 -parameter of attraction between two atoms, (eV);