An Algorithm for Computing Gradients of Functionals in Structural Materials Science Yuri G. Evtushenko Vladimir I. Zubov Alla F. Albu Federal Research Center Federal Research Center Federal Research Center ”Computer Science and ”Computer Science and ”Computer Science and Control” of RAS Control” of RAS Control” of RAS Moscow, Russia 119333 Moscow, Russia 119333 Moscow, Russia 119333 yuri-evtushenko@yandex.ru vladimir.zubov@mail.ru alla.albu@yandex.ru Abstract When describing and modeling the crystal structure of a material, in- teratomic interaction potentials are used. To solve the problems of parametric identification of potentials, which consists of the selection of potential parameters, the use of various optimization methods has recently become increasingly important. In this case, it becomes nec- essary to determine the gradients of some quantities characterizing the substance with respect to parameters of potential. Based on the Fast Automatic Differentiation-technique, an algorithm has been built to determine the exact values of the gradients of bulk modulus and shear modulus of the test substance with respect to potential parameters in the case when the total interatomic energy of the atoms’ system is determined using the Tersoff’s potential. 1 Introduction Recently, optimization problems from the field of structural materials science have become more relevant, among which the well-known problems of identifying the parameters of mathematical models stand out. Various mathematical models are used to study the atomic structures of materials. Some parameters of these models are unknown. They should be identified from the condition that the calculated properties of the modeled material are close to its properties, which were found experimentally. When describing and modeling the crystal structure of a material, interatomic interaction potentials are used. To describe the properties of crystals with a covalent bond (for example, carbon, silicon, germanium, etc.), Tersoff’s potential is often used ([Tersoff1988]). Structural identification of potentials for a specific crystalline material is one of the important stages of molecular dynamics modeling. The problem of parametric identification of potentials is often reduced to some optimization problem (see, for example, [Abgaryan2014]), which can be formally written as a problem of minimizing the following functional: ∑ m f (ξ) = ωi (yi (ξ) − yei )2 (1) i=1 Copyright ⃝ c by the paper’s authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0). In: Sergei S. Goncharov, Yuri G. Evtushenko (eds.): Proceedings of the Workshop on Applied Mathematics and Fundamental Computer Science 2021 , Omsk, Russia, 24-04-2021, published at http://ceur-ws.org 1 where ωi is the weight coefficient; yei is the value of the i-th material characteristic obtained experimentally, and yi (ξ) is the value of the same material characteristic calculated using Tersoff potential determined by the components of the vector ξ (ξ ∈ Rm is a vector to be identified). The solution to the problem is looked for on the set X ⊆ Rm , which is a parallelepiped. Its boundaries are chosen so that it obviously contained the admissible range of parameters. The number of items in the formula (1) varies depending on the studied material. For the numerical solution of this problem, gradient minimization methods are often used. There exists the need to efficiently calculate the exact value of the cost function gradient with respect to parameters of Tersoff potential. These derivatives are often calculated (in particular, see [Abgaryan2014]) using the finite difference method. Studies have shown that the finite difference method does not allow the calculation of the gradient of a cost function with acceptable accuracy and requires (m + 1) times to calculate the value of the function. One of the terms in the formula (1) is the total energy of the system of atoms. As the interatomic potential energy, the Tersoff potential was chosen. The other two terms in the formula (1) are bulk modulus and shear modulus of a material. In [Albu2016], using the Fast Automatic Differentiation-technique (FAD-technique) (see [Evtushenko1998]), formulas to calculate the exact gradient of the total energy with respect to parameters of Tersoff potential (specific for modeled substance) were received. In this paper, with the help of FAD-technique, we derived formulas to calculate the gradients of bulk modulus and shear modulus of material with respect to Tersoff parameters with machine precision. 2 Formulation of the Problem ∑ I ∑ I The energy E of a atoms’ system is calculated with the help of expression E = Vij , where Vij is the i=1 k=1 k ̸= i, j interaction potential between atoms marked i and j (i-atom and j-atom). In present paper the Tersoff potential is used as interaction potential: Vij = fc (rij ) (VR (rij ) − bij VA (rij )) , ( √ ) De ( √ ) SDe 2 R Vij = VR (rij ) = exp −β 2S(rij − re ) , A Vij = VA (rij ) = exp −β (rij − re ) , S−1 S−1 S    1,  1( ( )) r < R − Rcut , π(r − R) fc (r) = 1 − sin , R − Rcut < r < R + Rcut ,   2 2Rcut  0, r > R + Rcut , − 2η 1 ∑ I bij = (1 + (γζij )η ) , ζij = fc (rik )gijk ωijk , ωijk = exp(λ3 τijk ), k=1 k ̸= i, j ( c )2 c2 τijk = (rij − rik )3 , gijk = 1 + − . d d2 + (h − cos Θijk )2 Here I is the number of atoms in considered system; rij is the distance between i-atom and j-atom; Θijk is the angle between two vectors, first vector begins at i-atom and finishes at j-atom, second vector begins at i-atom and finishes at k-atom; R and Rcut are known parameters, identified from experimental geometric properties of substance. Tersoff Potential depends on ten parameters (m = 10), specific to modeled substances: De , re , β, S, η, γ, λ, c, d, h. The distance between i-atom and j-atom is determined by the formula: √ rij = (x1i − x1j )2 + (x2i − x2j )2 + (x3i − x3j )2 , where x1i , x2i , x3i are the Cartesian coordinates of i-atom. If Θijk is the angle between two vectors, connecting i-atom with j-atom and k-atom respectively, then 2 cos Θijk = qijk = (rij 2 + rik − rjk 2 )/(2rij rik ). 2 . Let a be the initial length of the edges of the lattice of atoms; e a = αa (α ∈ R) — length of( the edges of a+ρ ρ) the lattice of atoms after deformation; ρ = e a − a — deformation parameter. Then = 1+ a. If a a rk = (xk1 , xk2 , xk3 ) are the coordinates of some lattice atom before deformation and e ek2 , xek3 ) are its ( ρ) ( ρ) ( rk = ρ) (exk1 , x coordinates after deformation, then x ek1 = 1 + xk1 , x ek2 = 1 + ek3 = 1 + xk2 , x xk3 . a a a Let [( E(r , r , ..., r ) is the total energy of atoms system before deformation. Then E(e r1 , e r2 , ..., e rI ) = ρ) ( ρ) ( ρ) ] 1 2 I E 1+ r1 , 1 + r2 , ..., 1 + rI is the total energy of atoms system after deformation. The bulk a a a modulus and the shear modulus of the material are proportional to B(E) , that can be calculated by formula: ∂2 [( ρ) ( ρ) ( ρ) ] B(E) = E 1 + r 1 , 1 + r 2 , ..., 1 + rI . ∂ρ2 a a a ρ=0 Therefore, to calculate the gradients of the bulk modulus and shear modulus of the material with respect to Tersoff parameters, it is sufficient to calculate the gradient of the function B(E) with respect to these parameters. 3 Algorithm for Calculating the Gradient of a Function B(E) Let us construct the multistep algorithm to calculate the total energy E of atoms’ system (interaction potential is Tersoff Potential). For compactness further in the study we introduce vectors u and z having the following T T coordinates: uT = [u1 , u2 , ..., u10 ] , z T = [z1 , z2 , ..., z10 ] , where u1 = De , u2 = re , u3 = β, u4 = S, u5 = η, u6 = γ, u7 = λ, u8 = c, u9 = d, u10 = h; { √ } z1 = z1ijk = (x1i − x1k )2 + (x2i − x2k )2 + (x3i − x3k )2 ≡ F (1, Z1 , U1 ), { √ } z2 = z2ijk = (x1j − x1k )2 + (x2j − x2k )2 + (x3j − x3k )2 ≡ F (2, Z2 , U2 ), { } ij 2 ijk (z13 ) + (z1ijk )2 − (z2ijk )2 z3 = z3 = qijk = ≡ F (3, Z3 , U3 ), 2z1ijk z13 ij { } z4 = z4ijk = fc (z1ijk ) ≡ F (4, Z4 , U4 ), { ( )2 } ijk u8 (u8 )2 z5 = z5 = gijk = 1 + − ≡ F (5, Z5 , U5 ), u9 (u9 )2 + (u10 − z3ijk )2 { } z6 = z6ijk = τijk = (z13 ij − z1ijk )3 ≡ F (6, Z6 , U6 ), { } z7 = z7ijk = ωijk = exp((u7 )3 z6ijk ) ≡ F (7, Z7 , U7 ), { } z8 = z8ijk = fc (rik )gijk ωijk = z4ijk z5ijk z7ijk ) ≡ F (8, Z8 , U8 ), { ∑I } z9 = z9ij = ζij = k=1;k̸=i,j z8ijk ≡ F (9, Z9 , U9 ), { } ij z10 = z10 = γζij = u6 z9ij ≡ F (10, Z10 , U10 ), { } ij z11 = z11 = (γζij )η = (z10 )u5 ≡ F (11, Z11 , U11 ), { } ij ij − 2u1 z12 = z12 = bij = (1 + z11 ) 5 ≡ F (12, Z12 , U12 ), { √ } ij z13 = z13 = (x1i − x1j )2 + (x2i − x2j )2 + (x3i − x3j )2 ≡ F (13, Z13 , U13 ), { ( √ )} ij ij z14 = z14 = VijR = u4u−1 1 exp −u3 2u4 (z13 − u2 ) ≡ F (14, Z14 , U14 ), { ( √ )} ij ij z15 = z15 = VijA = uu41−1 u4 exp −u3 u24 (z13 − u2 ) ≡ F (15, Z15 , U15 ), { } ij ij z16 = z16 = fc (z13 ) ≡ F (16, Z16 , U16 ), { } ij ij ij ij ij z17 = z17 = Vij = z16 (z14 − z12 z15 ) ≡ F (17, Z17 , U17 ), ( ) i = 1, I, j = 1, I, j ̸= i, k = 1, I, k ̸= i, j . 3 The energy E of the atoms’ system with the help of new variables may be rewritten as follows: E(z(u)) = ∑ I ∑I ij z17 (u). Variables z1 , z2 , ..., z17 (the phase variables) are determined by the specified above multistep i=1 j=1 j ̸= i algorithm zl = F (l, Zl , Ul ), (l = 17), where Zl is the set of elements zn in the right part of the equation zl = F (l, Zl , Ul ), and Ul is the set of elements un that appear in the right side of this equation. Note that each component zl depends on a number of other components (zlij or zlijk ). Let us introduce also the following designations: ze1 , ze2 , ..., ze17 and e ze1 , e ze2 , ..., e ze17 , where    (( ρ ) ( ρ ) ( ρ ) )  ∂z ijk (e r , e r , ij ik jk e r ) ∂z ijk n 1 + rij , 1 + rik , 1 + r jk zen = zenijk : zenijk = n = a a a , n = 1, 8,  ∂ρ ρ=0 ∂ρ  ρ=0  (( ) )   ∂znij (e rij ) ∂znij 1 + aρ rij  zen = zenij : zenij = = , n = 9, 17,  ∂ρ ρ=0 ∂ρ  ρ=0  (( ) ( ) ( ) )   ijk ijk rij , reik , rejk ) ∂ 2 znijk (e ∂ 2 znijk 1 + aρ rij , 1 + aρ rik , 1 + aρ rjk  e e e zen = zen : zen = = , n = 1, 8,  ∂ρ2 ρ=0 ∂ρ2  ρ=0  (( ) )   ij ij ∂ 2 znij (e rij ) ∂ 2 znij 1 + aρ rij  e e e zen = zen : zen = = , n = 9, 17,  ∂ρ2 ρ=0 ∂ρ2  ρ=0 ( ) i = 1, I, j = 1, I, j ̸= i, k = 1, I, k ̸= i, j . The above values are calculated by the formulas: (( ) ) ∂fc 1 + aρ rik ze1ijk = rik /a; ze2ijk = rjk /a; ze3ijk = 0; ze4ijk = ; ∂ρ ρ=0 ze5ijk = 0; ze6ijk = 3z6ijk /a; ze7ijk = 3z6ijk z7ijk (u7 )3 /a; ∑ I ze8ijk = z5ijk (e z4ijk z7ijk + z4ijk ze7ijk ); ze9ij = ze8ijk ; ij ze10 = ze9ij u6 ; k=1 k ̸= i, j ij ij ij u5 −1 ij 1 ij ij − 2u1 −1 ij ij ze11 = ze10 u5 (z10 ) ; ze12 = − u5 ze11 (1 + z11 ) 5 ; ze13 = z13 /a; 2 √ √ (( ) ) ij ij ij u3 2u4 z13 z14 ij ij ij u3 2/u4 z13 z15 ij ∂fc 1 + aρ rij ze14 =− ; ze15 =− ; ze16 = ; a a ∂ρ ρ=0 ij ij ij ij ij ij ij ij ij ij ij ij ij ij ze17 = ze16 z14 − ze16 z12 z15 + ze14 z16 − ze12 z16 z15 − ze15 z16 z12 ; (( ) ) e ijk ijk ijk ijk ijk ∂ fc 2 1 + aρ rik ze1 = e ze2 = e ze3 = e ze5 = 0; e ze4 = ; ∂ρ2 ρ=0 ijk ijk e ze6 = 6z6ijk /a2 ; e 3 ijk ijk ijk ze7 = a2 z6 z7 (u7 ) (3z6 (u7 )3 + 2); 3 ijk ijk ijk ij ∑ I ijk e ze8 = z5ijk (e z4ijk ze7ijk + z4ijk e ze4 z7ijk + 2e ze7 ); e ze9 = e ze8 ; k=1 k ̸= i, j 4 ij ij ij ij e ze10 = e ze9 u6 ; e ze11 = u5 (u5 − 1)(eij 2 ij u5 −2 z10 ) (z10 ) + u5 (e ij u5 −1 ze10 )(z10 ) ; e ij 1 + 2u5 ij 2 ij − 2u 1 −2 1 eij ij − 2u 1 −1 ij ze12 = (e z11 ) (1 + z11 ) 5 − ze11 (1 + z11 ) 5 ; ze13 = 0; 4u5 2u5 (( ρ) ) ij ij 2 ij 2(u3 )2 u4 (z13 ) z14 ij ij 2 ij 2(u3 )2 (z13 ) z15 ij ∂ 2 fc 1+ rij e ze14 = ; e ze15 = ; e ze16 = a ; 2 a 2 a u4 ∂ρ2 ρ=0 ij ij ij ij ij ij e ze17 = e ze16 z14 + 2e z16 ze14 − e ij ij ze16 z12 z15 − 2eij ij ij z16 ze12 z15 − 2eij ij ij z16 z12 ze15 + ij ij ij ij ij ij ij ij +e ze14 z16 −e ze12 z16 z15 − 2e z15 z16 ze12 − e ij ij ij ze15 z16 z12 . To compute the second derivative of a function fc (r) there is a need for smoothing this function. It is proposed to replace the function fc (r) as follows:    0, r ≥ R + Rcut ,  1, r ≤ R − Rcut , fc (r) =   C · (f ( ∗) φ(r) , ) R ≤ r < R + Rcut ,  C · 2f∗ − (f∗ )ψ(r) , R − Rcut < r ≤ R, 2 2 1 3 Rcut Rcut where C = , f∗ = exp(− ), φ(r) = , ψ(r) = . 2f∗ 2 (r − R − Rcut ) 2 (r − R + Rcut )2 ∑ I ∑ I e ij Thus, B(E) is calculated by the formula B(E) = ze17 (u), where the variables z1 , z2 , ..., z17 , i=1 j=1 j ̸= i ze1 , ze2 , ..., ze17 , e ze1 , e ze2 , ..., e ze17 are determined by mentioned above multistep algorithm. According to the formulas of FAD-technique, the conjugate variables corresponding to the phase variables z1 , z2 , ..., z17 , ze1 , ze2 , ..., ze17 , eze1 , e ze2 , ..., e ze17 are determined by the formulas: ijk ijk ijk pijk e z4ijk ze7ijk + z4ijk e e4 z7ijk + 2e ze7 )e z4ijk z7ijk + z4ijk ze7ijk )e pe8 + (e pijk ijk ijk ijk 5 = (z 8 + z4 z7 p8 ; e ijk ijk 3 ijk ijk 3 pijk 7 =ze4 z5ijk e pe8 + ze4ijk z5ijk peijk ijk ijk ijk 8 + z4 z5 p8 + 2 z6 (u7 )3 (3z6ijk (u7 )3 + 2)e pe7 + z6ijk (u7 )3 peijk 7 ; a a [ ] ij u5 −3 ij 2 ij u5 −2 eij eij pijk ij 8 = p9 ; pij ij 9 = u6 p10 ; pij 10 = u5 (u5 − 1)(u5 − 2)(z10 ) z10 ) + u5 (u5 − 1)(z10 (e ) ze10 pe11 + ij u5 −2 ij ij ij u5 −1 ij +u5 (u5 − 1)(z10 ) ze10 pe11 + u5 (z10 ) p11 ; (1 + 4u5 )(1 + 2u5 ) ij − 2u 1 −3 ij 2 eij (1 + 2u5 ) ij − 2u 1 −2 pij 11 = − (1 + z11 ) 5 (z11 ) pe12 + (1 + z11 ) 5 × 8(u5 )3 4(u5 )2 ij ij (1 + 2u5 ) 1 ×e ze11e ij − 2u −2 ij ij ij − 2u −1 ij 1 1 pe12 + (1 + z11 ) 5 ze11 pe12 − (1 + z11 ) 5 p12 ; 4(u5 )2 2u5 e ij ij ijk eij eij e ij ij ij 2 2(u3 )2 u4 (z13 ) eij √ ij z13 pij e16 z15 12 = (−z − 2eij ij z15 ze16 − z16 ze15 )pe17 ; pij e16e 14 = z pe17 + e p14 − u3 2u 4 peij ; a2 a 14 e ij ij ij eij eij 2(u3 )2 (z13 ) eij √ ij 2 ij u3 z13 pij 15 = (− e z z 16 12 − 2e z ij ij e z 12 16 − z e z ) 16 12 17 e p + e p15 − 2/u 4 peij 15 ; a2 u4 a ij ij ij ij eij eij pij e e14 − e ze12 z15 − 2eij ij ze15 − z12 ze15 )pe17 ; pij 16 = (z z12 17 = 0; ijk peijk z4ijk z5ijk e 7 = 2e pe8 + z4ijk z5ijk peijk 8 ; peijk 8 =peij 9 ; peij eij 9 = u6 p 10 ; ij u5 −2 e ij ij u5 −1 ij ij peij 10 = 2u5 (u5 − 1)(z10 ) pe11 + u5 (z10 ) ze10 pe11 ; 5 ij −1/(2u5 )−1 (1 + 2u5 ) ij −1/(2u5 )−2 ij eij (1 + z11 ) ij ij e ij peij 11 = 2 (1 + z11 ) e z e p 11 12 − peij 12 ; peij 12 = −2(e ij ij z16 z15 + ze15 z16 )pe17 ; 2(u5 ) 2u5 ij eij ij ij eij ij ij eij peij 14 = 2e z16 pe17 ; peij 15 = −2(e ij ij z16 z12 + ze12 z16 )pe17 ; peij 16 = 2(e ij z14 ij ij − ze12 z15 − ze15 z12 )pe17 ; peij 17 = 0; ijk ijk ijk ij ij ij ij ij u5 −1 eij e pe7 = z4ijk z5ijk e pe8 ; e pe8 = e pe9 ; e pe9 = u6e pe10 ; e pe10 = u5 (z10 ) pe11 ; ij −1/(2u5 )−1 e ij (1 + z11 ) e ij e ij ij ij eij pe11 = − pe12 ; pe12 = −z15 z16 pe17 ; 2u5 ij ij eij ij ij ij eij ij ij ij eij ij e pe14 = z16 pe17 e pe15 = −z12 z16 pe17 ; e ij pe16 = (z14 − z12 z15 )pe17 ; e pe17 = 1; ij ij ij The adjoint variables are calculate in the following order: e pe17 , e pe16 , ..., e pe7 , peij eij 17 , ..., p ij ij 7 , p17 , ..., p5 . Those adjoint variables, which formulas for calculation aren’t provided above, aren’t used for calculation of the components of the gradient. According to the formulas of FAD-technique, the partial derivatives of a complex function Ω(u) = B(E(u)) = ∑I ∑I e ij ze17 (u) with respect to independent variables um , (m = 1, 10) are determined by the relations: i=1 j=1 j ̸= i ( ) ∂Ω ∑I ∑ I ij z14 ij ij z15 ij = p + p ; ∂u1 i=1 u1 14 u1 15 j=1 j ̸= i ∂Ω ∑I ∑ I ( √ √ ) ij = z14 u3 2u4 pij 14 + z ij 15 u 3 2/u4 pij 15 ; ∂u2 i=1 j=1 j ̸= i ∂Ω ∑ I ∑ I ( ( √ ) ( √ ) ) ij ij = z14 − 2u4 (z13 − u2 ) pij 14 + z ij 15 − 2/u4 (z ij 13 − u2 ) p ij 15 + ∂u3 i=1 j=1 j ̸= i ( √ ) ∑ I ∑ I ij 2 ij 4u3 u4 (z13 ) z14 eij ij ij 2u4 z13 z14 ij + pe14 − pe14 + i=1 a2 a j=1 j ̸= i ( ) ∑ I ∑ I ij 2 ij 4u3 (z13 ) z15 eij ij ij √ z13 z15 + e p 15 − 2/u4 peij 15 ; i=1 a2 u4 a j=1 j ̸= i (( ) ) ∂Ω ∑I ∑ I z ij √ ij ij = − 14 − 0.5u3 2/u4 (z13 − u2 )z14 pij + ∂u4 i=1 u4 − 1 14 j=1 j ̸= i (( ) ) ∑ I ∑ I ij z15 √ ij ij + − + 0.5(u3 /u4 ) 2/u4 (z13 − u2 )z15 pij + i=1 u4 (u4 − 1) 15 j=1 j ̸= i 6 ( ( ) ) ∑ I ∑ I ( u √ ) ij 2 2(u3 )2 (z13 ) ij eij 3 ij ij + − 2/u4 z13 z14 peij 14 + z14 pe14 + i=1 2a a2 j=1 j ̸= i (( ) ( ) ) ∑ I ∑ I u3 √ ij ij ij ij 2 2(u3 )2 (z13 ) ij eij + 2/u4 z13 z15 pe14 − z15 pe15 ; i=1 2au4 a2 (u4 )2 j=1 j ̸= i ( ) ∂Ω ∑I ∑ I ij u5 ij ij −1/(2u5 ) (1 + z11 ) = (z10 ) ln(z10 )pij 11 + ij ln(1 + z11 )pij 12 + ∂u5 i=1 2(u5 )2 j=1 j ̸= i ∑ I ∑ I (( ) ) ij u5 −1 ij ij u5 −1 ij ij + (z10 ) ze10 + u5 (z10 ) ze10 ln(z10 ) peij 11 + i=1 j=1 j ̸= i ∑ I ∑ I (( ) ) ij u5 −2 ij u5 −2 ij ij 2 eij + (2u5 − 1)(z10 ) + ((u5 )2 − u5 )(z10 ) ln(z10 ) (e z10 ) pe11 + i=1 j=1 j ̸= i ∑ I ∑ I (( ) ij ij ) + ij u5 −1 (z10 ) ij u5 −1 + u5 (z10 ) ij ln(z10 ) e ze10e pe11 + i=1 j=1 j ̸= i (( ) ) ∑ I ∑ I ij −1/(2u5 )−1 (1 + z11 ) ij ln(1 + z11 ) ij −1/(2u5 )−1 ij ij + 2 − 3 (1 + z11 ) ze11 pe12 + i=1 2(u5 ) 4(u5 ) j=1 j ̸= i ∑ I ∑ I (( ) ) (1 + u5 ) ij −1/(2u5 )−2 ij 2 eij + − (1 + z11 ) z11 ) pe12 + (e i=1 2(u5 )3 j=1 j ̸= i (( ) ) ∑I ∑I ij (1 + 2u5 ) ln(1 + z11 ) ij −1/(2u5 )−2 ij 2 eij + (1 + z11 ) z11 ) pe12 + (e i=1 8(u5 )4 j=1 j ̸= i (( ) ) ∑I ∑ I (1 + z11ij −1/(2u5 )−1 ) ln(1 + z11 ij ) ij −1/(2u5 )−1 eij eij + − (1 + z11 ) ze11 pe12 ; i=1 2(u5 )2 4(u5 )3 j=1 j ̸= i ∂Ω ∑I ∑ I ( ij ) = z9ij pij e9ij peij eeij e e10 ; 10 + z 10 + z 9 p ∂u6 i=1 j=1 j ̸= i ∑I ∑ I ∑ I ( ) ∂Ω ijk ijk 2 ijk 9 ijk ijk 2 ijk = 3z6 z7 (u7 ) p7 + 6 z7 (u7 ) pe7 z + ∂u7 i=1 a j=1 k=1 j ̸= i k ̸= i, j 7 (( ) ) ∑ ∑ ∑ 3z6ijk z7ijk ( ) ijk I I I + ijk 18z6 (u7 )5 + 6(u7 )2 e pe7 ; i=1 a2 j=1 k=1 j ̸= i k ̸= i, j (( ) ) ∂Ω ∑ ∑ I I ∑ I 2u8 2u8 = − pijk 5 ; ∂u8 i=1 (u9 )2 (u9 )2 + (u10 − z3ijk )2 j=1 k=1 j ̸= i k ̸= i, j    ∂Ω ∑I ∑ I ∑I  −2(u8 ) 2 2 2(u8 ) u9  ijk  =  +( )2  p5  ; ∂u9 i=1 (u9 )3 ijk j=1 k=1 (u9 )2 + (u10 − z3 )2 j ̸= i k ̸= i, j    ∂Ω ∑ I ∑ I ∑I  ijk 2(u8 ) (u10 − z3 ) 2  ijk  =  ( )2  p5  . ∂u10 ijk i=1 j=1 k=1 (u9 )2 + (u10 − z3 )2 j ̸= i k ̸= i, j The received formulas for calculation of the gradient of function B(E(u)) outwardly are represented quite difficult and bulky. Therefore, there is a natural question: whether to use simpler approaches, for example, finite difference method, to calculate the gradient functions B(E(u)). In [Albu2016] the comparison of function gradients calculated by the finite differences and by using Fast Automatic Differentiation formulas was presented. The results of comparison are the following: 1) when computing the gradient of complicated function using finite differences, one must conduct researches related to the choice of suitable increments of each parameter; 2) for different parameters, the researches must be carried out independently; 3) for the same parameter, the researches must be carried out if its value changed; 4) to calculate the gradient of complicated function using finite differences one must (m + 1) times calculate the value of function itself. In contrary to it, the FAD-technique enables us to calculate gradients of any complicated function with the machine accuracy for arbitrary parameters. The machine time that is needed to calculate the gradient does not exceed three times of calculation of the function itself. 3.0.1 Acknowledgements The research was supported by Russian Science Foundation (project No. 21-71- 30005). References [Tersoff1988] J. Tersoff. Empirical interatomic potential for silicon with improved elastic properties. Phys. Rev. B., 38:9902–9905, 1988. [Abgaryan2014] K.K. Abgaryan and M.A. Posypkin. Optimization methods as applied to parametric identifica- tion of interatomic potentials. Comput. Math. Math. Phys., 54(12):1929–1935, 2014. [Albu2016] A.F. Albu. Application of the Fast Automatic Differentiation to the Computation of the Gradient of the Tersoff Potential. Informacionnye tekhnologii i vychislitel’nye sistemy, 1:43–49, 2016. [in Russian] [Evtushenko1998] Y.G. Evtushenko. Computation of Exact Gradients in Distributed Dynamic Systems. Opti- mizat. Methods and Software, 9:45–75, 1998. 8