Question about pair_sw.cpp for assistance with Garofalini implementation

Hello everyone,

I am a graduate student at the University of California, Merced.

I have been working on an implementation of the Garofalini interatomic potential (amorphous silica) for possible addition to the growing list of potentials that have been made available for use with LAMMPS. I have been using the Stillinger-Weber potential as the closest relevant starting point, as it addresses 3-body interactions. In attempting to completely understand the functioning of the Stillinger-Weber code, I’m trying to grasp more thoroughly the decomposition of forces as represented by the lines:

fj[0] = delr1[0]*(frad1+csfac1)-delr2[0]facang12; [pair_sw.cpp (591)]
fj[1] = delr1[1]
(frad1+csfac1)-delr2[1]facang12; [pair_sw.cpp (592)]
fj[2] = delr1[2]
(frad1+csfac1)-delr2[2]*facang12; [pair_sw.cpp (593)]

I have been studying the Stillinger-Weber code intensively without being able to fully understand the elements of the force components above. Please let me know if Dr. Thompson or anyone else would be willing to expand on these expressions for the benefit of this project.

Thank you,
Benjamin Doblack

Hi Benjamin,

This piece of code is more or less a transcription of the mathematical derivation of the forces. It’s hard to read, but it is really fairly straightforward. Starting with the threebody energy expression E (that I called facrad), the force Fj on atom j is obtained by taking the derivative of E with respect to vector Rj. By the chain rule, this gives two terms, a radial term and an angular term. The first is proportional to dE/d|Rij| (variable frad1 below) and is oriented along vector Rij. The second is proprtional to dE/dCosThetaijk and has components along vectors Rij (variable csfac1 below) and Rik (variable facang12 below). By permuting j and k, you can get the same expressions for the force Fk on atom k. Because E is invariant w.r.t. translation of all three atoms, Fi = -Fj - Fk.


p.s. a good warm-up exercise is to convince yourself that d|Rij|/dRj = Rij/|Rij|, where vector Rij = Rj - Ri.

Hi Aidan,

Thank you kindly for your time. This helped me figure out some invalid assumptions I made in generating the partials for Fj and Fk.