[lammps-users] Question about neighbor list, PBC and distance dependence rule

Dear LAMMPS users,

I’m using LAMMPS (3 Mar 2020) with rigid, dipole, molecule and misc package. I have compiled LAMMPS with serial by now, not with mpi (eventually I would like to do it for mpi)

  1. Which flags do I have to use to find the full neighbor list when I use PBC conditions? As you can see in the code below I already use the request for the full list. I’m using in the LAMMPS script neighbor 0.0 bin and no other neighbors modify by now

  2. Also what should be the proper way to set the neighbor list cutoff? Would the command neighbor bin be the correct one for this purpose?

If you need LAMMPS script, Initial Snapshot or anything else, please let me know

I send you more details in case that you need it:

I’m trying to implement active matter but with local momentum conservation with a specific rule. Imagine that I have a system with particles of type i (input, let say 1) and one particle of type j (input, let say 2. This is actually a rigid spherical colloid made of particles that have one angle bond with the type 2 particle in the center and two particles in the edge of the colloid to define an axis). I want to apply a force F_a (input, around 300 in lj units) in the particle of type 2 (in the direction of the axis defined by the angle bond) and apply a force F_p (computed in the program with a specific rule) on each neighbor (type 1) of this particle such

sum over all neighbors of type 1 (F_p) = - F_a

This way sum(F)=F_a-F_a =0 all time (except numerical fluctuations) and local momentum will be conserved with this redistribution of the force. I find a similar post in the mail list but it doesn’t solve my problem

https://lammps.sandia.gov/threads/msg26076.html

I already implemented a code to do that but use my own list of neighbors “by hand” (with two loops over all atoms, computing the distance and with a cutoff… not optimal…) and works properly. The problem is that I have tried to do the same using LAMMPS’ neighbor list and something happens when the list has to use PBC, but with no PBC the code works. When the particle of type 2 is close to the boundary (this is, when the cutoff for the list is more than the distance of the particle to the boundary) the sum over all neighbors (F_p) != - F_a as you can see in the figure below (I show the modulus of the force here)

Screenshot from 2021-04-19 18-15-46.png

I have compared the codes for the same cutoff (3.0 by hand and I suppose that 3.0 if lammps takes the neighbor list cutoff from the pair_style + skin distance) and it looks like the neighbors in the boundary are not exactly the same, and that is why the force is not zero, that’s mean that I’m doing something wrong with the PBC conditions…

Just to be sure I also set

neigh_modify every 1 delay 0 check no

and the neighbor list given by lammps had more neighbors but still not all of which were present in the “by hand” neighbor list

Thank you very much

ciao
José

/////////////////////////////////////////////////////////////////////////////////
//////////////////// fix_code.cpp //////////////////////////////////////

/////////////////////////////////////////////////////////////////////////////////

#include “domain.h”
#include “math.h”
#include “stdlib.h”
#include “string.h”
#include “fix_puller.h”
#include “atom.h”
#include “force.h”
#include “neighbor.h”
#include “neigh_list.h”
#include “neigh_request.h”
#include “group.h”
#include “update.h”
#include “error.h”

using namespace LAMMPS_NS;
using namespace FixConst;

/* ---------------------------------------------------------------------- */

FixPuller::FixPuller(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) {

if (narg == 9){
type1 = force->numeric(FLERR, arg[3]);
type2 = force->numeric(FLERR, arg[4]);
type_angle = force->inumeric(FLERR, arg[5]);
cutoff = force->numeric(FLERR, arg[6]);
Pe = force->numeric(FLERR, arg[7]);
Rc = force->inumeric(FLERR, arg[8]);
A=1;
B=0;
alpha=1;
beta=1;
}
else if (narg == 13){
type1 = force->numeric(FLERR, arg[3]);
type2 = force->numeric(FLERR, arg[4]);
type_angle = force->inumeric(FLERR, arg[5]);
cutoff = force->numeric(FLERR, arg[6]);
Pe = force->numeric(FLERR, arg[7]);
Rc = force->inumeric(FLERR, arg[8]);
A = force->inumeric(FLERR, arg[9]);
B = force->inumeric(FLERR, arg[10]);
alpha = force->inumeric(FLERR, arg[11]);
beta = force->inumeric(FLERR, arg[12]);
}
else{
error->all(FLERR, “Illegal fix hy/puller command”);
}

force_reneighbor = 1;
next_reneighbor = -1;

}

