Hybrid SW+ILP potential

The hybrid SW + ILP potential in GPUMD combines the Stillinger-Weber potential (SW) [Stillinger1985] for intralyer interactions and the interlayer potential (ILP) [Ouyang2018] [Ouyang2020] for interlayer interactions to simulate homostructures based on transition metal dichalcogenides layered materials. The SW term of this hybrid potential uses the modification version and more details are in [Jiang2015] and [Jiang2019].

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 modified SW can be written as:

\[\begin{split}U_i^{\mathrm{SW}} =& \sum_i\sum_{j\gt i}\phi_2\left(r_{ij}\right)+ \sum_i\sum_{j\neq i}\sum_{k\gt j}\phi_3 \left(r_{ij}, r_{ik}, \theta_{ijk}\right)\\ \phi_2\left(r_{ij}\right) =& A_{ij}\epsilon_{ij}\left[B_{ij} \left(\frac{\sigma_{ij}}{r_{ij}} \right)^{p_{ij}} - \left(\frac{\sigma_{ij}}{r_{ij}} \right)^{q_{ij}} \right] \exp \left(\frac{\sigma_{ij}}{r_{ij}-a_{ij}\sigma_{ij}} \right)\\ \phi_3\left(r_{ij}, r_{ik}, \theta_{ijk} \right) =& \lambda_{ijk} \epsilon_{ijk} \left[f_C(\delta) \delta \right]^2 \exp \left(\frac{\gamma_{ij}\sigma_{ij}}{r_{ij}-a_{ij}\sigma_{ij}} \right) \exp \left(\frac{\gamma_{ik}\sigma_{ik}}{r_{ik}-a_{ik}\sigma_{ik}} \right) \\ & \delta = \cos \theta_{ijk} - \cos \theta_{0ijk}\\\end{split}\]

where \(\phi_2\) and \(\phi_3\) are two-body and three-body terms. The cutoff of SW potential is \(a\cdot\sigma\). For some materials, such as borophene and transition metal dichalcogenides, some unnecessary angle types should be excluded in the three-body interaction by multiplying the cutoff function \(f_C(\delta)\). Here, for transition metal dichalcogenides, \(\delta_1\) is set to 0.25 and \(\delta_2\) is set to 0.35.

File format

This hybrid potential requires 2 potential files: ILP potential file and SW potential file. We have adopted the ILP file format that similar but not identical to that used by lammps. To identify the layers, it’s required to set two group_methods in model.xyz file. group_method 0 is used to identify the different layers and group_method 1 is used to identify different sublayers. One transition metal dichalcogenide layer has three sublayers, i.e., one \(MoS_2\) layer has one Mo sublayer and two S sublayers. For atoms in the same layer, the group_id of group_method 0 must be the same and for atoms in the same sublayer, the group_id of group_method 1 must be the same. Now this hybrid could be only used to simulate transition metal dichalcogenide homostructures (\(\mathrm{MX_2}\)), with M a transition metal atom (Mo, W, etc.) and X a chalcogen atom (S, Se, or Te).

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

potential <ilp file> <sw file>

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

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

sw_ilp is the name of this hybrid potential.
number of atom types is the number of atom types defined in the model.xyz. Here, this value must be set to 2 for transition metal dichalcogenide homostructures.
list of element is a list of all the elements in the potential.
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 \(\mathrm{MX_2}\), the ilp file is required to set as:

sw_ilp 2 M X
beta_MM alpha_MM delta_MM epsilon_MM C_MM d_MM sR_MM reff_MM C6_MM S_MM rcut1_MM rcut2_MM
beta_MX alpha_MX delta_MX epsilon_MX C_MX d_MX sR_MX reff_MX C6_MX S_MX rcut1_MX rcut2_MX
beta_XM alpha_XM delta_XM epsilon_XM C_XM d_XM sR_XM reff_XM C6_XM S_XM rcut1_XM rcut2_XM
beta_XX alpha_XX delta_XX epsilon_XX C_XX d_XX sR_XX reff_XX C6_XX S_XX rcut1_XX rcut2_XX

The sw file use the same atomic type list as the ilp file and just contains parameters of SW. The potential file reads, specifically:

A_MM B_MM a_MM sigma_MM gamma_MM
A_MX B_MX a_MX sigma_MX gamma_MX
A_XX B_XX a_XX sigma_XX gamma_XX
lambda_MMM cos0_MMM
lambda_MMX cos0_MMX
lambda_MXM cos0_MXM
lambda_MXX cos0_MXX
lambda_XMM cos0_XMM
lambda_XMX cos0_XMX
lambda_XXM cos0_XXM
lambda_XXX cos0_XXX