For most pair potentials, there are just 2 places in pair*.cpp where force

is computed. In the compute() function and the single() functiion. Pair lj/cut

has a few more versions, b/c it can be used by rRESPA in some additional modes.

You might ask, why not make every pair style have exactly one place where

force is computed, i.e. a single() routine that everything else calls. This

would be inefficient when you want to compute forces for a long list of

atom pairs at once, which is what compute() does.

For pair SW, it has a 2-body loop and a 3-body loop in compute(). If that

isn't similar to what you need, then you can write something else. Pair SW

uses standard neighbor lists to find 2-body and 3-body interactions. If you

need something different, you'll have to explain what it is.

Steve

'Very much appreciate your reply. Yes I recognized the two and three body

loops.

Okay, specifics: for my two-body I want V=(qi*qj)/r + A exp(-r/B) where A

and B are factors; and q is a charge.

For the 3-body part V= L*exp[g_ji/(r_ij-r0ij)+g_ik/(r_ik-r0ik)] *

[(cos(tjik)-cos(t0jik))*sin(tjik)*cos(tjik)]^2

and

V=L*exp[g_ji/(r_ij-r0ij)+g_ik/(r_ik-r0ik)] * ((cos(tjik)-cos(t0jik))^2,

where L and g are factors

Here a different potential function applies for different triplets.

Upon inspecting the two-body SW code I found that there is an "extra" 1/r in

fforce (what does ff.. mean?) when compared to the negative of the derivate

of the potential function. I don't know how to explain that. There is also

another computation for "eflag" which I don't know what it's for but I'll

probably need to change /that/ as well.

I haven't analyzed the three-body code as closely yet but, again, I don't

know how it works. It looks like I'll need two modifications of the SW code

to accommodate both 3-body functions.