I try to add a new manybody potential into LAMMPS. I have to calculate the coordination number before calculating the atom force. If I directly define the coordination number to be a value, the code has no problem. But if I calculate it through a neighbor list. Some of j will be bigger than the total atom number and cause some problems during calculating the coordination number(I use ppp boundary condition). I think it related to the periodic conditions. But I don’t know how to avoid this situation.
Here is part of my code:
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// 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];
// printf("i %d f\n",i,zi);
// two-body interactions, skip half of them
jlist = firstneigh[i];
jnum = numneigh[i];
zi=0;
//////////////////////////////////The coordination number of atomi
for (jj = 0; jj < jnum; jj++)
{
j = jlist[jj];
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
r1 =sqrt( delx*delx + dely*dely + delz*delz);
if(r1<=(R-D)) zi+=1.0;
else if(r1<(R+D))
{ c1=(r1-R+D)/D;
zi+=1.0-0.5*c1+sin(PI*c1)/(2.0*PI); }
else zi+=0.0;
}
/////////////////////////////////////
for (jj = 0; jj < jnum; jj++)
{
j = jlist[jj];
j &= NEIGHMASK;
jtag = tag[j];
if (itag > jtag) {if ((itag+jtag) 2 == 0) continue;}
else if (itag < jtag) {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;
}
jtype = map[type[j]];
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delxdelx + delydely + delzdelz;
ijparam = elem2param[itype][jtype][jtype];
if ((rsq - params[ijparam].cut2sq)>(-1e-6)) continue;
///////////////////////////////////////The coordination number of atomj
llist = firstneigh[j];
lnum = numneigh[j];
zj=0;
for (ll = 0; ll < lnum; ll++)
{
l = llist[ll];
delx = x[j][0] - x[l][0];
dely = x[j][1] - x[l][1];
delz = x[j][2] - x[l][2];
r1 =sqrt( delxdelx + delydely + delzdelz);
if(r1<=(R-D)) zj+=1.0;
else if(r1<(R+D))
{ c1=(r1-R+D)/D;
zj+=1.0-0.5c1+sin(PIc1)/(2.0*PI);
}
else zj+=0.0;
}
Someone with more energy than me may want to
look at your code and debug it. But you are
the best person to do that. If it is not calculating
what you want, put in print statements, run a simple
geometry, etc. And figure out what it is doing
wrong. There is no simpler path to enlightenment.
Steve
I have debugged it for many times and the problem is that sometimes j (the atom id of the neighor list of i ) is bigger than the total atom number because of the periodic boundary conditions and it causes problems while calculating the coordination number of j atom. My question is that how to avoid j to be bigger than the total atom number in periodic boundary conditions. Thank you.
Lisa
2011/5/18 Steve Plimpton <[email protected]>
I have debugged it for many times and the problem is that sometimes j (the
atom id of the neighor list of i ) is bigger than the total atom number
because of the periodic boundary conditions and it causes problems while
calculating the coordination number of j atom. My question is that how to
avoid j to be bigger than the total atom number in periodic boundary
conditions. Thank you.
lisa,
i suspect that you don't really understand
your problem. what you describe is just the
symptom (and i suspect that you may be
even wrong with that).
for example, i don't see an provisions in your
code whether you have newton's third law enabled
or not, which is an essential factor for how atoms
are listed in as neighbors.
this is actually not overly complicated, but requires
an approach different from "i modify the code until
my problems go away".
there are not many alternatives to single stepping
your execution with a debugger and/or printing out
all details of all variables (you may want to try out
the DDD debugger frontend to gdb, which is good
for watching arrays).
axel.
I am using gdb to debug my program. I only know how to check the place
that causes segment fault. How can I set the breakpoint in the
potential file or check the neighbor list formation using gdb? (I am
using a remote cluster and does not support graphics)
I am using gdb to debug my program. I only know how to check the place
that causes segment fault. How can I set the breakpoint in the
potential file or check the neighbor list formation using gdb? (I am
gdb allows to inspect your data in many ways.
neither of us here can give you a tutorial in good
programming or debugging and how individual
tools work. it is too time consuming and over
e-mail also very complicated. it should be, however,
identify somebody local that can help you, in
case you cannot teach new skills to yourself.
using a remote cluster and does not support graphics)
program development is _much_ easier on a local machine.
that all being said, the key problem seems to be
that you are guessing what what is what and that
you are guessing wrong and thus your code is wrong.
if my comment on newton's 3rd law, for example,
doesn't ring a bell, then you are missing some fundamental
knowledge on neighbor lists (or verlet lists). that is
something that requires that you stick your nose
into the literature (textbooks, original papers) discussing
these for a while and then perhaps also read through
steve's paper on the parallelization of lammps and
_then_ get back to writing your pair style again.
if i understand correctly what you want to do, then you'll
need to either add some MPI communication into your
pair style or require that newton's third law is turned off.
axel.
J atoms bigger than nlocal are ghost atoms on that proc.
You presumably don't want to calculate the coord # of
ghost atoms, so your code logic should reflect that.
If you need the coord # of ghost atoms later in your
potential, then you should probably communicate it,
similar to the way the EAM potentials communicate
the accumlulated density on ghost atoms, via calls
to communicate.
Steve
Thank you. I realize that I was thinking on the wrong way. I am not sure if I understand it right now.
On one processor, the atom i has a neighbor j, which is on another processor. This processor can get the position information of the atom j, but not the neighbor list information. So if I just calculate the diatance between i and j and calculate the force, it can work thoughly, but if I calculate the coordination number of j, then the code has a problem. I need to first communicate with another processor to get its neighbor list information, then calculate the coordination number of j.
Am I right? I have to calculate the force between atom i and j through the coordination number of i and j. So I think I need to calculate the coordination number of j. I will study the code pair_eam.cpp. Thank you very much.
Lisa
2011/5/19 Steve Plimpton <[email protected]>
Thank you. I realize that I was thinking on the wrong way. I am not sure if
I understand it right now.
On one processor, the atom i has a neighbor j, which is on another
processor. This processor can get the position information of the atom j,
but not the neighbor list information. So if I just calculate the diatance
between i and j and calculate the force, it can work thoughly, but if I
calculate the coordination number of j, then the code has a problem. I need
to first communicate with another processor to get its neighbor list
information, then calculate the coordination number of j.
it is not as simple as that.
the neighbor list indices for each processor a different
as they point to a local list. the only way to identify atoms
are through their global tag value, and that can be very
time consuming.
i think the easier way is to request the proper neighbor list.
lammps can produce different types. there are half and full
neighbor lists. a half neighbor list has every i,j pair for local
atoms only once, and that is convenient for pairwise additive
potentials or contributions, since newton's third law requires
F_ij to be equal to -F_ji. in a full neighbor list both pairs are
contained. this may be more convenient in your case to get
the coordination number. the other important option is with
respect to the ghost particles (coordinates "owned" by other
processors). for ghost atoms, you only have interactions
with i being a local atom and j being a ghost. now, with the
flag newton on, you apply newtons 3rd law to these forces,
too, and then communicated the forces to the appropriate
neighboring processors. for your case it would be easier
to require that to be turned off. this way you will have for
every i particle _all_ particles in the j list that are within
the cutoff, regardless of whether they are local or remote.
if i am not totally mistaken, you should be able to avoid
the communication to compute the coordination number.
rather than looking at the eam potential, you may want
to check out stillinger weber (sw) instead. it also shows you
how to traverse the neighborlists for the force calculation
in a pair wise additive style and avoid the redundant force
calculation.
cheers,
axel.
Ok, thank you very much. I will try it.
Lisa
2011/5/19 Axel Kohlmeyer <[email protected]>
You're going down the wrong path. You cannot easily compute the coord #
of a ghost atom even if the neighboring proc gives you neighbor
info, since the central proc doesn't know about those neighbors.
You're better off letting each proc compute the coord # of the
atoms it owns, and then communicate the coord #. That is
what pair eam does with the per-atom density. Just replace
density with coord # and it's the same thing.
Steve
Thank you very much. I did as what you suggest and there is no segmentation fault wrong. I will do further simulations to test the code. I have one confusion. What is the difference between
comm->reverse_comm_pair(this) and comm->forward_comm_pair(this)?
I just want to communicate with other procs and do not sum the coordination up.
Lisa
2011/5/20 Steve Plimpton <[email protected]>
Forward acquires ghost atom properties from other proc's owned
atoms. Reverse send my ghost atom values to other proc's
owned atoms, typically so it can sum the contributions.
Steve
Sorry to bother you again. I forgot one part of the force before and now I need add the derivative of coordination number on the r into force. The two body potential has the form
f2=coord[i]*coord[j]*f2w.
coord[i]=sum of the f(rij)
So when I calculate the force of pair ij.
for (i=0;i<inum;i++)
{
j=neighbor list of i
{…calculate the energy and force of ij pair.
{
k=neighbor list of i
force[i]+=-df(rik)/drik;
}
///this will be also done for j neighbor list llist
{
l=neighbor list of j
force[j]+=df(rjl)/drjl;
}
How can I do the j part? I can not find the fimilar part in pair_eam.
Lan
2011/5/21 Steve Plimpton <[email protected]>
If you are asking how to get a quantity stored
with the j neighbors of i, some of which are ghost
atoms, then pair_eam.cpp does exactly that with
the comm->forward/reverse calls in the middle
of the compute() method.
Steve