/* ---------------------------------------------------------------------- */

int FixPuller::setmask(){
int mask = 0;
mask |= POST_FORCE;
return mask;
}

/* ---------------------------------------------------------------------- */

void FixPuller::init_list(int /id/, NeighList *ptr)
{
list = ptr;
}

/* ---------------------------------------------------------------------- */
void FixPuller::init() {

int irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->pair = 0;
neighbor->requests[irequest]->fix = 1;
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->full = 1;

}

/* ---------------------------------------------------------------------- */

void FixPuller::setup(int vflag)
{
if (strstr(update->integrate_style,“verlet”))
post_force(vflag);
}

/* ---------------------------------------------------------------------- */

void FixPuller::post_force(int vflag){

////////////// For the PBC in the colloid axis
double *boxlo = domain->boxlo;
double *boxhi = domain->boxhi;
double x_size = boxhi[0] - boxlo[0];
double y_size = boxhi[1] - boxlo[1];
double z_size = boxhi[2] - boxlo[2];
/////////////////////////////////////////////

int i,j,ii,jj,inum,jnum,itype,jtype;
double xtmp, ytmp, ztmp;
int *ilist, *jlist,*numneigh, **firstneigh;

int i1,i2,i3,n;
double delx,dely,delz;
double rsq,r;

double **x = atom->x;
double **f = atom->f;
double **mu = atom->mu; // Used as auxiliar variable
int *type = atom->type;
int atype;
int **anglelist = neighbor->anglelist;
int nanglelist = neighbor->nanglelist;
double Fa;
double Fa_x, Fa_y, Fa_z;
int nlocal = atom->nlocal;
int newton_bond = force->newton_bond;

int max_num_neigh=20000;
int list_neigh[max_num_neigh], counter; //Auxiliar list
double norma, weight_factor[max_num_neigh], prod_scal, cos_theta_2;

for (n = 0; n < nanglelist; n++) { //[FOR in nanglelist]
atype = anglelist[n][3];
if(atype == type_angle){ // [IF angle types]
i1 = anglelist[n][0];
i2 = anglelist[n][1];
i3 = anglelist[n][2];

delx = x[i3][0] - x[i1][0];
if (delx > x_size * 0.5) delx = delx - x_size;
if (delx <= -x_size * 0.5) delx = delx + x_size;
dely = x[i3][1] - x[i1][1];
if (dely > y_size * 0.5) dely = dely - y_size;
if (dely <= -y_size * 0.5) dely = dely + y_size;
delz = x[i3][2] - x[i1][2];
if (delz > z_size * 0.5) delz = delz - z_size;
if (delz <= -z_size * 0.5) delz = delz + z_size;

rsq = delxdelx + delydely + delz*delz;
r = sqrt(rsq);

Fa = Pe/1.;

Fa_x = Fadelx/r;
Fa_y = Fa
dely/r;
Fa_z = Fa*delz/r;

if (newton_bond || i2 < nlocal) {/// [IF bonds]
f[i2][0] += Fa_x;
f[i2][1] += Fa_y;
f[i2][2] += Fa_z;
mu[i2][0] = Fa_x;
mu[i2][1] = Fa_y;
mu[i2][2] = Fa_z;
}/// [IF bonds]

}/// [IF angle types]

} /// [FOR in nanglelist]

inum = list->inum; // # of I atoms neighbors are stored for
ilist = list->ilist; // local indices of I atoms
numneigh = list->numneigh; // # of J neighbors for each I atom
firstneigh = list->firstneigh;// ptr to 1st J int value of each I atom

for (ii = 0; ii < inum; ii++){ // [FOR ii]
i = ilist[ii];
itype = type[i];
if(itype == type1){ //[If typei and type1]
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
jlist = firstneigh[i];
jnum = numneigh[i];
norma = 0;
counter = 0;
for (jj = 0; jj < jnum; jj++) { // [FOR jj]
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];

if (jtype == type2){ // [If typej and type2]
delx = x[j][0]-xtmp;
dely = x[j][1]-ytmp;
delz = x[j][2]-ztmp;

Fa_x=mu[i][0];
Fa_y=mu[i][1];
Fa_z=mu[i][2];

rsq = delxdelx + delydely + delz*delz;
r=sqrt(rsq);
list_neigh[counter]=j;
/// Redistribution Rule

prod_scal = delxFa_x + delyFa_y + delzFa_z;
cos_theta_2 = sqrt(0.5
prod_scal/(r*Fa) + 0.5);
weight_factor[counter] = (cutoff - r)*cos_theta_2;
norma += weight_factor[counter];
counter++;
} // [If typej and type2]
} // [FOR jj]

for (jj = 0; jj < counter; jj++) { // [Redistribution Loop]
j=list_neigh[jj];
f[j][0] -= Fa_xweight_factor[jj]/norma;
f[j][1] -= Fa_y
weight_factor[jj]/norma;
f[j][2] -= Fa_z*weight_factor[jj]/norma;
} // [Redistribution Loop]

} //[If typei and type1]

}// [FOR ii]

}

