Optimal interatomic potentials using modified method of least squares: Optimal form of interatomic potentials

1 Introduction

Phenomenological potential functions [1] are well known to describe the interactions (or forces) between neighbour and/or adjacent atoms [2]. It is therefore necessary to select a potential function that is most appropriate towards the intended aim of the particular experiment(s) or computational software. Although, in quantum mechanics, the Lennard–Jones potential seems to be the most preferred potential. We will be confirming this hypothesis by classifying selected potentials, and after performing a one-dimensional search, the most appropriate potential will be proposed. Interatomic potentials theory is a field of study teeming with possibilities due to its modern applications in quantum mechanics [3], nanotechnology, and nanoengineering [1]. In this modern times, a lot of scientific effort has been channelled towards proposing and modifying interatomic potentials for greater computational efficiency. Between interatomic potentials, we often meet with phenomenological potentials having the generalized mathematical representation

where U ˜ ( ρ ) = U ˜ is the interatomic potential depending on ρ = ρ ( r ) -function of the interatomic distance, r . A , − B > 0 are scalars as well as α > β > 0 , which are scalar exponents. The interatomic distance, ρ = ρ ( r ) , largely determines which potential will result from the generalized representation in Eq. (1). For example, at ρ ( r ) = r , the classical Morse potential

will be obtained. In potential (2), A = ε e 2 β r m , B = − 2 ε e β r m , ε is the depth of the energy minimum, r m is the corresponding minimum value of equilibrium distance between two atoms, and β is the steepness of the negative exponent term. At the same ρ ( r ) = r , the generalized Morse potential, which is also mentioned in the study by Teik-Cheng [4] as the Biswas-Hamann potential [5,6],

will be obtained. Finally, the classical Rydberg potential [5,6]

is also obtainable. The original form of this potential was proposed in the study by Rydberg [7]. Potentials presented in Eqs. (2)–(4) are referred to as the second-order potentials because they can be treated as solutions of the second-order homogeneous ordinary differential equation (ODE)

where ρ = ρ ( r ) , e.g., ρ = r , ρ = ln r , ρ = r 2 . Values of a and b are dependent on Eqs. (2), (3), or (4).

Eigenvalues of ODE (5) are obtained from the characteristic equation, λ 2 + a λ + b = 0 and are

where a = α + β > 0 , b = α ⋅ β > 0 . If a 2 − 4 b > 0 , and then we obtain the generalized Morse potential (3), with

In a particular case where α = 2 β , we will obtain the classical Morse potential, (2). This situation is realized at 2 a 2 = 9 b in ODE (5). Another particular situation in ODE (5), arises when a 2 = 4 b . In this case, the interatomic potential has the classical Rydberg form (4).

In the previously analyzed cases, we considered only instances where a 2 − 4 b ≥ 0 . It is worthwhile to consider a situation with a 2 − 4 b < 0 . In this case, the interatomic potential is described in a real form as follows:

where δ = − a 2 and ω = − 4 b − a 2 2 . Potential (7) was referred to as the modified generalized Morse potential in [8]. Other different forms of the function, ρ = ρ ( r ) , will also give some other classical interatomic potentials. For example, at ρ ( r ) = ln r , we will obtain

which is the generalized Lennard–Jones potential [4,9], originally proposed by Jones [9]. In particular case, at α = 12 and β = 6 , we have the well-known classical 12 − 6 Lennard–Jones potential. Considering another example, at ρ ( r ) = r 2 , we have the two-body portion of the Kaxiras–Pandey potential function [10,11]

Elaborate discussions on the theory of interatomic potentials are extensively detailed in previous studies [1,12–15]. In tandem with formula (8), let us propose analogues of the Lennard–Jones-Rydberg potential

and along with formula (7), let us also propose the analogue of the modified Morse–Lennard–Jones potential

Analogously, we propose further generalizations of the Kaxiras–Pandey potential (9) in the following forms as the Kaxiras–Pandey–Rydberg potential

and the modified Morse–Kaxiras–Pandey potential

All these potentials were previously considered as functions with real values of parameters [8,16– 25,27]. The first generalization of the previously mentioned potentials is being presented in this article, where potentials (3), (8), and (9) were considered for both real and complex conjugate pairs α , β , A , and B . This point of view is a substantial breakthrough because it increases the domain in which the values of parameters are being searched. This generalization directly follows from the proposed approach in the frames of which the functional forms of the potentials are considered as possible solutions of the second-order ODE (5) with constant (real) values a , b . The research question we seek to answer in this article is: What is the most appropriate form of ρ = ρ ( r ) − function of the interatomic distance, r ? Such that this optimal potential will enable nano-scientists to adequately observe interactions at the atomic level, and make sound judgments or conclusions from their studies. To answer this question, we will formulate a generalized form for each of the three forms of ρ = ρ ( r ) previously discussed and use the method of least squares to estimate the unknown parameter, γ , of the generalized potentials.

