In the Tersoff potential for alloys, the bond order term requires the calculation of the cutoff function using the i-k parameters and separation (of question here is interactions like this: i-j-k = Si-Si-C). However, the way the potential input file is written, the values for the cutoffs (R and D) provided should be the i-j ones, and there is no separate value for the i-k set (and the code doesn’t seem to hunt these down in the parameter file.

The provided files seem to just use the same cutoff for the Si-Si-C type interaction as would be used for Si-C-Si - but this doesn’t seem to hold based on the way that Tersoff wrote the formulation. The Si-Si cutoff values are 2.85 +/- .15 and the Si-C cutoff is 2.36 +/- .15 - so I would figure this could lead to some differences when the structure begins to deform a fair bit. In my past simulations I had been using the 2.85 value in the Si-Si-C interactions. Curious if anyone knows if this makes a big difference in behavior - or if it may be worth adding in the extra cutoff parameters to the parameter file line?

The way code is written, for a cluster, the two body parameters come from i-j-j set and 3 body parameters come from i-j-k set.

For eg. Si-Si-C interaction, 2-body parameters come from Si-Si-Si set and 3-body parameters come from Si-Si-C. So, for Si-Si-C cluster, R (2.36) for i-k comes form Si-Si-C set, but for R i-j (2.85) it comes from Si-Si-Si set.

For Si-C-Si interaction, 2-body parameters will come from Si-C-C set and 3-body parameters will come from Si-C-Si. So, for this cluster, R i-k will come from Si-C-Si which is 2.85. But 2-body interaction Ri-j will come from Si-C-C set which is 2.36.

I think your confusion may be due to the fact that R and D are both two-body
and three-body parameters, but they can have different values in both cases.
The statements that Rutupurna made are all correct. The key lines of code
confirming this are:

This behavior matches the equation for zeta_ij given in the pair_tersoff
documentation, which in turn are consistent with the Tersoff alloy paper,
Tersoff_2 in the documentation.

I will add that it is possible to think up modified expressions for zeta_ij
that can not be handled by the current code, e.g. combining the ij and ijk
cut-offs in a non-linear way. Something like this happened for
Stillinger-Weber. The solution is quite simple; just add any extra
parameters to the argument list for zeta() as needed: