Molecular dynamics simulations for mechanical properties of borophene: parameterization of valence force field model and Stillinger-Weber potential

While most existing theoretical studies on the borophene are based on first-principles calculations, the present work presents molecular dynamics simulations for the lattice dynamical and mechanical properties in borophene. The obtained mechanical quantities are in good agreement with previous first-principles calculations. The key ingredients for these molecular dynamics simulations are the two efficient empirical potentials developed in the present work for the interaction of borophene with low-energy triangular structure. The first one is the valence force field model, which is developed with the assistance of the phonon dispersion of borophene. The valence force field model is a linear potential, so it is rather efficient for the calculation of linear quantities in borophene. The second one is the Stillinger-Weber potential, whose parameters are derived based on the valence force field model. The Stillinger-Weber potential is applicable in molecular dynamics simulations of nonlinear physical or mechanical quantities in borophene.

From the above literature survey, we find that most existing theoretical works are based on the first-principles calculations. These calculations are accurate, but are limited to small (or bulk) structures due to high computation requirements. As an alternative approach, the molecular dynamics simulation can be utilized to investigate very large systems typically containing more than 10 4 atoms. The key ingredient in the molecular dynamics simulation is the interatomic interaction. The only one interaction potential available for the borophene is the ReaxFF force field model 46 . The present work aims to develop efficient linear and nonlinear empirical potentials for the borophene, which can assist further theoretical investigations for borophene of large size.
In this paper, we provide the valence force field (VFF) model and the Stillinger-Weber (SW) potential for the description of the interaction in borophene with low-energy triangular structure. The VFF model is a linear potential, which can be used to calculate elastic quantities such as the phonon dispersion. The SW potential is a nonlinear potential, which is derived based on the VFF model. The SW potential can be applied in the molecular dynamics simulation of nonlinear physical or mechanical properties for borophene. We demonstrate the usage of the SW potential with the publicly available LAMMPS package. Figure 1 shows the structure of borophene, in which structural parameters are from the ab initio calculation 28 . Figure 1(a) shows the top view of the borophene with the rectangular unit cell in the xy plane. The two lattice bases are a 1 = 2.866 Å and a 2 = 1.614 Å. The first Brillouin zone is shown by red rectangle on the left. The two basic vectors for the reciprocal lattice are b 1 = 2.192 Å −1 and b 1 = 3.893 Å −1 . There are two inequivalent boron atoms in the unit cell. Boron atoms are categorized into the top chain and the bottom chain. The top chain includes atoms like 1, 4, and 7. The bottom chain includes atoms like 2, 3, 5, and 6. Figure 1(b) discloses the puckered configuration of the borophene. The pucker is The top chain includes atoms like 1, 4, and 7. The bottom chain includes atoms like 2, 3, 5, and 6. The unit cell is shown by blue rectangle. The first Brillouin zone is shown by red rectangle on the left. (b) Perspective view illustrates the puckered configuration, with h as the distance between the top and bottom chains along the outof-plane z-direction. The pucker is perpendicular to the x-axis and is parallel with the y-axis. perpendicular to the x-direction, while it is parallel with the y-direction. The height of the pucker is h = 0.911 Å, which is the distance between the top chain and the bottom chain along the out-of-plane z-direction. Valence force field model. The VFF model is a useful linear model for the description of interatomic interactions in covalent materials, in which interactions are decomposed into some characteristic bond stretching and angle bending components 47 . These interaction components are in close relation with the vibration morphology of phonon modes. It is thus a proper approach to determine the VFF model for a covalent material based on its phonon dispersion. There are two typical terms for the VFF model, i.e., the bond stretching V r and the angle bending V θ , 2 where Δ r and Δ θ are the small variations for the bond length and the angle. The two force constant parameters are K r and K θ . Table 1 shows five VFF terms for the borophene, four of which are the bond stretching interactions shown by Eq. (1) while the other two terms are the angle bending interactions shown by Eq. (2). The VFF model for the borophene is determined with the assistance of the phonon dispersion. The phonon dispersion for the borophene has been obtained by previous ab initio calculations. Parameters for the VFF model are determined by fitting the phonon dispersion to the ab initio calculations 28 . GULP 48 is used for the computation of the phonon dispersion and the fitting of parameters in the VFF model. Figure 2 shows that the phonon dispersion calculated from the VFF model is comparable with ab initio results. Good agreements are achieved for the three acoustic branches in the long wavelength limit. Similar phonon dispersion can be found from other ab initio calculations 26,29,38 . There is no imaginary mode from the VFF model, which appears in the ab initio calculations for the flexure mode in the long wave limit 26,28,29,38 . This imaginary mode is due to the long-range interaction in the borophene, which is not considered by the VFF model. Hence, there is no imaginary mode here. Figure 3 shows the vibration morphology of the three optical phonons at the Γ point. Each boron atom chain vibrates as a rigid chain in these phonons. The three phonon modes at the Y point are doubly degenerate, due to the inverse reflection symmetry between the top and bottom boron chains in the borophene. We have presented the vibration morphology for the three phonons at the Y point in Fig. 4 for the bond stretching interactions, and in the unit of eV for the angle bending interaction. The fourth line gives the initial bond length (in unit of Å) for the bond stretching interaction and the initial angle (in unit of degrees) for the angle bending interaction. The angle θ ijk has atom i as the apex. 1200 28 . Phonon branches are doubly degenerate along YX due to the inverse symmetry between the top and bottom boron chains.