2 Formulation of the optimization problem

Synthesizing the previous mathematical representations, let us now formulate the problem of optimization of the phenomenological interatomic potential representation in either Lennard–Jones forms (8), (10), or (11). We can rather consider the forms

where we considered ρ = ρ ( r ) as ρ = r γ and γ is the parameter subjected to optimization. Alternatively, we can also consider the synthesized form

or

Let us remark that (14)–(16) can be treated as solutions of the ODE

or after change of variables ρ = r γ by ODE (5). In this case, we use form (14) if α and β are real and different; form (15) if α ≈ β > 0 , and form (16) if eigenvalues in ODE (5) are complex conjugates, λ 1 , 2 = δ ± i ω ( i 2 = − 1 ) , where δ = Re ( λ 1 , 2 ) and ω = ∣ Im ( λ 1 , 2 ) ∣ .

We formulate the problem of optimization of the parameter representation of the interatomic potential as follows:

where U i exp = U exp ( r i ) are the experimental values of the potential at interatomic distance, r = r i . To solve this problem, we select interval of values of parameter γ ∈ [ γ min , γ max ] , divide it by M subintervals, calculate, and then we minimize the function G at every value

for m = 0 , 1 , … , M . Further, we recalculate r i ⟶ ρ i = r i γ m and integrate ODE (5) from ρ min to ρ

where U ˜ min = U ˜ ( ρ = ρ min ) and U ˜ ′ ( ρ ) = d U ˜ ( ρ ) d ρ . Integrated ODE (19) can also be expressed as follows:

where I ˜ 1 ( ρ ) = ∫ ρ min ρ U ˜ ( σ ) d σ , and c = − ( U ˜ ′ min + a U ˜ min ) are the new unknown parameters. After the second integration, we obtain

where Δ U ˜ ( ρ ) = U ˜ ( ρ ) − U ˜ min , I ˜ 2 ( ρ ) = ∫ ρ min ρ I ˜ 1 ( σ ) d σ , Δ ρ ˜ = ρ − ρ min . Model (21) is linear with respect to unknown parameters a , b , and c , which are calculated using the least squares method through the goal function minimization:

where I 1 , i = I ˜ 1 ( ρ i ) , I 2 , i = I ˜ 2 ( ρ i ) , Δ ρ i = ρ i − ρ min , Δ U i = Δ U ˜ ( ρ i ) . The solution of this problem is

where M 1 = [ I 1 I 2 Δ ρ ] is ( N × 3 ) matrix, I 1 = ( I 1 , i ) , I 2 = ( I 2 , i ) , Δ ρ = ( Δ ρ i ) and Δ U = ( Δ U i ) are ( N × 1 ) − columns. System (23) was obtained by partially differentiating goal function (22) with respect to a , b , c , setting the derivatives to 0 and writing the resulting equations in the matrix form.

By calculating eigenvalues λ 1 , 2 using formula (6), we make a choice between the representations (14)–(16). Namely, if λ 1 and λ 2 are real and well separated, we use (14); if they are complex conjugates, we use (16). If λ 1 and λ 2 are close to each other, it is reasonable to use representation (15). To find the proper value of α in (15), it is worthwhile to repeat the procedure of the algorithm explained by Kikawa et al. [22], for the classical Rydberg potential identification.

Next step involves definition of parameters A and B in (14)–(16). In this case, we use basic functions f ˜ 1 ( r ) and f ˜ 2 ( r ) , which are equal to

In this case, the general representation of formulas (14)–(16) reduces to:

and the second goal function, which is subjected to minimization is:

The solution of this conventional problem is

where M 2 = [ f 1 f 2 ] is ( N × 2 ) matrix, and f 1 = ( f 1 , i ) , f 2 = ( f 2 , i ) , and U = ( U i ) are ( N × 1 ) -vector-columns.

Repeating this algorithm M times for the selected γ m ( m = 1 , 2 , … , M ) , we find γ ˜ min = min ( m ) ( γ m ) . This value can be used as guess value for further search of the optimal value of γ . It is also possible to round γ ˜ min to the closest integer for faster calculations of values of the interpolated interatomic potential.

As we have explained earlier, γ can be 1 (like in the case of generalized Morse, classical Morse, and classical Rydberg potentials); can be 2 (like in the case of Kaxiras–Pandey potential and it’s analogues); or can also be 3 for some modified potentials. The treatment of all potentials as solutions of the corresponding ODE (5) immediately defines the form of the potential and domain (real or complex) for its parameters. Moreover, ODE (5) is linear with respect to parameters a and b , which guarantees unique solution to the problem of determination of exponential parameters α and β . This is the novelty of the presented research instead of the conventional “brute force” fitting of the original transcendent potentials. It is necessary to mention that this fitting is performed on the last stage of fitting when the guess values of the parameters are determined by means of the subtle method of minimization of two quadratic goal functions, which was previously formulated in [26] and has been extensively used by the authors in previous studies [8,16– 25,27].

3 Numerical simulations

In what follows, we will be performing numerical simulations of goal functions (18), (22), and (26). We will start by first assuming numerical values for the parameter, γ , to be γ = 0.5 or 3; then we obtain estimates for α and β , and use these estimates to calculate the estimates for A and B . Finally, we carry out the one-dimensional search in order to determine the most appropriate value of γ . The process will be carried out for experimental data sets of metals such as gold [28], copper [29], aluminium [30], titanium [30], and copper–silver alloy [29]. The experimental data sets of the metals were downloaded from the online repository, Interatomic Potentials Repository website^[1]. The form of the interatomic potentials considered has already been presented in Eqs. (14)–(16). It should be noted that the closer the values of the goal function to 0, the more accurate the parameter estimates are. The closer the goal function values are to 1, the more inaccurate the parameter estimates are. For gold atoms [28], and γ = 0.5 , the following values were obtained

hence, α = 12.36613 and β = 5.79498 . The form (14) is preferred since both values are real, and distinct. Numerical simulations then yield A = 2.00117 × 1 0 8 , and B = − 4.94973 × 1 0 3 . The one-dimensional search for parameter, γ (steepest point of the presented curve), yielded γ = 0.38847 (estimated from Figure 1)

Figure 1

One-dimensional search for parameter, γ , using gold atoms.

When we considered the form, ρ ( r ) = ln r , for the Lennard–Jones potential and its analogues, the following values were obtained

The estimated value for parameter γ was then used to construct and reconstruct potential energy curves (PECs) using downloaded experimental data sets of gold atoms and parameter values fitted to potential (14).

The error curves between the two PECs in Figure 2 shows a dimension of magnitude 4 × 1 0 − 3 . In Figure 2 (and other subsequent figures), U n represents the fitted experimental data, U U U n represents the fitted estimated parameters, while ρ n represents the interatomic distance. For copper atoms [29], and γ = 3 , the following values were obtained

therefore, α = 0.35283 and β = 0.04893 . The form (14) is preferred since α and β are real, and distinct. Hence, A = 18.66203 and B = − 0.50888 . The one-dimensional search for parameter, γ (steepest point of the presented curve), yielded γ = 0.01659

Figure 2

Comparison of PECs for original data sets and estimated parameters.

From Figure 3, we can infer that the minimum of the optimized goal function has two values, 0 or 3 . This means we can use ρ ( r ) = r γ , where γ = 3 , or we could also use ρ ( r ) = ln r for the optimal potential. In principle, it will be more efficient to use ρ ( r ) = ln r . When we considered the form, ρ ( r ) = ln r , for the Lennard–Jones potential and its analogues, the following values were obtained

The estimated values using ρ ( r ) = ln r was then used to construct and reconstruct, PECs using experimental data sets of copper atoms and parameter values.

Figure 3

One-dimensional search for parameter, γ , using copper atoms.

The error curves between the two PECs in Figure 4 shows a dimension of magnitude 1 × 1 0 − 3 . For aluminium atoms [30], and γ = 3 , the following values were obtained

therefore, α = 0.19517 , β = 0.04052 . The form (14) is preferred since α and β are real, and distinct. Hence, A = 6.06236 and B = 0.92577 . The one-dimensional search for parameter, γ , yielded γ = 1.43507 .

Figure 4

Comparison of PECs for original data sets and estimated parameters.

In the case of Figure 5, the minimum lies between 1 and 2. Hence, γ = 1 or 2 can be used. This means for aluminium atoms, the optimal potential is either the generalized Morse or the Kaxiras–Pandey potential. When we considered the form, ρ ( r ) = ln r , for the Lennard–Jones potential and its analogues, the following values were obtained

