Hybrid NEP+ILP potential

The hybrid NEP + ILP potential [Bu2025] in GPUMD combines the neuroevolution potential (NEP), [Fan2022b] (NEP3), and [Song2024] (NEP4), for intralyer interactions and the interlayer potential (ILP) [Ouyang2018] [Ouyang2020] for interlayer interactions to simulate van der Waals materials. Now this hybrid potential supports to simulate homo- and heterostructures based on graphene, \(h\)-BN and transition metal dichalcogenides (TMDs) layered materials. The nep potential here doesn’t support USE_TABLE flag to accelerate now.

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

\[U_i^{\mathrm{NEP}} = \sum_{\mu=1}^{N_\mathrm{neu}}w^{(1)}_{\mu}\tanh\left(\sum_{\nu=1}^{N_\mathrm{des}} w^{(0)}_{\mu\nu} q^i_{\nu} - b^{(0)}_{\mu}\right) - b^{(1)},\]

More details of NEP potential are in Neuroevolution potential. Note that in hybrid NEP + ILP potential, the NEP potential will just compute the intralayer interactions.

File format

This hybrid potential requires 3 kinds of files: one for ILP potential, one for NEP potential and the other for mapping NEP potential to groups in model file. We have adopted the ILP file format that similar but not identical to that used by lammps. The NEP potential file is not required to modify, while to make the ILP and NEP potentials identify the layers, it’s required to set some groups in model.xyz file.

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

potential <ilp file> <nep map file>

where ilp file and nep map file are the filenames of the ILP potential file and NEP mapping file.

ilp file is similar to other empirical potential files in GPUMD. But in addition, ILP uses different group_ids to identify the different layers, so you need to add two group_methods in ilp file:

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

  • nep_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 (can be in any order).

  • group_method for layers is the group_method set in model.xyz to identify different layers. For example, monolayer graphene and monolayer \(\mathrm{MoS}_2\) are both single layer so for the atoms in each layer the group_id of group_method for layers are the same.

  • group_method for sublayers is used to identify the different sublayers. For example, monolayer graphene contains one sublayer while monolayer \(\mathrm{MoS}_2\) contains three sublayers, one Mo sublayer and two S sublayers. For the atoms in each sublayer the group_id of group_method for sublayers 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Å.

nep_map_file can map one or more NEP potential files to different layers. The setting is as:

<group_method for layers> <number of NEP files> <list of NEP files>
<number of groups>
<NEP_id for group_0>
<NEP_id for group_1>
...

  • group_method for layers is the same as the setting in ilp file.

  • number of NEP files is the number of NEP files used in your simulation.

  • list of NEP files is a list of all the NEP filenames. Note that the first file will be identified as NEP_0 and then NEP_1 and so on.

  • number of groups is the number of groups in group_method for layers.

  • The last number of groups lines map the NEP to each group. If NEP_id for group_0 is set to 0, the intralayer interactions between atoms within group_id 0 are computed by the first NEP file (NEP_0) in list of NEP files. If set to 1, then computed by the second NEP file (NEP_1) and so on.

Examples

Example 1: bilayer graphene

Assume your have three files: ILP potential file (C.ilp), NEP potential file (C.nep) and NEP mapping file (map.nep). The potential setting in run.in file is as:

potential C.ilp map.nep

Assume that the first line in C.nep is:

nep3 1 C

and group_method 0 is used to identify the different layers. Then C.ilp is required to set as:

nep_ilp 1 C
0 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 first 0 in the second line represents ILP potential uses group_method 0 to identify different layers. The second 0 represents group_method 0 is used to identify the sublayers. For the system with only graphene and \(h\)-BN, just set it the same as the previous number.

Then, map.nep file required to set as:

0 1 C.nep
2
0
0

The first 0 in the first line represents NEP potential uses group_method 0 to identify different layers. The next 1 represents there is just one NEP potential file. The number in the second line represents there are two groups in the group_method 0. The last two lines represent the group_0 and group_1 in group_method 0 will use C.nep potential file (NEP_0).

Example 2: bilayer \(h\)-BN / \(\mathrm{MoS}_2\)

Assume your have four files: ILP potential file (BNMoS.ilp), NEP potential files (BN.nep, MoS.nep) and NEP mapping file (map.nep). The potential setting in run.in file is as:

potential BNMoS.ilp map.nep

Assume the first line in BN.nep is:

nep4 2 B N

and in MoS.nep is:

nep4 2 Mo S

We also assume the group_method 0 is used to identify the different layers and group_method 1 is used to identify the different sublayers for ILP. In group_method 1, atoms in the sublayers of Mo and S should be set as the different group_id. Then BNMoS.ilp is required to set as:

nep_ilp 4 B N Mo S
0 1
beta_BB alpha_BB delta_BB epsilon_BB C_BB d_BB sR_BB reff_BB C6_BB S_BB rcut1_BB rcut2_BB
beta_BN alpha_BN delta_BN epsilon_BN C_BN d_BN sR_BN reff_BN C6_BN S_BN rcut1_BN rcut2_BN
beta_BMo alpha_BMo delta_BMo epsilon_BMo C_BMo d_BMo sR_BMo reff_BMo C6_BMo S_BMo rcut1_BMo rcut2_BMo
beta_BS alpha_BS delta_BS epsilon_BS C_BS d_BS sR_BS reff_BS C6_BS S_BS rcut1_BS rcut2_BS
...
beta_SS alpha_SS delta_SS epsilon_SS C_SS d_SS sR_SS reff_SS C6_SS S_SS rcut1_SS rcut2_SS

Assume group_id of \(\mathrm{MoS}_2\) is 0 and of \(h\)-BN is 1. Then map.nep file is set as:

0 2 BN.nep MoS.nep
2
1
0

The 1 in the third line means group_0 (\(\mathrm{MoS}_2\)) uses MoS.nep potential file (NEP_1) and the last 0 means group_1 (\(h\)-BN) uses BN.nep potential file (NEP_0).