Embedding Tersoff interatomic potential parameters in Gulp

GULP
1

Hi, I’m very new to Gulp and am trying to run some simulations on W and WC, employing the Tersoff potential. I had a look at the help doc of this package to insert the parameters required i.e.:
boonebody
: atom1 alphaA alphaB nA nB
botwobody
: atom1 atom2 A B za zb rtaper rmax
borepulsive
: atom1 <omega_ik> lambda
boattractive theta
: atom1 <omega_ik> lambda c d h
By comparison between the equation used in the Gulp and the one provided by Juslin et al. (Journal of applied physics 98, 123520 (2005)), I considered za = λ and zb =µ.
Now I have these questions:
1: What are rtaper and rmax exactly? In example 77 they are in fact the values of R and S respectively, reported by Munetoh et al (Comp. Mat. Sci., 39, 334 (2007)), but in example 28 they don’t match with the values of R and S reported for Si. In this example rmax=cutd which makes sense but not sure how the value of rtaper is set here.
2: What does lambda exactly stand for in borepulsive and boattractive theta? I was expecting to have the value of λ in the repulsive term and the value of µ in the attractive term for this flag, but in example 28, it’s got the value of µ for both borepulsive and boattractive theta and in example 77, it’s got a zero value for both borepulsive and boattractive theta.

Any advice is highly appreciated!!!

Regards,

2

The best way to understand the input for Tersoff potentials is to take one of the library files that comes with GULP (such as tersoff.lib) and compare to the original paper. This way any variation in how the equations may be written (in an equivalent but different form) can be seen, as well as how the parameters were translated. I’d recommend starting with the original Tersoff papers since this is the reference point.

  1. rmax is the cutoff of the potential and rtaper is the distance at which the taper function begins.
  2. The lambda values are as per the equations in the help text:

For case without “theta” sub-option;

zeta = Sum(rik) [f(rik).omega_ik.exp(lambda**m.(rij-rik)**m)]

else;

zeta = Sum(rik) [f(rik).g(theta).omega_ik.exp(lambda**m.(rij-rik)**m)]

So lambda is an exponent in inverse Angstroms that controls the exponential term (raised to the power m).
Regards,
Julian

1 Like