Clarification for the Tersoff potential (3-body terms)

Dear all,

I’m sorry about yet another (long) mail about the Tersoff potential parameters, as there have been so much already (most recently last week…). I write this mail because recently I found an error in my own potential file for BN, which was based on the work of Albe and Möller[1,2], which I used a couple of years ago. This made me look into the standard Tersoff potential files again to create a correct one, after which I was left confused. This comes mainly from subtle differences between ordering in SiC_Erhart-Albe.tersoff and BNC.tersoff . There are some other issues with the latter file I found, which I will mention below.

For the SiC_Erhart-Albe potential different coefficients are given for C-C, Si-Si and C-Si interactions (from [3]). For the two body interactions everything is fine and done according to the manual. For the three-body interactions the coefficients are taken like i,k,j. Example:
C Si C --> C-C coefficients
C C Si --> C-Si coefficients etc.
Meaning that k, the second element in the potential file is not really relevant for the 3-body interactions.

For BNC this is done differently, ignoring the C interactions one keeps interactions for B-N, B-B and N-N. Here the three-body terms are taken like i,j,k, which seems more logical. For example:
B N B --> B-N coefficients
B B N --> B-B coefficients etc.

To me it seems, maybe naively, that one of them should be wrong. Looking through the mail list and documentation I find some conflicting information. For example, the i,k,j scheme is mentioned in:
http://sourceforge.net/p/lammps/mailman/message/25688926/
while here it seems that it should be i,j,k:
http://sourceforge.net/p/lammps/mailman/message/32382114/
The documentation seems to be ambiguous about this fact, as it mentions that SiCSi and SiSiC will (in general) not be the same. This statement is true for both i,j,k ordering and i,k,j. For most cases, eg. perfect solids, it goes well either way. When one has defects like two interchanged atoms, these differences become important.

If somebody could clarify this it would be greatly appreciated.

Aidan can likely answer this.

Steve

The SiC potential format as well as the i, k, j schemes are incorrect. The third element is the atom that affects the “bonded” first and second elements, so the parameters for entry “A B C” should follow those for “A B X”.

However, the errors in the SiC potential does not matter since the bigA and bigB of the questioned entries are zero so they really have no effect at all.

Ray

Dear Ray,

Thank you for your answer. I agree with you that the errors do not matter for two-body interactions, as indeed bigA and bigB are zero (according to the manual they are always taken to be zero/ignored for two-body interactions when j is not equal to k). However, other parameters such as gamma, c, d etc. are non-zero and should be taken into account for three-body interactions.

I’ve been trying all evening to get sensible results for SiC with the parameters changed, but so far I’ve failed.

Guus

They should follow the first and second elements as I replied previously. But they have no effect if the bigA and bigB are zeros.

Ray

Dear Ray and others,

Below I attach, what then should be, the correct SiC parameters according to the scheme: Si-Si when i=j=Si, C-C when i=j=C-C and Si-C when i!=j. Values for n, beta, lambda2, B, lambda1, and A are zero for cases when j != k (note that setting these to non-zero values will not make any difference).

I did some test calculations with SiC and there are some differences with the default lammps version:
-For simple Energy vs lattice constant calculations an extra valley appears at a<a_0, where the default potential file already starts to repel. That the two versions of the potential file differ for smaller distances seems logical, as more atoms lie within the interaction range of each other, the differences in the 3-body terms then start to play a role.

-For a = a_0 both versions give the same energy, as indeed for this the terms for which A and B are zero are not relevant. However, when one considers a simple antisite defect (such as a C atom changed to a Si atom) the energy starts to differ with the two methods. Again, this is logical, as the relevant 3-body terms are different. Attached below a small input script to test this (note: one has to run this twice with different potential files to compare).

To conclude, if I interpret this correctly, I believe that the ordering matters and that the non-zero terms for j!=k do play a role in specific 3-body interactions when dealing with point defects. That A and B are zero for these entries is not important as those values are ignored by the code anyway. The other terms in those entries are not ignored, as can be seen by the defect example given above.

Best,

Guus

---------------------------------new SiC file--------------------------------

Si Si Si 1 0.114354 0 2.00494 0.81472 -0.259 1 1 1.538104933 219.5216243 2.82 0.14 2.833189287 2145.7128

Si Si C 1 0.114354 0 2.00494 0.81472 -0.259 0 0 0 0 2.82 0.14 0 0

Si C Si 1 0.011877 0 273987 180.314 -0.68 0 0 0 0 2.4 0.2 0 0

Si C C 1 0.011877 0 273987 180.314 -0.68 1 1 1.768074213 225.1894805 2.4 0.2 3.265633071 1779.361441

C C C 1 0.11233 0 181.91 6.28433 -0.5556 1 1 1.930900933 175.4266506 2 0.15 4.184262322 2019.844902

C C Si 1 0.11233 0 181.91 6.28433 -0.5556 0 0 0 0 2 0.15 0 0

C Si C 1 0.011877 0 273987 180.314 -0.68 0 0 0 0 2.4 0.2 0 0

C Si Si 1 0.011877 0 273987 180.314 -0.68 1 1 1.768074213 225.1894805 2.4 0.2 3.265633071 1779.361441

------------------------------test simulation for defects --------------------------

units metal
boundary p p p
atom_style atomic

variable b equal 4.36/2.0

lattice custom ${b} a1 1.0 1.0 0.0 a2 0.0 1.0 1.0 a3 1.0 0.0 1.0 &
basis 0.0 0.0 0.0 basis 0.25 0.25 0.25

region mybox block 0.0 10.0 0.0 10 0 1 units lattice
create_box 2 mybox

region myatomsbox block 0.0 10.0 0.0 10 0 1 units lattice
create_atoms 1 region myatomsbox basis 1 1 basis 2 2

pair_style tersoff
#change the potential file to your favourite
pair_coeff * * SiC_Erhart-Albe.tersoff Si C
mass 1 28.085
mass 2 12.011

neighbor 0.3 bin
neigh_modify delay 10
thermo_style custom step v_b temp pe ke etotal vol press pxx pyy pzz

run 0

set atom 142 type 1

run 0

set atom 143 type 2

run 0

set atom 142 type 2

run 0

The question you are asking is not really about LAMMPS, but rather it
is about the choices made by various developers of Tersoff potential
parameterizations. Ultimately the only way to answer these types of
questions is to examine the original computer code that was used to
develop the particular parameterization in question. With that said,
as far as I can tell from his papers, Karsten Albe always chose to
have the threebody parameters depend on i and k, but not j.* Hence
the current LAMMPS version of the file SiC_Erhart-Albe.tersoff is
correct and your new version is wrong. You correctly note that in the
file BNC.tersoff, a different choice is made. I have to assume this
is consistent with Tahir Cagin's original code, since the BNC.tersoff
file was contributed by his group. I hope this helps.

Aidan

* I am using the standard Tersoff notation here, where i is the
central atom (atom 1), j is a bonded neighbor (atom 2), and k is
another neighbor that contributes to the i-j bond-order parameter Bij.

Guus,

Different developers of Tersoff potentials make different choices. The
LAMMPS implementation supports many of them. The Tersoff potential
developed in the group of Carsten Albe, including the one defined by
SiC_Erhart-Albe.tersoff, have the threebody parameters depending on
the types of atoms i and k, but not j. Here i refers to the central
atom, j its bonded neighbor, and k another neighbor of i that
contributes to the bond-order parameter. For that reason, your
proposed change to SiC_Erhart-Albe.tersoff is wrong. In the case of
the potential defined by BNC.tersoff, threebody parameters depend on
the types of atoms i and j, but not k. I do not know if this is
consistent with the original code used by Tahir Cagin and co-workers
when they developed the BNC Tersoff potential, but I have to assume it
is, since BNC.tersoff was contributed by another researcher from Tahir
Cagin's group. I hope this helps.

Aidan