The estimated value for parameter γ was then used to construct and reconstruct PECs using experimental data sets of aluminium atoms and parameter values.

Figure 5

One-dimensional search for parameter, γ , using aluminium atoms.

The error curves between the two PECs in Figure 6 shows a dimension of magnitude 2 × 1 0 − 2 . For titanium atoms [30], and γ = 3 , the following values were obtained

therefore, α = 0.25529 and β = 0.03544 . The form (14) is preferred since α and β are real, and distinct. Hence, A = 22.14025 and B = − 2.47284 . The one-dimensional search for parameter, γ , yielded γ = 0.16645 .

Figure 6

Comparison of PECs for original data sets and estimated parameters.

Figure 7 shows the optimized goal function having a minimum of γ ≈ 0 , which means the optimal potential is obtained when ρ ( r ) = ln r . When we considered the form, ρ ( r ) = ln r , for the Lennard–Jones potential and its analogues, the following values were obtained

The estimated value for parameter γ was then used to construct and reconstruct PECs using experimental data sets of titanium atoms and parameter values fitted to potential.

Figure 7

One-dimensional search for parameter, γ , using titanium atoms.

The error curves between the two PECs in Figure 8 shows a dimension of magnitude 2 × 1 0 − 2 . Finally, when we analyze the silver–copper alloy [29], with γ = 3 , the following values were obtained

therefore, α = 0.31588 , and β = 0.02722 . The form (14) is preferred since α and β are real, and distinct. Hence, A = 24.46179 and B = − 0.28071 . The one-dimensional search for parameter, γ , yielded γ = 2.98832 (estimated from Figure 9).

Figure 8

Comparison of PECs for original data sets and estimated parameters.

Figure 9

One-dimensional search for parameter, γ , using silver–copper alloy.

Without a doubt, the optimal potential for silver–copper alloy is realized when γ = 3 . This potential will be an analogues or a modification of the generalized Morse potential. When we considered the form, ρ ( r ) = ln r , for the Lennard–Jones potential and its analogues, the following values were obtained

The estimated value for parameter γ was then used to construct and reconstruct PECs using experimental data sets and parameter values.

The error curves between the two PECs in Figure 10 shows a dimension of magnitude 2 × 1 0 − 3 . As was graphically illustrated in Figures 2, 4, 6, 8, and 10, the obtained parameter estimates for α , β , A , and B used to plot PECs, all showed reasonable agreement with the experimental data sets. The obtained PECs were graphically indistinguishable through a large range of the interatomic distance, r . Table 1 illustrates the values of the goal function for both forms of ρ ( r ) considered and for each metal atoms used for experimental values. It should be recalled that the closer the numerical values are to 0, the more accurate the parameter estimates calculated.

Figure 10

Comparison of PECs for original data sets and estimated parameters.

The goal function values were obtained by minimizing the goal functions (26) and (18) through built in functions in Mathcad^®.

Table 2 concisely summarizes the results obtained from numerical simulations. The preferred choice of potential used for numerical simulation was based on the agreement of the reconstructed PECs with experimental data sets.

4 Discussion and conclusion

In this article, interatomic potentials that can be treated as solutions of some second-order ODE were classified and identified. A generalization of three forms of potentials were presented, and the most appropriate form of a generalized potential will be based on the estimated value of parameter, γ . The goal function values were also stated in Table 1. Reconstructed PECs showed good agreement for a relatively large vicinity of the potential minimum, between the estimates and experimental data sets. It should be noted that there are two ways of optimization: the first one relates to the selection of the “ γ ” exponent. (In the available literature, the authors used γ = 1 for the generalized Morse potential [4] and γ = 2 for the Kaxiras–Pandey potential [10,11]. The question is why not γ = 3 2 , for example, or some other value?.) The second optimization relates to either real or complex conjugate exponents and the corresponding factors (in the direct brute force optimisation, it is impossible to change the field of real valued parameters to the field of complex conjugate ones, it is possible only in the frames of the authors approach using the multiple (double) goal functions method).

In general, we can infer that the form with Lennard–Jones potential has the lower goal function values and hence is the optimal interatomic potential (most preferable potential), for many cases. The one-dimensional search for the most appropriate value of the parameter, γ , fails when copper, titanium, and aluminium were analyzed. However, this defect does not affect the reconstructed PECs. Having performed numerical simulations using different values of γ , we can also infer that when γ = 3 , the optimal potential is obtained and is more preferable as compared to other cases when γ = 1 or 2.