reverse ilist

Dear lammps-users,

Is there a way to get the index value ii from a ilist value i? Is there
any drawback to making a reverse ilist array? (e.g. reverseilist[i]=ii;)

Thanks for the advice,
John

Dear lammps-users,

Is there a way to get the index value ii from a ilist value i? Is there
any drawback to making a reverse ilist array? (e.g. reverseilist[i]=ii;)

What are you trying to achieve? I suspect there may be a better way.

Axel

Hi Axel,

I have a problem where I am trying to count up the number of 3 coordinated
atoms as part of a pair potential. I want to do this on the fly when I am
finding local neighbors. I am using a single array of size nmax to store
the count of 3 coordinated atoms. Later on I want to call back this
information with threecoord[i] or threecoord[j] for example. I am getting
a memory leak when I do this.

    threecoord[i]=0.0;
    for (jj = 0; jj < jnum; jj++) {
      j = jlist[jj];
      j &= NEIGHMASK;
      jtag = tag[j];
      jtype = map[type[j]];
      iparam_ij = elem2param[itype][jtype];

      dxij = x[j][0]-xtmp;
      dyij = x[j][1]-ytmp;
      dzij = x[j][2]-ztmp;
      rijsq = dxij*dxij + dyij*dyij + dzij*dzij;
      rij = sqrt(rijsq);

      if (rijsq > params[iparam_ij].lcutsq) continue;
        neighptrj[nj++] = j;
        threecoord[i] += fc(iparam_ij,rij);
    }
    localneigh[i] = neighptrj;
    numlocal[i] = nj;
    npntj += nj;

Later on I call this back for a conditional operation.
    spisum = threecoord[i];
    spjsum = threecoord[j];
    if (spisum > 2.0 && spisum < 4.0 && spjsum > 2.0 && spjsum < 4.0) {

And that's where I am getting a leak. Do I have to use a two dimensional
array to do this? For example do I have to do a threecoord[i][j] like we
do for j=jlist[jj] where jlist=firstneigh[i]? Or is there someway to
connect j back to the ilist, such as reverseilist[j]=jj of ilist?

Thanks for the advice,
John

Hi Axel,

I have a problem where I am trying to count up the number of 3 coordinated
atoms as part of a pair potential. I want to do this on the fly when I am
finding local neighbors. I am using a single array of size nmax to store
the count of 3 coordinated atoms. Later on I want to call back this
information with threecoord[i] or threecoord[j] for example. I am getting
a memory leak when I do this.

why do you make that list double precision?
also, it is not really clear what you want.
do you want a flag that tells "this atom has 3 neighbors within XX"?
or do you want a counter "this atoms has X neighbors within XX"?
or do you want a list "these are the 3 neighbors of atom I within XX"?

what happens, if atom I has more than 3 neighbors?
does this count as 3, too? ...and then do you look at
the 3 closest? or all permutations of 3? or ???

i am completely lost on the "backreference" thing. this doesn't
make any sense at all. once you count, you lose the information
about identity. you obviously have to do two passes on the
data. same as the other 3-body potentials, you can just do a
special "derived" neighborlist for atoms that are within the
"reduced"(?) cutoff and then process that to count and process
it again to apply the resulting forces.

i think you need to be more clear about these things before
somebody can really help you. it almost looks like it would
be easier to get somebody to do the implementation for you
than trying to iterate your through the entire process.

threecoord[i]=0.0;
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtag = tag[j];
jtype = map[type[j]];
iparam_ij = elem2param[itype][jtype];

 dxij = x\[j\]\[0\]\-xtmp;
 dyij = x\[j\]\[1\]\-ytmp;
 dzij = x\[j\]\[2\]\-ztmp;
 rijsq = dxij\*dxij \+ dyij\*dyij \+ dzij\*dzij;
 rij = sqrt\(rijsq\);

 if \(rijsq &gt; params\[iparam\_ij\]\.lcutsq\) continue;
   neighptrj\[nj\+\+\] = j;
   threecoord\[i\] \+= fc\(iparam\_ij,rij\);

}
localneigh[i] = neighptrj;
numlocal[i] = nj;
npntj += nj;

Later on I call this back for a conditional operation.
spisum = threecoord[i];
spjsum = threecoord[j];
if (spisum > 2.0 && spisum < 4.0 && spjsum > 2.0 && spjsum < 4.0) {

And that's where I am getting a leak. Do I have to use a two dimensional
array to do this? For example do I have to do a threecoord[i][j] like we
do for j=jlist[jj] where jlist=firstneigh[i]? Or is there someway to
connect j back to the ilist, such as reverseilist[j]=jj of ilist?

i suspect you are just stabbing into the dark and hoping
to hit your target by accident. that won't only work in extremely
rare cases. i would suggest to stop thinking like a cognitive
person and breaking down what you need into small steps
that a computer can handle. there are different ways to store
"associations" and they have different costs and benefits.
you have to find out what exactly you need and what exactly
it is that certain code does before you start implementing
"complicated" things.

sorry. there ain't no escape from the blues.

axel.

Hi Russell,

One possible problem with the threecoord[] array is that it's not properly
allocated. Try "memory->create(threecoord,atom->nmax,"SOMENAME")"

Also another problem is that you are accessing threecoord[j], where j can
be a ghost image. So you have to add one more step after the top piece,
which is the forward communication "comm->forward_comm_pair" to pass
threecoord information from local to ghosts. Otherwise ghost j does not
have info on threecoord[].

Best,
Ray

Hi Russell,

One possible problem with the threecoord array is that it's not properly
allocated. Try "memory->create(threecoord,atom->nmax,"SOMENAME")"

Also another problem is that you are accessing threecoord[j], where j can
be a ghost image. So you have to add one more step after the top piece,
which is the forward communication "comm->forward_comm_pair" to pass
threecoord information from local to ghosts. Otherwise ghost j does not
have info on threecoord.

that is a *very* good point.

i would recommend to have a look at pair_eam.cpp
since what eam does has some similarity: your first
collect the contributions of the neighboring atoms to
the "electron density" at atom i and then need to
communicate/collect the missing pieces (part of it
depends on whether you have newton_pair enabled
or not) and then compute the embedding term
based on that. in a second round.

axel.