Hybrid Tersoff+ILP potential

The hybrid Tersoff + ILP potential in GPUMD combines the Tersoff (1988) potential [Tersoff1988] for intralyer interactions and the interlayer potential (ILP) [Ouyang2018] [Ouyang2020] for interlayer interactions to simulate homo- and heterostructures based on graphene and \(h\)-BN layered materials. The Tersoff term of this hybrid potential uses the same implementation as Tersoff (1988) potential.

Potential form

The site potential of ILP can be written as:

\[U_i^{\mathrm{ILP}}= \mathrm{Tap}(r_{ij}) \left[U^{\mathrm{att}}(r_{ij})+U^{\mathrm{rep}}(r_{ij}, \boldsymbol{n}_i, \boldsymbol{n}_j)\right]\]

The function \(\mathrm{Tap}(r_{ij})\) is a cutoff function and takes the following form in the intermediate region:

\[\mathrm{Tap}(r_{ij})=20{\left(\frac{r_{ij}}{R_{\mathrm{cut}}}\right)}^7- 70{\left(\frac{r_{ij}}{R_{\mathrm{cut}}}\right)}^6+84{\left(\frac{r_{ij}}{R_{\mathrm{cut}}}\right)}^5- 35{\left(\frac{r_{ij}}{R_{\mathrm{cut}}}\right)}^4+1\]

The repulsive term and the attractive term take the following forms:

\[ \begin{align}\begin{aligned}U^{\mathrm{att}}(r_{ij})&=-\frac{1}{1+e^{-d_{ij}\left[r_{ij}/(S_{\mathrm{R},ij}\cdot r_{ij}^{\mathrm{eff}})-1\right]}}\frac{C_{6,ij}}{r_{ij}^{6}}.\\U^{\mathrm{rep}}(r_{ij}, \boldsymbol{n}_i, \boldsymbol{n}_j)&=e^{\alpha_{ij}\left(1-\frac{\gamma_{ij}}{\beta_{ij}}\right)} \left\{\epsilon_{ij}+C_{ij}\left[e^{-{(\frac{\rho_{ij}}{\gamma_{ij}})}^2}+e^{-{\left(\frac{\rho_{ji}}{\gamma_{ij}}\right)}^2}\right]\right\}.\end{aligned}\end{align} \]

where \(\boldsymbol n_i\) and \(\boldsymbol n_j\) are normal vectors of atom \(i\) and \(j\), \(\rho_{ij}\) and \(\rho_{ji}\) are the lateral interatomic distance, which can be expressed as:

\[\begin{split}\rho_{ij}^{2}&= r_{ij}^2-{(\boldsymbol r_{ij} \cdot \boldsymbol n_i)}^2\\ \rho_{ji}^{2}&= r_{ij}^2-{(\boldsymbol r_{ij} \cdot \boldsymbol n_j)}^2\end{split}\]

Other variables are all fitting parameters.

The site potential of Tersoff can be written as:

\[U_i^{\mathrm{Tersoff}} = \frac{1}{2} \sum_{j \neq i} f_C(r_{ij}) \left[ f_R(r_{ij}) - b_{ij} f_A(r_{ij}) \right].\]

The function \(f_{C}\) is a cutoff function, which is 1 when \(r_{ij}<R\) and 0 when \(r_{ij}>S\) and takes the following form in the intermediate region:

\[f_{C}(r_{ij}) = \frac{1}{2} \left[ 1 + \cos \left( \pi \frac{r_{ij} - R}{S - R} \right) \right].\]

The repulsive function \(f_{R}\) and the attractive function \(f_{A}\) take the following forms:

\[\begin{split}f_{R}(r) &= A e^{-\lambda r_{ij}} \\ f_{A}(r) &= B e^{-\mu r_{ij}}.\end{split}\]

The bond-order is

\[b_{ij} = \left(1 + \beta^{n} \zeta^{n}_{ij}\right)^{-\frac{1}{2n}},\]

where

\[\begin{split}\zeta_{ij} &= \sum_{k\neq i, j}f_C(r_{ik}) g_{ijk} e^{\alpha(r_{ij} - r_{ik})^{m}} \\ g_{ijk} &= \gamma\left( 1 + \frac{c^2}{d^2} - \frac{c^2}{d^2+(h-\cos\theta_{ijk})^2} \right).\end{split}\]

File format

This hybrid potential requires 2 potential files: ILP potential file and Tersoff potential file. We have adopted the ILP file format that similar but not identical to that used by lammps. To identify the different layers, it’s required to set one group_methods in model.xyz file. Now this hybrid potential could be only used to simulate homo- and heterostructures of graphene and \(h\)-BN.

In run.in file, the potential setting is as:

potential <ilp file> <tersoff file>

where ilp file and tersoff file are the filenames of the ILP potential file and Tersoff potential file. ilp file is similar to other empirical potential files in GPUMD:

tersoff_ilp <number of atom types> <list of elements>
<group_method for layers>
beta alpha delta epsilon C d sR reff C6 S rcut1 rcut2
...

tersoff_ilp is the name of this hybrid potential.
number of atom types is the number of atom types defined in the model.xyz.
list of element is a list of all the elements in the potential.
group_method for layers is the group_method set in model.xyz to identify different layers. For example, monolayer graphene and monolayer \(h\)-BN are both single layer so for the atoms in each layer the group_id of group_method for layers are the same.
The last line(s) is(are) parameters of ILP. rcut1 is used for calculating the normal vectors and rcut2 is the cutoff of ILP, usually 16Å.

More specifically, for graphene, if group_method 0 is used for different layers, the ilp file is required to set as:

tersoff_ilp 1 C
0
beta_CC alpha_CC delta_CC epsilon_CC C_CC d_CC sR_CC reff_CC C6_CC S_CC rcut1_CC rcut2_CC

The tersoff file use the same atomic type list as the ilp file and just contains parameters of Tersoff potential. The potential file reads, specifically for graphene:

A_CCC B_CCC lambda_CCC mu_CCC beta_CCC n_CCC c_CCC d_CCC h_CCC R_CCC S_CCC m_CCC alpha_CCC gamma_CCC

More parameter details are in Tersoff (1988) potential and NEP+ILP potential.