# How do LAMMPS read in Stillinger-Weber force field parameters?

Hi,

According to the description on manual https://docs.lammps.org/pair_sw.html, I think LAMMPS only use \lambda for 3-body interaction, and it should read \lambda only from the 3-body line (e.g. Si Si C) in the input *.sw file. Any change of \lambda in the 2-body line (e.g. Si C C) shouldn’t affect the final result.

Is my understanding correct? How does LAMMPS read in the *.sw parameters actually? Since I have tried different values of \lambda in 2-body line and it changes the final result.

Thanks,

Zui

All lines in the potential files are describing 2-body and 3-body interactions. When looking up parameters for 2-body interactions, some entries are ignored, because they are redundant, but you cannot do that for the 3-body parameters.

That is what I would expect.

\begin{aligned} E & =\sum_i \sum_{j>i} \phi_2\left(r_{i j}\right)+\sum_i \sum_{j \neq i} \sum_{k>j} \phi_3\left(r_{i j}, r_{i k}, \theta_{i j k}\right) (1)\\ \phi_2\left(r_{i j}\right) & =A_{i j} \epsilon_{i j}\left[B_{i j}\left(\frac{\sigma_{i j}}{r_{i j}}\right)^{p_{i j}}-\left(\frac{\sigma_{i j}}{r_{i j}}\right)^{q_{i j}}\right] \exp \left(\frac{\sigma_{i j}}{r_{i j}-a_{i j} \sigma_{i j}}\right) (2)\\ \phi_3\left(r_{i j}, r_{i k}, \theta_{i j k}\right) & =\lambda_{i j k} \epsilon_{i j k}\left[\cos \theta_{i j k}-\cos \theta_{0 i j k}\right]^2 \exp \left(\frac{\gamma_{i j} \sigma_{i j}}{r_{i j}-a_{i j} \sigma_{i j}}\right) \exp \left(\frac{\gamma_{i k} \sigma_{i k}}{r_{i k}-a_{i k} \sigma_{i k}}\right) (3) \end{aligned}

Then I wonder how to set up \lambda in 2-body line, since I only fitted the \lambda for 3 body, as shown in the above formula (3). Seems that we only need this parameter in 3-body term? How to set up \lambda for 2-body line is correct in this case?

I understand that for \sigma, a and \gamma of 3-body interaction, I should define them in 2-body line, because there are 2 pairs of {\sigma, a, \gamma} for 3-body ineractions, while we only have one pair of parameter in one 3-body line. So LAMMPS read these parameters from the 2-body line. But I don’t understand how does \lambda work.

Thanks again,
Zui

I repeat, there are no 2-body lines. All lines are required for 3-body interactions.

You have one value of \lambda for each permutation of atom types for the triple “i” “j” “k”. For a two element system, there are 8 of them. The potential file correspondingly has to have 8 entries.

Hi @Zui,

I think there is a misunderstanding of what you mean by “3-body line” compared to what the equations are about ans what you understand from @akohlmey’s answer.

The Si C C line needs parameters for equation (2) and equation (3) for pair and triplet of atoms. The format assumes that the 1st atom is of type Si and the 2nd and 3rd of type C. The Si Si C line does exactly the same for when the type of the second atom is Si etc. However in the second case, the 2 body parameters are not read because that would be redundant with the 1st case, which already defined them.

So only the 2 body parameters of certain lines are used to define the internal 2 body parameters because it would just be redundant otherwise. The fact that the two body parameters comes from 122 elements is just an arbitrary convention. This does not define “2 body lines” and “3 body lines”, it is just that some parameters are skipped on some lines to avoid unnecessary operations and remove sources of confusion. This is emphasised by the fact that, as stated in the documentation, since the input order in the file is arbitrary, parameters still need to be consistent with regard to elements permutation (2 body 122 and 211 should be equal, same for 3 body 121 and 112 etc.)

1 Like

Hi @Germain,

Thanks, it is more clear now. But the force field still doesn’t work for my project. Here I want to describe my problem in a more detailed way. I guess there’s still sth wrong with how I enter the parameters.