Dear LAMMPS users,

I’m using LAMMPS (3 Mar 2020) with rigid, dipole, molecule and misc package. I have compiled LAMMPS with serial by now, not with mpi (eventually I would like to do it for mpi)

  1. Which flags do I have to use to find the full neighbor list when I use PBC conditions? As you can see in the code below I already use the request for the full list. I’m using in the LAMMPS script neighbor 0.0 bin and no other neighbors modify by now

your request is sufficient to obtain a full neighbor list from a fix. however, your choice of neighbor list skin is not a good idea, since it will require to rebuild the neighbor list in every MD step to be correct.

  1. Also what should be the proper way to set the neighbor list cutoff?

the neighborlist cutoff is by default the largest cutoff of all pair style interactions. this is most of the time a sufficient choice and has the benefit that it allows efficient reuse of existing neighbor lists. you can see in compute rdf and compute adf how to request a different cutoff in a neighbor list request.

Would the command neighbor bin be the correct one for this purpose?

no.

[…]

I have compared the codes for the same cutoff (3.0 by hand and I suppose that 3.0 if lammps takes the neighbor list cutoff from the pair_style + skin distance) and it looks like the neighbors in the boundary are not exactly the same, and that is why the force is not zero, that’s mean that I’m doing something wrong with the PBC conditions…

you should not need to worry about PBC, these are taken care of by LAMMPS as part of the domain decomposition. the neighbor list and bond/angle/dihedral lists will always contain the closest neighbors for a given set of atoms based on the atom that “owns” the interaction. please note that these may be different periodic copies for the non-owned atoms, if you have a bonded interaction straddling a subdomain boundary.

axel.

Dear Axel,

the neighborlist cutoff is by default the largest cutoff of all pair style interactions. this is most of the time a sufficient choice and has the benefit that it allows efficient reuse of existing neighbor lists. you can see in compute rdf and compute adf how to request a different cutoff in a neighbor list request.

Thank you, it was useful to find the way to use a new cutoff in the fix using the request (I will write it just for users that maybe need this help too)
neighbor->requests[irequest]->cut = 1;
neighbor->requests[irequest]->cutoff = cutoff_value;

you should not need to worry about PBC, these are taken care of by LAMMPS as part of the domain decomposition. the neighbor list and bond/angle/dihedral lists will always contain the closest neighbors for a given set of atoms based on the atom that “owns” the interaction. please note that these may be different periodic copies for the non-owned atoms, if you have a bonded interaction straddling a subdomain boundary.

The point is that when I run the simulations and check every time the number of neighbors for the central particle where I’m applying the force (particle of type 2), it should have 32 neighbors of type 1 (the particles that are part of the rigid colloid which center is particle of type 2) and X neighbors of type 3 (that in the model are dpd particles), but this is what I get (see table below). As you can see at some point when the cutoff is close enough to the boundary and the colloid star to go through PBC I lost some of my atoms of type 1 from the neighbor list (last column on the table) and this atoms should be always there even for PBC (because I have a rigid object and type 1 atoms have to be fixed at the same distance). It makes me feel that there is something that I miss in the explanations for the lists and the images of atoms

I don’t want to waste your time, so if you are not sure what is going on, could you please recommend any LAMMPS function that uses a full neighbor list and some distance-dependence rule to calculate? This way I will try to find by myself

Thank you very much,
ciao

José

/// Step Total force in system Total Number of Neighbors Neighbors type 3 Neighbors type 1
1 0.000000 567 535 32
2 0.000000 569 537 32

