Looping over all atoms

Hey all,

I’m trying to make a pair style that loops over all atoms (two sums over “i” and “j”, "i ~= j ") but I think most pair styles in LAMMPS loop over half the atoms ("two sums over “i” and “j”, “i < j”). Does LAMMPS do this? Here is the loop code for the Morse potential:

// loop over neighbors of my atoms

for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];

for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_lj = special_lj[sbmask(j)];
j &= NEIGHMASK;

delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delxdelx + delydely + delz*delz;
jtype = type[j];

I can’t find anywhere that implies this is looping over only half the atoms - so are all the atoms being considered?

Hey all,

I'm trying to make a pair style that loops over all atoms (two sums over "i"
and "j", "i ~= j ") but I think most pair styles in LAMMPS loop over half
the atoms ("two sums over "i" and "j", "i < j"). Does LAMMPS do this? Here
is the loop code for the Morse potential:

unless you are programming some kind of manybody potential, where
you'll need to look at neighbors of neighbors, there is no real need
to loop over all atoms. just take pair i-j and apply the result to
both atoms, i.e. the pair j-i is effectively the same interaction with
the opposite sign on the forces. that is what is happening in the
example below.

as you also see below, LAMMPS doesn't loop over all atoms, but just
over neighbors of atoms and for that it uses a list of neighbors. thus
the way to change this behavior is to request a different kind of
neighbor list. please have a look at src/MANYBODY/pair_sw.cpp for one
of the least complex manybody potentials in LAMMPS, that uses a full
neighbor list, instead of the (default) half neighbor list. please
also note, that the pairwise additive term is still computed over half
of the neighbors (newton's third law is a beautiful thing, it gives
you a 2x speedup (almost) for free).

axel.

I’m actually coding a Taylor expansion potential, so the displacement of every single atom must be considered even for the 2 body term. In that sense, it isn’t even a pairwise potential and that’s why I need to loop over all atoms.

Where in these examples is the half neighbor list being enforced? Here’s the SW loop:

// loop over full neighbor list of my atoms

for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
itag = tag[i];
itype = map[type[i]];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];

// two-body interactions, skip half of them

jlist = firstneigh[i];
jnum = numneigh[i];

for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtag = tag[j];

So by setting jnum = numneigh[i], we only consider half the neighbors? Later on in the code:

jnumm1 = jnum - 1;

for (jj = 0; jj < jnumm1; jj++) {
j = jlist[jj];
j &= NEIGHMASK;

This time, the loop is over “jnum - 1”. So the difference between looping over the half neighbor list and full neighbor list is looping over “numneigh[i]” and “numneigh[i] - 1”. It looks like firstneigh[i] and numneigh[i] come from the following pointers:

numneigh = list->numneigh;
firstneigh = list->firstneigh;

I can’t seem to find the “list” class… Sorry for the confusion, I’m just trying to understand how to loop over a full neighbor list versus only the half neighbor list, and what that correlates to in the source code.

Drew Rohskopf | Graduate Research Assistant

Atomistic Simulation & Energy Group

Georgia Institute of Technology

(404)403-0313

I'm actually coding a Taylor expansion potential, so the displacement of
every single atom must be considered even for the 2 body term. In that
sense, it isn't even a pairwise potential and that's why I need to loop over
all atoms.

your description sounds quite inconsistent and confused to me. in part
i don't have any experience with a "Taylor expansion potential" and
what it encompasses. most likely your description is not using the
right "MD slang". regardless, with *any* neighbor list you are *not*
looping over *all* atoms, but only over neighbors within the given
cutoff plus skin margin. looping over all atoms would be really tricky
in LAMMPS, since LAMMPS uses domain decomposition and each processor
only owns a subset of the total system plus some "ghost atoms" from
its neighbors.

Where in these examples is the half neighbor list being enforced? Here's the
SW loop:

  // loop over full neighbor list of my atoms

  for (ii = 0; ii < inum; ii++) {
    i = ilist[ii];
    itag = tag[i];
    itype = map[type[i]];
    xtmp = x[i][0];
    ytmp = x[i][1];
    ztmp = x[i][2];

    // two-body interactions, skip half of them

    jlist = firstneigh[i];
    jnum = numneigh[i];

    for (jj = 0; jj < jnum; jj++) {
      j = jlist[jj];
      j &= NEIGHMASK;
      jtag = tag[j];

So by setting jnum = numneigh[i], we only consider half the neighbors? Later

no. keep on reading. there is a section where it does:

      jtag = tag[j];
      if (itag > jtag) {
        if ((itag+jtag) 2 == 0\) continue; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;\} else if \(itag &lt; jtag\) \{ &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;if \(\(itag\+jtag\) 2 == 1) continue;
      } else {
        if (x[j][2] < ztmp) continue;
        if (x[j][2] == ztmp && x[j][1] < ytmp) continue;
        if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue;
      }

this reduces the requested full neighbor list to a half neighbor list.

on in the code:

    jnumm1 = jnum - 1;

    for (jj = 0; jj < jnumm1; jj++) {
      j = jlist[jj];
      j &= NEIGHMASK;

This time, the loop is over "jnum - 1". So the difference between looping
over the half neighbor list and full neighbor list is looping over
"numneigh[i]" and "numneigh[i] - 1". It looks like firstneigh[i] and
numneigh[i] come from the following pointers:

no. if you look more carefully you should see that this is to do a
loop over all i-j-k atom triples without double counting.

  numneigh = list->numneigh;
  firstneigh = list->firstneigh;

I can't seem to find the "list" class... Sorry for the confusion, I'm just
trying to understand how to loop over a full neighbor list versus only the
half neighbor list, and what that correlates to in the source code.

you probably also need to update your knowledge of C++ and practice
reading C++ code some more.
the loop construct is pretty much the same, regardless of the kind of
neighbor list.

as mentioned in my previous response, the kind of neighbor list is set
up in the neighbor list request. for pair styles, this is done in the
Pair::init_style() (virtual) method, i.e. by default
Pair::init_style() will request a half neighbor list with the selected
treatment for interactions between local and ghost atoms as set by the
newton command.
in PairSW::init_style() this function is overridden and the request
changed from a half list to a full list (for which the newton pair
setting is - for obvious reasons - irrelevant).

keep on reading some more code, experiment with a simple pair style
and follow the flow of control. remember that in C++ this can be
scattered across the class hierarchy, i.e. a lot of the "bookkeeping"
and "management" as well as default settings are handled in the Pair
base class.

i suggest as an exercise, you should try and implement a pair style
morse/full, that uses a full neighbor list instead of a half list and
do this by making a copy of the existing morse and then changing it as
necessary, while making certain to get the same result. that exercise
should hopefully provide you with the necessary insight.

axel.