Figure 2. Phonon dispersion for borophene along ΓYSXΓ from the VFF model (blue lines) is compared to the data from the ab initio calculation (black dots)
here), only atoms from the bottom chains are involved in the vibration. Actually, the phonon branches are doubly degenerate along the whole YX line, due to the inverse reflection symmetry in the borophene.
An obvious feature in the phonon dispersion of the borophene shown in Fig. 2 is that the highest-frequency branch (around 1300 cm −1 ) along SX has much higher frequency than the other branches. Figure 5 shows the vibration morphology of these phonons at the X point. The highest-frequency phonon is the intra-chain optical vibration; i.e., neighboring boron atoms within the same chain are vibrating in an out-of-phase manner. It means that the intra-chain bonds like r 14 in Fig. 1(a) are much stronger than the inter-chain bonds like r 12 . This anisotropy is reflected in the force constant parameters for the VFF model listed in Table 1. The intra-chain parameter K 14 is larger than the inter-chain parameter K 12 by a factor of four, so the resultant frequency for the intra-chain optical phonons (around 1300 cm −1 ) is about twice of the frequency for the inter-chain optical phonons (around 700 cm −1 ) shown in Fig. 2.
Parameterization of the SW potential. The VFF model presented in the previous section describes the linear component of the interaction within borophene, but it does not provide any information for the nonlinear interaction. The SW potential is one of the most efficient nonlinear potentials. The SW potential includes both linear and nonlinear interactions. One of the present authors (J.W.J.) has shown that the linear component of the SW potential can be derived analytically based on the VFF model 49 . A distinct feature of this analytic approach is that the parameters of the SW potential are governed by a constraint, which guarantees that all bonds and angles are fully relaxed in the initial configuration. The constraint is shown as follows, where d is the equilibrium bond length, and other quantities are potential parameters. The nonlinear component of the SW potential is determined by the parameter B, which can be determined by the third-order nonlinear constant.
The parameters for the two-body SW potential can be derived from the bond stretching terms in the VFF model. The obtained parameters are shown in Table 2. The parameters for the three-body SW potential are obtained from the angle bending terms in the VFF model. The parameters are shown in Table 3. Using this SW potential, we calculate the phonon dispersion for the borophene, which is the same as the phonon dispersion calculated using the VFF model as shown in Fig. 6. It indicates that the SW potential has inherited the linear interaction of the VFF model.
The publicly available package LAMMPS is widely used for MD simulations. The SW potential has a slightly different form in LAMMPS. The parameters for the SW potential can be determined by comparing the potential forms in GULP and LAMMPS. Table 4 shows some representative parameters for the SW potential used by LAMMPS. It should be noted that we have introduced eight atom types for boron atoms in the borophene as shown in Fig. 7. It is because the cut off (4.0 Å) for the two-body SW potential between atoms like 3-5 is larger than the distance between atoms 3-6 (3.289 Å). However, there is no interaction between atoms like 3-6. It is thus necessary to differentiate bonds like 3-5 and 3-6 with atoms 3 and 6 denoted by different atom types, so that the two-body SW potential is not imposed on the bonds like 3-6. There are 512 lines (categorized into six classes) in the SW potential script used by LAMMPS. We have shown in Table 4 one representative term for each of these six classes. These files for LAMMPS can be found in the Supplemental materials, including the SW potential script (borophene.sw), an input file for LAMMPS (in.tension), and a structure generation code (xyz.f90).
The quantity (r ij ) in the first line lists one representative term for the two-body SW potential between atoms i and j. Atom indexes are from Fig. 1(a).  θ θ − (cos cos ) 0 2 . The first line (θ ijk ) presents one representative term for the three-body SW potential. The angle θ ijk has the atom i as the apex. Atom indexes are from Fig. 1(a).  Table 4. SW potential parameters used by LAMMPS 52 . The two-body potential expression is Scientific RepoRts | 7:45516 | DOI: 10.1038/srep45516 The SW potential is not able to describe accurately the complicate chemical bonding processes. However, it contains reasonably accurate nonlinear components, so it is able to provide some qualitative descriptions for the thermal-induced nonlinear phenomena. As an example, Fig. 10 shows the total potential energy in the borophene at varying temperatures. A sudden jump at the critical temperature T C = 550 ± 50 K reveals that the borophene is stable at temperatures bellow the critical temperature T C . However, the structure becomes unstable for temperatures above the critical temperature. Figure 11 shows the destruction process for the borophene at 600 K. Some boron atoms have obviously larger thermal vibration amplitudes at t = 290 ps. These boron atoms will be    We note that the VFF model and the SW potential parameterized in the present work are only applicable to describe the interaction for the triangular phase of the borophene, but can not be applied to other phases for the borophene.