135 0.000000 539 507 32
136 0.000000 542 510 32
137 0.000000 543 511 32
138 0.011048 534 502 32
139 0.234134 534 502 32
140 0.283695 531 499 32
141 0.424807 529 497 32
142 0.729579 532 500 32
143 1.235670 537 505 32
144 1.908956 544 512 32
145 7.364590 542 510 32
146 8.486852 544 512 32
147 9.815999 546 514 32
148 13.618118 545 513 32
149 26.076434 543 511 32
150 29.629257 534 502 32
151 32.788468 535 503 32
152 54.075279 540 508 32
153 62.671367 545 513 32
154 67.085843 540 508 32
155 74.160670 539 507 32
156 81.173937 538 506 32
157 82.643093 537 505 32
158 86.631725 531 499 32
159 98.556941 525 493 32
160 106.829399 522 490 32
161 110.248625 517 485 32
162 122.019600 506 474 32
163 130.706923 502 470 32
164 138.250838 497 465 32
165 145.203924 496 464 32
166 152.010099 496 464 32
167 153.289178 485 454 31
168 161.283072 482 451 31
169 167.650213 475 444 31
170 173.163933 468 437 31
171 179.127490 458 427 31
172 184.639553 448 417 31
173 188.245499 437 406 31
174 195.317210 426 395 31
175 197.876553 421 392 29
176 203.353505 407 379 28
177 122.122855 418 391 27
178 114.694889 429 401 28
179 109.427010 442 412 30
180 103.481637 450 419 31
181 97.494409 455 424 31
182 91.727639 462 431 31
183 81.829542 468 437 31
184 77.957204 486 455 31
185 71.725189 494 463 31
186 67.480463 505 474 31
187 63.829778 510 479 31
188 59.705530 517 485 32
189 56.257255 524 492 32
190 51.956494 532 500 32
191 48.221088 532 500 32
192 40.713060 539 507 32
193 35.431268 542 510 32
194 28.723456 552 520 32
195 26.146749 562 530 32
196 19.648395 558 526 32
197 17.783220 560 528 32
198 15.845565 565 533 32
199 12.789800 570 538 32
200 8.293480 570 538 32
201 5.992944 566 534 32
202 4.926889 569 537 32
203 2.581642 573 541 32
204 2.236681 576 544 32
205 1.783662 576 544 32
206 1.545694 576 544 32
207 0.955972 569 537 32
208 0.210595 570 538 32
209 0.147619 571 539 32
210 0.089014 577 545 32
211 0.058276 570 538 32
212 0.027785 567 535 32
213 0.009079 565 533 32
214 0.008049 564 532 32
215 0.000000 559 527 32
216 0.000000 560 528 32
217 0.000000 565 533 32
218 0.000000 562 530 32
219 0.000000 571 539 32
220 0.000000 566 534 32
221 0.000000 567 535 32

300 0.000000 557 525 32

[…]

The point is that when I run the simulations and check every time the number of neighbors for the central particle where I’m applying the force (particle of type 2), it should have 32 neighbors of type 1 (the particles that are part of the rigid colloid which center is particle of type 2) and X neighbors of type 3 (that in the model are dpd particles), but this is what I get (see table below). As you can see at some point when the cutoff is close enough to the boundary and the colloid star to go through PBC I lost some of my atoms of type 1 from the neighbor list (last column on the table) and this atoms should be always there even for PBC (because I have a rigid object and type 1 atoms have to be fixed at the same distance). It makes me feel that there is something that I miss in the explanations for the lists and the images of atoms

I don’t want to waste your time, so if you are not sure what is going on, could you please recommend any LAMMPS function that uses a full neighbor list and some distance-dependence rule to calculate? This way I will try to find by myself

since the two computes I mentioned are the only cases setting a custom neighbor list cutoff, those are your only options for code to compare to.
since you didn’t share much information about the specifics of your setup and settings I can only guess. One possible explanation for the behavior you see is that your neighbor list cutoff is larger than the communication cutoff. The latter determines up to what distance ghost atoms are created (since this is overhead). The default for this is the neighbor list cutoff, i.e. the largest pairwise cutoff plus neighbor list skin.

Ultimately, there are 3 possible reasons why you won’t find a specific (local) index as a neighbor:

  • it is a ghost atom and the communication cutoff is too short
  • the neighbor is excluded from the neighbor list. that can be done explicitly or through the special_bond setting of exactly 0.0
  • the neighbor is flagged with extra bits because it is a special_bonds neighbor where the scaling factor is neither 0.0 nor 1.0. that is why you should be doing a bitwise and with NEIGHMASK.

from those 3 conditions only the first one matches your observation.

axel.