I want to use Stillinger-Weber + Coulomb force field to relax an InAs quantum dot. The origin force field fitted in paper is written in GULP format. Here’s the link for this paper: https://journals.aps.org/prb/abstract/10.1103/PhysRevB.85.041306 and the origin file for GULP https://www.chemie.uni-hamburg.de/institute/pc/arbeitsgruppen/bester/downloads.html.

The form of this force field:

The potential proposed in this work consists of short-range two-body terms V_2(i,j), a three-body term V_3(i,j,k), and a long-range Coulomb term V_C(i,j).

\begin{aligned}V= & \sum_{i<j}^{j \in 1 N N} V_2(i, j)+\sum_{i<k}^{k \in 2 N N} V_2(i, k) +\sum_{i<j<k} V_3(i, j, k)+\sum_{i<j} V_C(i, j)\end{aligned}

The two body terms, which describe the short-range bond-stretching interactions between the first nearest neighbor (1NN) and second nearest neighbor (2NN) atoms, have the following form:

V_{2}(i,j)=\left\{ \begin{array}{l}A\left(\frac{B}{r_{ij}^{4}}-1\right)\exp\left(\frac{\rho}{r_{ij}-b}\right),r_{ij}<b\\0,r_{ij}\geqslant b\end{array}\right.

where r_{ij} denotes the distance between atom i and j , and b represents the cutoff distance of the interaction.

The three-body term corresponds to the contribution of an angle distortion with respect to the ideal tetrahedral angle and is given by:

V_3(i, j, k)=h\left(r_{i j}, r_{i k}\right)+h\left(r_{j i}, r_{j k}\right)+h\left(r_{k i}, r_{k j}\right)

where,

h\left(r_{i j}, r_{i k}\right)=\lambda \exp \left(\frac{\eta}{r_{i j}-b}+\frac{\eta}{r_{i k}-b}\right)\left(\cos \theta_{j i k}+\frac{1}{3}\right)^2
V_{C}(i,j)=\frac{Z_{i}^{*}Z_{j}^{*}}{\left|r_{ij}\right|}

The fitted parameters in paper:

\rho = 2.2819 for InAs read from GULP input.

Compared this with LAMMPS form:

\begin{aligned} E & =\sum_i \sum_{j>i} \phi_2\left(r_{i j}\right)+\sum_i \sum_{j \neq i} \sum_{k>j} \phi_3\left(r_{i j}, r_{i k}, \theta_{i j k}\right) \\ \phi_2\left(r_{i j}\right) & =A_{i j} \epsilon_{i j}\left[B_{i j}\left(\frac{\sigma_{i j}}{r_{i j}}\right)^{p_{i j}}-\left(\frac{\sigma_{i j}}{r_{i j}}\right)^{q_{i j}}\right] \exp \left(\frac{\sigma_{i j}}{r_{i j}-a_{i j} \sigma_{i j}}\right) \\ \phi_3\left(r_{i j}, r_{i k}, \theta_{i j k}\right) & =\lambda_{i j k} \epsilon_{i j k}\left[\cos \theta_{i j k}-\cos \theta_{0 i j k}\right]^2 \exp \left(\frac{\gamma_{i j} \sigma_{i j}}{r_{i j}-a_{i j} \sigma_{i j}}\right) \exp \left(\frac{\gamma_{i k} \sigma_{i k}}{r_{i k}-a_{i k} \sigma_{i k}}\right) \end{aligned}

We can get the parameter file for LAMMPS:

# Parameters for zincblende InAs from
# Table II of P. Han, G. Bester, PRB, 83, 174304 (2011)
#
# Note that there is a mistake in this reference (i.e. rho is
# missing in the definition of the potential and parameter table).
# The value for rho is found in the GULP input files at:
#
# All other parameters defined in the reference are correct.
#
# The twobody ik pair parameters are entered on the ikk line.
# The threebody ijk parameters (where i is the central atom)
# are defined on the ijk line.
#
# These entries are in LAMMPS "metal" units: epsilon = eV;
# sigma = Angstroms; other quantities are unitless
#
# for all interactions:
# eps = 1
# cos(theta) = -1/3
# p = 4
# q = 0
# other LAMMPS parameters are given in the reference as
# sigma = rho
# a = b / rho
# lambda = lambda
# gamma = eta / rho
# A = A
# B = B / rho^4
#
#               eps          sigma            a            lambda          gamma       cos(theta)         A              B              p              q             tol
In In In   1.000000e+00   2.281900e+00   2.471625e+00   2.794020e+01   1.257592e+00  -3.333333e-01   1.019500e+00   1.034534e+00   4.000000e+00   0.000000e+00   0.000000e+00
As As As   1.000000e+00   2.281900e+00   2.471625e+00   2.794020e+01   1.257592e+00  -3.333333e-01   1.019500e+00   1.034534e+00   4.000000e+00   0.000000e+00   0.000000e+00

As In In   1.000000e+00   2.281900e+00   1.799991e+00   2.794020e+01   1.257592e+00  -3.333333e-01   4.407400e+00   1.782599e+00   4.000000e+00   0.000000e+00   0.000000e+00
In As As   1.000000e+00   2.281900e+00   1.799991e+00   2.794020e+01   1.257592e+00  -3.333333e-01   4.407400e+00   1.782599e+00   4.000000e+00   0.000000e+00   0.000000e+00

In In As   1.000000e+00   2.281900e+00   1.799991e+00   2.794020e+01   1.257592e+00  -3.333333e-01   4.407400e+00   1.782599e+00   4.000000e+00   0.000000e+00   0.000000e+00
In As In   1.000000e+00   2.281900e+00   2.471625e+00   2.794020e+01   1.257592e+00  -3.333333e-01   1.019500e+00   1.034534e+00   4.000000e+00   0.000000e+00   0.000000e+00
As In As   1.000000e+00   2.281900e+00   2.471625e+00   2.794020e+01   1.257592e+00  -3.333333e-01   1.019500e+00   1.034534e+00   4.000000e+00   0.000000e+00   0.000000e+00
As As In   1.000000e+00   2.281900e+00   1.799991e+00   2.794020e+01   1.257592e+00  -3.333333e-01   4.407400e+00   1.782599e+00   4.000000e+00   0.000000e+00   0.000000e+00


Is this correct? I think the problem may come from \epsilon and cutoff tol.

Is tol in LAMMPS defined by a \times \sigma ? (It is written in manual of LAMMPS, but the symbol is not clear…) And how do LAMMPS read different \epsilon_{ij} and \epsilon_{ijk}, since these two parameters are defined in the same line at the same place ?

Thanks,

Zui

In LAMMPS, the parameters for the ij pair are taken from the ijj entry.

1 Like

It is very hard for people on a forum to tell on-the-fly if your conversion from GULP inputs to LAMMPS is correct without redoing the process themselves. And to do so seriously, one basically needs to do the reading of the papers, the codes documentation and check every parameters and equations equivalence, which is tedious. More so when making test case and input files to compare outputs. This is your part of the work and largely outside the scope of the forum. There is no workaround, we’ve all been there. Note that “It does not work for my project” is a meaningless sentence to me since I don’t know what your project is and what makes you say things do not work, because you never said those things.

The meaning of the tol parameter and the way LAMMPS reads the inputs are already described in details in the documentation page you linked in your first message. Best I can advise you is to read it more carefully take notes and even rewrite it with your own words until you get a proper grasp of what is going on. Most questions you asked so far find their answer there.

Finally, as a general advice, you also have to compare the implementation of Stillinger-Weber equations and notation between the paper, and LAMMPS, but also from GULP. Because that is the code your input file was made for. You will get more valuable insights by asking people with which you can interact directly and know what you’re up to than people on a forum.

Side note: The paper you provided a link to is not the one publishing the model.

I see， thank you.