Multi-fidelity Gaussian process based empirical potential development for Si:H nanowires

Moonseop Kim, Huyi Yin, Gung Lin,c,*

a School of Mechanical Engineering, Purdue University, West Lafayette, IN 47906-2045, USA

b School of Computer and Information Engineering, Xiamen University of Technology, Xiamen 361024, China

c Department of Mathematics, Purdue University, West Lafayette, IN 47906-2045, USA

Keywords:Multi-fidelity Gaussian process regression Inter-atomic potential and forces Elasticity

ABSTRACT In material modeling, the calculation speed using the empirical potentials is fast compared to the first principle calculations, but the results are not as accurate as of the first principle calculations.First principle calculations are accurate but slow and very expensive to calculate. In this work, first,the H-H binding energy and H2-H2 interaction energy are calculated using the first principle calculations which can be applied to the Tersoff empirical potential. Second, the H-H parameters are estimated. After fitting H-H parameters, the mechanical properties are obtained. Finally, to integrate both the low-fidelity empirical potential data and the data from the high-fidelity firstprinciple calculations, the multi-fidelity Gaussian process regression is employed to predict the HH binding energy and the H2-H2 interaction energy. Numerical results demonstrate the accuracy of the developed empirical potentials.

In the last three decades, empirical potentials have been advanced. With the advance of supercomputers, these potentials are anticipated to be widely used for the next three decades [1].Atomistic calculations by empirical potentials can be utilized in understanding the structural aspects of Si or Si-H systems that are found in many important areas such as the surface of nanopatterning Si [2, 3], nano-electro-mechanical systems (NEMS)[4], superconductivity of silane [5], optical modulators [6], and applications of α-Si:H materials [7]. In the past, empirical potentials for Si [8–11] and for Si-H [12–14] have been developed. But the bulk elastic properties of Si cannot be resolved using such empirical potentials. In Ref. [12], it has been shown that at the hydrogen-induced reconstruction of the silicon surface, the distance between hydrogen and hydrogen is 1.64 ? (1 ? = 1×10-10m) and bond angle H-Si-H is 106° using existing empirical potential. However, when H-H distance and the bond angle are compared with the results from the first-principle calculations,H-H distance and the bond angle are 2.1638 ? and 104.805° respectively. In this situation, the bond angle is distinguished from 1.195°, which means that the difference of the bond angles can be ignored, however, the biggest issue is that H-H distance is distinguished from 0.5238 ?. Hence, if the existing empirical potential is used for Si nanowires, the computation speed is fast, but the results obtained from the existing empirical potentials are not accurate compared to the results from the first-principle calculations. It is critical to fix such errors. In this paper, we propose two novel techniques to construct the empirical potentials with an emphasis on parameter fitting and multi-fidelity modeling, in which the relationship between material properties and potential parameters is explained. The input database has been obtained from the density functional theory (DFT) calculations with the Vienna ab initio simulation package (VASP) [15]. This paper is constituted as follows. First, the structure of silicon nanowires passivated hydrogen is introduced with some specific shapes. Second, governing equations of Tersoff empirical potentials are presented to explain which parameters can be obtained from H-H binding energy and H2-H2interaction energy.Third, we give a brief explanation of the multi-fidelity Gaussian process regression for prediction of the results of H-H binding energy and H2-H2interaction energy. Fourth, we represent three optimization methods for H-H parameter fitting. The rootmean-square-error obtained by the Nelder–Mead simplex method is compared with the results from the other two optimization methods. Lastly, we evaluate the mechanical properties (Young's modulus and equilibrium elongation) using the estimated parameters obtained from the H-H parameter fitting.

In this study, Si nanowires passivated hydrogen model is chosen. If Si nanowires have dangling bonds, it will oxidize in the air circumstance. By passivating hydrogen to the surface of Si nanowires, it can be stabilized from the oxidization. Figure 1 has expressed a cross-section of silicon nanowires passivated hydrogen [16]. Green dots and blue dots represent silicon atoms and hydrogen atoms, respectively. Cross-section of Si nanowires represents by the Wulff structure selected by minimizing the surface energy. It can be stabilized by passivating hydrogen to the surface of Si nanowires. In the mechanical property, Young's modulus is calculated after the H-H parameter fitting by increasing the size of cross-section compared to the results of Young's modulus using the existing empirical potentials and the firstprinciples calculations.

The atomistic computer simulations based on the empirical potential is fast for calculation. In this system, the number of atoms is not limited compared to the first principles, however,the accuracy of calculation is not adequate, therefore, reliability of the empirical potential presented so far is needed to verify.Various empirical potentials depend on the material, for instance, Ni and Ti are calculated through embedded atom method (EAM) [17], and Si is calculated through Tersoff empirical potential [9] and Stillinger-Weber empirical potential [8]. In this study, silicon nanowires passivated hydrogen model, and Tersoff empirical potential are used to verify the accuracy of existing Tersoff empirical potentials. Tersoff empirical potential is based on the concept of bond order, the force of bonds between atoms is not consistent and depends on the local environment.The total energy function is given as [12]

Fig. 1. Cross-section of silicon nanowires passivated hydrogen<0 01> .

where V is total energy function, fR(r) and fA(r) are repulsive energy and attractive energy respectively. These functions are defined as function of distance between i and j atoms, r is interatomic distance and bijis bond order. In this study, A, B, λ1and λ2are decided as H-H fitting parameters.

where ζijis the function of effective coordination number, H(N)is the function of bond number, cos θijkis the bond angle, rijand rikare the distance between i and j atoms and between i and k respectively andandare equilibrium distance between i and j atoms and between i and k respectively. Lastly, in this study, α, β, η, δ , and c are determined from H2- H2parameter fitting.

where fc(n) is a cutoff function determining whether there is coherence or not between the atom and its neighbor atom. r is the interatomic distance, R and D, the influence is 1. If the interatomic distance is larger than R–D, the influence is 0.Finally, if the interatomic distance is between R–D and R+D, it is influenced by Eq. (5).

Here, we provide the steps for multi-fidelity modeling with Gaussian processes (GP). The steps on multi-fidelity are given as

where u1(x) and u2(x) are independent. In the Gaussian process regression, it is assumed that the mean of GP is zero and k(x,x′;θ ) is the covariance matrix between all possible pairs ( x,x′) in the set of vectors of hyper-parameters θ. As shown in Ref.[18], the basic idea is that we begin with two independent GP u1(x) and u2(x); then we define the low-fidelity and the highfidelity models [19–24]

This demonstrates the "relationship" between the low- and highfidelity models since both include the G P u1(x). In particular,setting k1= cov[u1, u1] and k2= cov[u2, u2] we have:

cov[u1, u2] = 0 and cov[u2, u1] = 0 by independence and to sum up KLL, KLH, KHH:

This gives us a complete model that incorporates both the lowand high-fidelity. In particular, we model the column vector[ fL(x); fH(x)] using a zero-mean prior and the covariance matrix defined block-wise by [ KLL, KLH; KHL, KHH]. Since the mean and covariance are known, the whole Gaussian process model is specified, and the training can be performed using the standard procedure.

To represent the uncertainty (or noise) in the observation data,the covariance of the noise for both the low- and high-fidelity data is added on the diagonal of the covariance matrix in Eq.(16). The level of the noise in the observation data will affect the prediction accuracy as shown in Figs. 2 and 3.

where K is the covariance matrix and the negative log marginal likelihood (NLML) is used as the "cost function" which should be minimized to get the best-fit model by using hyper-parameters θ. In prediction, if we consider a Gaussian likelihood and the posterior distribution is easy to apply and can be used to involve predictive deduction for a new output fH, given a new input x?as

In numerical results, to reduce the computational cost for expensive calculations (DFT), multi-fidelity Gaussian process regression for prediction [19–24] is used. H-H binding energy and H2-H2interaction energy obtained from the empirical potential are applied to the low-fidelity model. Results of H-H binding energy and H2-H2interaction energy obtained from DFT are implied to the high-fidelity model with a limited number of samples due to high computational cost. In Figs. 2 and 3, the number of high-fidelity samples ( NH= 4) is fixed and we compare the standard deviation of the high-fidelity prediction by increasing the number of low-fidelity samples. As we can see,the standard deviation of the high-fidelity is decreased when the number of low-fidelity samples is increased.

Fig. 2. Multi-fidelity prediction results of H-H binding energy with high-fidelity samples ( NH = 4) and three different numbers of low-fidelity samples ( NL = 3, 5, 7). If we increase N L, the accuracy is increased (the standard deviation of the high-fidelity decreases).

Fig. 3. Multi-fidelity prediction results of H 2- H2 interaction energy with high-fidelity samples ( NH = 4) and three different numbers of low-fidelity samples ( NL = 3, 5, 7). If we increase N L, the accuracy is increased (the standard deviation of the high-fidelity decreases).

The structure of Si nanowires is passivated by hydrogen to prevent oxidation. We divide the Si nanowires into three parts for parameter fitting. First, Si-Si parameter fitting is needed for the internal structure of Si nanowires, Second, H-H parameter fitting is required for the surface computation. Lastly, Si-H parameters for the interface calculation are needed. In this study, the H-H parameters fitting is focused on using H-H binding energy and H2-H2interaction energy for surface computation. In the surface of Si nanowires, there are two forces for hydrogen relations that are attractive and repulsive forces. This H-H parameter fitting is complicated so we have to design systematically.When we fit H-H parameters using H-H binding energy, we only use Eqs. (1)–(3) except for Eqs. (4) and (5). This is because H-H binding energy is calculated by two hydrogen atoms which can be neglected ζijterm that needs more than three atoms,however, H2-H2parameter fitting by using interaction energy is applied to more than three hydrogen atoms so it should be applied to ζijterm. In Table. 1, three methods of optimization were used for H-H parameters fitting. Finally, we investigate the error of each method between DFT results and the fitting line through the root-mean-square-error in Table 1. The Nelder–Mead simplex method is the best compared to the other two optimization methods. The parameters using the Nelder–Mead simplex method were chosen, which are listed in Table 1.

In the study, we compute the binding energy using DFT for reference results that are adjusted to the results of the empiricalpotentials objective function. Parameters can be obtained after fitting between hydrogen and hydrogen. In Fig. 4, black dots represent DFT results for H-H binding energy and the red line is the fitting line. A, B, λ1, λ2, and R(e)values are listed on Table 1.

Table 1 To fit H-H parameters for molecular hydrogen, three optimization methods are compared which are Nelder–Mead simplex method (N-M) [25], Broyden–Fletcher–Goldfarb–Shanno quasi–Newton method (BFGS) [26–29], trust-region method (T-R)[30, 31].

We compute the interaction energy using DFT for reference results that are adjusted to the results of the empirical potentials objective function and we obtained parameters after fitting hydrogen inter-molecular. We use the Nelder–Mead simplex methods of optimization for the H-H parameter fitting using the interaction energy. In Fig. 5, the black dots represent the DFT results for H2-H2interaction energy, and the red line is the fitting line. In this H-H fitting, A, B, λ1, λ2, and R(e)are obtained from H-H parameter fitting using the binding energy applied to H2-H2parameter fitting and α, β, η, δ, and c values are presented in Table 1.

Figures 6 and 7 represent the H-H binding energy and H2- H2interaction energy compared to DFT, empirical potential (this work), and existing empirical potential. In Fig. 6 and 7, crossdata represents the H-H binding energy and H2- H2interaction energy using the existing Tersoff empirical potential parameters.It is critical to note that these data's curve shape is not smooth because, in existing classical molecular dynamics (MD), the cutoff range is too short to calculation, so it does not calculate for the long-range interaction. In this study, we fixed cut off range much longer and calculate the long-range of hydrogen molecules. As we can see, DFT and empirical potential (this work)result matches each other after H-H parameters fitting. We have developed a systematic process to construct empirical potential for Si nanowires' passivated hydrogen.

In Figs. 8 and 9, mechanical properties were performed using H-H parameters which are obtained after H-H part parameter fitting. To evaluate the new fit of the H-H part, Young's modulus and equilibrium elongation of Si nanowires are calculated by increasing the wire width of Si nanowires. We compared the numerical results with the DFT results and existing empirical results in Figs. 8 and 9. The size reliance of Young's modulus and equilibrium elongation shows critical improvement compared to the DFT results and the existing potential results. Until now,the surface part of Si nanowires is fitted using our systematic fitting method and shows mechanical properties to prove enhancement. However, the improvement of the irregular mechanical properties can be observed by the H-H parameter fitting of the surface. But the perfect result of matching could not be observed. The reason for this is that not only the surface but also the silicon-hydrogen parameter fitting between the silicon and the surface must be performed. In addition, the hydrogen parameters obtained so far are limited to the calculation of the mechanical properties of nanowires with hydrogen and silicon. Furthermore, it is necessary to derive parameter fitting that can be applied to various types of materials, and potential errors of existing empirical potentials must be corrected for calculations using various materials as well as silicon nanowires.

Fig. 4. H-H parameter fitting using the Nelder-Mead simplex method.

Fig. 5. H2-H2 parameter fitting using the Nelder-Mead simplex method.

Fig. 6. Binding energy after H-H parameter fitting.

Fig. 7. Interaction energy after H-H parameter fitting.

Fig. 8. Young's modulus increasing as wire width of Si nanowires increasing.

Fig. 9. Equilibrium elongation increasing as wire width of Si nanowires increasing.

Acknowledgement

We gratefully acknowledge the support from the National Science Foundation of USA (Grants DMS-1555072 and DMS-1736364).