Discussion
To summarize, we have developed two empirical potentials to describe the interaction between boron atoms within borophene. The first one is the linear VFF model, which is determined based on the phonon dispersion of the borophene. The VFF model is suitable for the calculation of linear quantities in borophene. Starting from the VFF model, we derive the second empirical potential, the SW potential, to describe the interaction for borophene. The SW potential includes both linear and nonlinear couplings, so it is applicable for the simulation of some nonlinear quantities. We also provide some necessary files for the application of the SW potential by using LAMMPS.   Method One minor improvement for the implementation of three-body SW potential in LAMMPS. We discuss a practical improvement for the three-body SW potential implemented in LAMMPS. Let's take angles θ 143 and θ 137 in Fig. 1 as an explicit example to demonstrate this improvement.
It is obvious that these two angles should be quite different. As shown in Table 3, there is a three-body SW potential for θ 143 , but there is no interaction for θ 137 . The values of these two angles are quite different in the equilibrium configuration; i.e., θ = .  64 581 143 0 and θ = .  115 419 137 0 . However, their arm lengths are the same; i.e., r 13 = r 13 and r 14 = r 17 . In the three-body SW potential implemented in LAMMPS, angles are distinguished by their two arm lengths, so angles θ 143 and θ 137 become indistinguishable. Consequently, the three-body SW potential (intent solely) for θ 143 will also be applied to θ 137 , which causes a technical issue.
In general, it will become important to distinguish two angles constructed by the same type of atoms, when there are more than two different types angles around each atom. Similar situation also occurs in the simulation of MoS 2 using the three-body SW potential implemented in LAMMPS. We have proposed to distinguish each angle by its opposite arm length 49,51 . That is, the bonds r 34 and r 37 are quite different in length, so they can be used to distinguish angles θ 143 and θ 137 . However, this bond length based criterion is an ad hoc technique, as the bond length is a structural dependent parameter. We note that this criterion is also adopted by GULP 48 .
Here, we provide a more universal criterion to distinguish two angles. It is based on the input initial value θ 0 for the angle of the three-body SW potential. For example, in Table 3, the input initial value for the three-body SW potential regarding θ 143 is θ = .  64 581 143 0 . For a given angle θ, if the following condition is satisfied, then the three-body SW potential will be applied to this angle, It can be checked that θ 137 does not obey the criterion in Eq. (4) and it is thus excluded from the three-body SW potential. The angle θ 143 obeys the criterion in Eq. (4) under moderate deformations. Hence, the three-body SW potential is applied to the angle θ 143 , while this three-body SW potential is not applied to the other angle θ 137 .
The criterion in Eq. (4) can be implemented into LAMMPS by modifying its pair_sw.cpp source file in the following two steps. A continuous cut-off function in the range of [0.25, 0.35] has been introduced to avoid possible boundary effects in molecular dynamics simulations. (1). First, find the following line in the pair_sw.cpp source file, Then recompile the LAMMPS package. The recompiled LAMMPS executable file can be used to simulate borophene with the SW potential parameterized in the present work. It can also be used to simulate MoS 2 correctly using its SW potential. We note that this criterion does not affect other normal simulations using LAMMPS, where it is not necessary to distinguish angles. We expect this criterion to be implemented in LAMMPS in its future versions.