When describing and modeling the crystal structure of a material, interatomic 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 necessary 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.
Introduction (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
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 )) , bij = (1 + (γζij )η)− 21 , ζij = fc(rik)gijkωijk,
ωijk = exp(λ3τijk),
De S − 1
R − Rcut < r < R + Rcut, ViRj = VR(rij ) = exp (
√ ) −β 2S(rij − re) ,
ViAj = VA(rij ) =
SDe exp S − 1 ( −β √ 2
S
) (rij − re) , τijk = (rij − rik)3, gijk = 1 + (c )2 d
c2 − 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,
I ∑ k = 1 k ̸= i, j 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
cos Θijk = qijk = (ri2j + ri2k − rj2k)/(2rij rik).
Let a be the initial length of the edges of the lattice of atoms; ea = αa (α ∈ R) — length of the edges of the lattice of atoms after deformation; ρ = ea − a — deformation parameter. Then a + ρ = (1 + ρ )a. If a a rk = (xk1, xk2, xk3) are the coordinates of some lattice atom before deformation and rek = (xek1, xek2, xek3) are its coordinates after deformation, then xek1 = (1 + aρ )xk1, xek2 = (1 + aρ )xk2, xek3 = (1 + aρ )xk3.
Let E(r1, r2, ..., rI ) is the total energy of atoms system before deformation. Then E(er1, er2, ..., erI ) = E [(1 + aρ )r1, (1 + aρ )r2, ..., (1 + aρ )rI ] is the total energy of atoms system after deformation. The bulk modulus and the shear modulus of the material are proportional to B(E) , that can be calculated by formula: B(E) = ∂2 ∂ρ2
E [(1 + ρ ) a
( r1, 1 + ρ ) a
( r2, ..., 1 + ρ ) a rI ] ρ=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 coordinates: uT = [u1, u2, ..., u10]T , zT = [z1, z2, ..., z10]T , where u1 = De, u2 = re, u3 = β, u4 = S, u5 = η, u6 = γ, u7 = λ, u8 = c, u9 = d, u10 = h; z1 = {zijk = √(x1i − x1k)2 + (x2i − x2k)2 + (x3i − x3k)2}
1 z2 = {zijk = √(x1j − x1k)2 + (x2j − x2k)2 + (x3j − x3k)2} 2 ≡ F (1, Z1, U1), ≡ F (2, Z2, U2), (u8 )2 u9 −
(u8)2 (u9)2 + (u10 − z3ijk)2 }
≡ F (5, Z5, U5), ≡ F (3, Z3, U3), { { z3 = zijk = qijk = 3 z4 = {z4ijk = fc(z1ijk)} ≡ F (4, Z4, U4), (z1ij3)2 + (z1ijk)2 − (z2ijk)2 }
2z1ijkz1ij3 z5 = zijk = gijk = 1 + 5 z6 = {z6ijk = τijk = (z1ij3 − z1ijk)3} ≡ F (6, Z6, U6), z7 = {z7ijk = ωijk = exp((u7)3z6ijk)} ≡ F (7, Z7, U7), z8 = {z8ijk = fc(rik)gijkωijk = z4ijkz5ijkz7ijk)} ≡ F (8, Z8, U8), z9 = {z9ij = ζij = ∑Ik=1;k̸=i,j z8ijk} ≡ F (9, Z9, U9), z10 = {z1ij0 = γζij = u6z9ij } ≡ F (10, Z10, U10), z11 = {z1ij1 = (γζij )η = (z10)u5 } ≡ F (11, Z11, U11), z12 = {z1ij2 = bij = (1 + z1ij1)− 2u15 } ≡ F (12, Z12, U12), z13 = {z1ij3 = √(x1i − x1j )2 + (x2i − x2j )2 + (x3i − x3j )2 } z14 = {z1ij4 = ViRj = u4u−11 exp z15 = {z1ij5 = ViAj = uu41−u41 exp ( ( −u3√2u4(z1ij3 − u2) −u3√ u24 (z1ij3 − u2) )} )} z16 = {z1ij6 = fc(z1ij3)} ≡ F (16, Z16, U16), z17 = {z1ij7 = Vij = z1ij6(z1ij4 − z1ij2z1ij5)} ≡ F (17, Z17, U17),
≡ F (13, Z13, U13), ≡ F (14, Z14, U14), ≡ F (15, Z15, U15), I ∑ i=1
The energy E of the atoms’ system with the help of new variables may be rewritten as follows: E(z(u)) = I ∑
z1ij7(u). Variables z1, z2, ..., z17 (the phase variables) are determined by the specified above multistep component zl depends on a number of other components (zlij or zijk). l 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 Let us introduce also the following designations: ze1, ze2, ..., ze17 and ez1, ez2, ..., eze17, where e e zn = zijk : zijk = e en en ∂znijk(reij , reik, rejk) ∂ρ ρ=0 = ∂znijk ( (1+ aρ )rij , (1 + aρ )rik, (1 + aρ )rjk
) ∂ρ zn = zenij : zenij = e ∂znij (reij ) ∂ρ ∂znij ( (1+ aρ )rij
) ∂ρ
, ρ=0 ijk : ezen = ∂2zijk(reij , reik, rejk) n ∂ρ2 ∂2znijk ( (1+ aρ )rij , (1 + aρ )rik, (1 + aρ )rjk
) ∂ρ2 ezn = e ij ij een : ezen = z
ρ=0 ∂2znij (reij )
∂ρ2 The above values are calculated by the formulas: ρ=0 = ρ=0 = = ijk z ee8
ijkzijk + 2ze4ijkzijk + zijkzijk = zijk(eze4 5 7 e7 4 ee7 );
ij eze9 = ijk z ee7
ijk z ee4 = ∂2fc ( (1+ aρ )rik
) ∂ρ2 = a32 z6ijkzijk(u7)3(3z6ijk(u7)3 + 2); 7
I ∑ ze1ij7 = ze1ij6z1ij4 − ze1ij6z1ij2z1ij5 + ze1ij4z1ij6 − ze1ij2z1ij6z1ij5 − ze1ij5z1ij6z1ij2; ∂2znij ( (1+ aρ )rij
) ∂ρ2 k ̸= i, j). zijk = 0; e3 e9 zij = zijk = 3z6ijkzijk(u7)3/a; e7 7
I ∑ ;
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: Thus, B(E) is calculated by the formula B(E) = ∑ I ij ∑ eze17(u), where the variables z1,z2,...,z17, ze1,ze2,...,ze17, ez1,ez2,...,eze17 are determined by mentioned above multistep algorithm.
e e
According to the formulas of FAD-technique, the conjugate variables corresponding to the phase variables z1,z2,...,z17, ze1,ze2,...,ze17, ez1,ez2,...,eze17 are determined by the formulas:
e e ijkzijk + 2ze4ijkzijk + z4ijkezei7jk pi5jk = (eze4 7 e7
ijk )epe8 + (ze4ijkz7ijk + z4ijkze7ijk)pei8jk + z4ijkz7ijkpi8jk; pi7jk = ezei4jkz5ijkepei8jk + ze4ijkz5ijkpei8jk + z4ijkz5ijkpi8jk + a32z6ijk(u7)3(3z6ijk(u7)3 + 2)epei7jk + a3z6ijk(u7)3pei7jk; pi8jk = pi9j; pi9j = u6pi1j0; pi1j0 = u5(u5 − 1)(u5 − 2)(z1ij0)u5−3(ze1ij0)2 + u5(u5 − 1)(z1ij0)u5−2ezei1j0]epei1j1+ [ +u5(u5 − 1)(z1ij0)u5−2zij e10pei1j1 + u5(z1ij0)u5−1pi1j1; pi1j1 = − (1 + 4u5)(1 + 2u5)(1 + z1ij1)−2u15−3(z1ij1)2epei1j2 + (14+(u52)u25)(1 + z1ij1)−2u15−2×
8(u5)3 ×ezei1j1epei1j2 + (14+(u52)u25)(1 + z1ij1)−2u15−2ze1ij1pei1j2 − 2u15 (1 + z1ij1)−2u15−1pi1j2; pi1j2 = (−ezei1j6z1ij5 − 2ze1ij5ze1ij6 − z1ij6keei1j5)epe17;
ij z
ij 2(u3)2(z1ij3)2 ij pi1j5 = (−ezei1j6z1ij2 − 2ze1ij2ze1ij6 − z1ij6ezei1j2)epe17 + epe15 − √2/u4 a2u4 u3z1ij3pei1j5; a pi1j6 = (eze14 − ezei1j2z1ij5 − 2ze1ij2ze1ij5 − z1ij2ezei1j5)epe17;
ij ij pijk = 2ze4ijkzijkpijk + z4ijkzijkpei8jk; e7 5 ee8 5 pei1j0 = 2u5(u5 − 1)(z1ij0)u5−2epei1j1 + u5(z1ij0)u5−1zij e10pei1j1;
pi1j7 = 0; pijk = pei9j; e8
pei9j = u6pei1j0; pi1j4 = eze16epe17 + ij ij 2(u3)2u4(z1ij3)2 ij
epe14 − u3√2u4za1ij3pei1j4; a2 (1 + z1ij1)−1/(2u5)−2zij ij e11epe12 −
ij ij pe12 = −2(ze1ij6z1ij5 + zij z1ij6)epe17; e15 ij ij pe15 = −2(ze1ij6z1ij2 + zij z1ij6)epe17; e12 I ij ∑ eze17(u) with respect to independent variables um, (m = 1, 10) are determined by the relations: u1 ( z1ij4 pi1j4 + z1ij5 pi1j5 )
; u1 (z1ij4 (−√2u4(z1ij3 − u2))pi1j4 + z1ij5 (−√2/u4(z1ij3 − u2))pi1j5)+ i=1 j = 1 ∂u6 i=1 j = 1 ((u9)2 + (u10 − z3ijk)2)2 pi5jk ; j ̸= i ∂Ω ∂Ω j ̸= i
I ∑
I ∑ k = 1 k ̸= i, j
I ∑ (( 3z6ijkzijk (
7
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).