Problem with using bonds with periodic boundary conditions

Dear all,

I am trying to implement an area constraint potential for cells (circular objects made of beads connected with springs and harmonic angle potentials) in LAMMPS. This potential is simply a quadratic potential that depends on the current area of the cell. To calculate the current area, I simply take the cross product between successive beads. (My version of LAMMPS is 30 July 2016).

I am trying to implement a fix for this calculation and the part of the code that does the area calculation is the following:

// zero out local per-mol areas
for (int i = 0; i < nmol; i++) {
currareaproc[i] = 0.0;
}

// compute current per-mol areas
for (int n = 0; n < nbondlist; n++) {

// get the succesive bonds of a molecule/cell
int i1 = bondlist[n][0];
int i2 = bondlist[n][1];
int type = bondlist[n][2];

// get the molecule index of the atom // to see which cell the atom is connected to
int index = molecule[i2]-1;
if (index < 0) continue;

// calculate the cross products locally over procs
currareaproc[index] += (x[i1][0]*x[i2][1] - x[i2][0]*x[i1][1])/2.0;

}

// sum over the procs and distribute to each proc
MPI_Allreduce(currareaproc,currarea,nmol,MPI_DOUBLE,MPI_SUM,world);
MPI_Barrier(world);

When I specifically don’t add a force based on this area calculation, everything is fine. However, when I look at the videos of a single cell or multiple cells while monitoring the area by printing the area values, I observe that the area of a cell becomes unreasonably large when it is crossing the periodic boundary. So any kind of force that depends on area blows up as cells are passing the boundaries.

I turned to interweb services for some help and it looks like a lot of people had similar problems with the periodic boundaries. I tried to apply the suggestions provided over pages and pages of Google. The main suggestion was to increase the communication cutoff, which I did by setting “comm_modify cutoff ${diameter_of_a_cell}”, but that didn’t help. I tried also increasing the neighbor skin distance, but it also didn’t work. It seems like somehow bonded atoms see themselves over the length of the box. I also confirmed this by setting the check in neigh_modify command. When I run over a single processor with this check on, then LAMMPS gives an error saying that bond extent is larger than 1/2 the box length. However, when I run over multiple processors, then I get a bond missing error.

I am looking for a solution for the last couple of days, I also consulted my colleagues but I am still out of ideas at this point. It would be great to get some suggestions.

Thank you and best regards,
Ozer

What is bondlist()? It looks your own
data structure. I don’t know if it
is parallel, or stores indices of ghost atoms.
If it stores indices of atoms that are on
opposite sides of a periodic box, then
you can’t compute a correct cross product
that way. You’d have to use a minimum
image convention to assure that x_i1 and x_i2
are coords close to each other.

Steve

Between successive beads of the circle, I have harmonic springs, which I implement by “bond_style harmonic”. I create a list of bonds and then give the list as an input to LAMMPS.

My idea is to use the same list that goes over to the bond_style in this fix that I am implementing. So I essentially copied and pasted code from bond_harmonic.cpp to achieve this. So a more complete version of the area calculation is the following:

void FixAreaCons::post_force(int vflag)
{
double **x = atom->x;
double **f = atom->f;
int **bondlist = neighbor->bondlist;
int nbondlist = neighbor->nbondlist;
tagint *molecule = atom->molecule;

// zero out local per-mol areas
for (int i = 0; i < nmol; i++) {
currareaproc[i] = 0.0;
}

// compute current per-mol areas
for (int n = 0; n < nbondlist; n++) {

// get the succesive bonds of a molecule/cell
int i1 = bondlist[n][0];
int i2 = bondlist[n][1];
int type = bondlist[n][2];

// get the molecule index of the atom // to see which cell the atom is connected to
int index = molecule[i2]-1;
if (index < 0) continue;

// calculate the cross products locally over procs
currareaproc[index] += (x[i1][0]*x[i2][1] - x[i2][0]*x[i1][1])/2.0;

}

// sum over the procs and distribute to each proc
MPI_Allreduce(currareaproc,
currarea,nmol,MPI_DOUBLE,MPI_SUM,world);
MPI_Barrier(world);

}

I made sure that the bonded atoms are the bonded atoms that I give as input by printing out the global tags of the atoms in succession and it looks correct. The first atom and the last atom of the circle are connected to each other and so on.

I can’t debug your code. I’m telling you
that if you are getting incorrect cross products,
it could well be b/c (in parallel or serial)
you are not using the coords of 2 atoms
that are close together. It’s up to you to
figure out why that is happening. Bond/harmonic
obviously works across PBC, so maybe you
are mis-using its data structures, or not
creating your data structure in the same way.
Either way, you have to debug.

Steve

Dear all,

I am trying to implement an area constraint potential for cells (circular
objects made of beads connected with springs and harmonic angle potentials)
in LAMMPS. This potential is simply a quadratic potential that depends on
the current area of the cell. To calculate the current area, I simply take
the cross product between successive beads. (My version of LAMMPS is 30 July
2016).

I am trying to implement a fix for this calculation and the part of the code
that does the area calculation is the following:

  // zero out local per-mol areas
  for (int i = 0; i < nmol; i++) {
    currareaproc[i] = 0.0;
  }

  // compute current per-mol areas
  for (int n = 0; n < nbondlist; n++) {

    // get the succesive bonds of a molecule/cell
    int i1 = bondlist[n][0];
    int i2 = bondlist[n][1];
    int type = bondlist[n][2];

    // get the molecule index of the atom // to see which cell the atom is
connected to
    int index = molecule[i2]-1;
    if (index < 0) continue;

    // calculate the cross products locally over procs
    currareaproc[index] += (x[i1][0]*x[i2][1] - x[i2][0]*x[i1][1])/2.0;

  }

  // sum over the procs and distribute to each proc
  MPI_Allreduce(currareaproc,currarea,nmol,MPI_DOUBLE,MPI_SUM,world);
  MPI_Barrier(world);

When I specifically don't add a force based on this area calculation,
everything is fine. However, when I look at the videos of a single cell or
multiple cells while monitoring the area by printing the area values, I
observe that the area of a cell becomes unreasonably large when it is
crossing the periodic boundary. So any kind of force that depends on area
blows up as cells are passing the boundaries.

I turned to interweb services for some help and it looks like a lot of
people had similar problems with the periodic boundaries. I tried to apply
the suggestions provided over pages and pages of Google. The main suggestion
was to increase the communication cutoff, which I did by setting
"comm_modify cutoff ${diameter_of_a_cell}", but that didn't help. I tried
also increasing the neighbor skin distance, but it also didn't work. It
seems like somehow bonded atoms see themselves over the length of the box. I
also confirmed this by setting the check in neigh_modify command. When I run
over a single processor with this check on, then LAMMPS gives an error
saying that bond extent is larger than 1/2 the box length. However, when I
run over multiple processors, then I get a bond missing error.

I am looking for a solution for the last couple of days, I also consulted my
colleagues but I am still out of ideas at this point. It would be great to
get some suggestions.

i don't think that walking the bondlist is such a great idea. for
serial execution, you can still make it work, when you make certain,
that you have consistent coordinates, e.g. by using unwrapped
positions or finding the minimum image positions with respect to the
center of your object. but it would still not be compatible with
newton force off and it will break with parallel calculations, as
parts of your object may be in different subdomains.

i would suggest something like the following:

1) use the molecule id to identify each object. require that atoms in
the objects have globally unique and consecutive atom ids (as per the
usual LAMMPS conventions, atoms with the molecule id should be
exempted from this test).

2) use compute chunk/atom molecule to get a chunk data structure and
use it in a new area/chunk compute

3) take compute com/chunk and extend it for your purposes to become
area/chunk. in all generality, you need to do two passes in your new
compute for each chunk you do.
   a) compute the center of mass of each object. this is your
reference coordinate. if you want your compute to be general, you
should also count the number of atoms per molecule id.
   b) loop over the local atoms with the same chunk id. for each atom
get its local index i and its atom id via atom->tag[i]. then get the
local id of the next atom via j = atom->map(atom->tag[i]+1) and check
if it has the same molecule id as atom i. (there are some
complications to be considered here, but more about this later). get
the coordinates for each of the two atoms and convert them to minimum
image coordinates relative to the center of mass. now you can compute
your triangular area and add it to the per chunk storage and finally
sum it all up.
   c) special cases to consider. if atom id j is >= nlocal, then the
slice is built from two atoms that are owned on two different
processors. as noted above, that can be a problem when using bonds, as
the bond could be listed twice. however, if you stick with adding the
area for all i < nlocal cases (which should be the case), it should
work out. the other special case to worry about is the triangle betwee
the first and the last atom in the chain (or rather the one from the
last to the first). this would be the case when atom j has a different
molecule id than i, in that case, you should take the previously
computed atoms_per_chunk value and use the coordinates of j =
atom->map(atom->tag[i]+1-atoms_per_chunk). if any of the atom->map()
calls returns a negative number, then something is amiss, most likely
your bond is too long and you need to check your model and/or increase
the communication cutoff accordingly.

axel.

Thank you very much for your help. Following the idea of Axel Kohlmeyer, I made it work now. As a reference to anyone who might read this in the future, I want to add two more points on the method proposed by Axel Kohlmeyer:

  • There is one more complication in connecting the first and the last atoms, and that is with the case when the atom is the last atom in the whole simulation (when tag[i] = natoms), as it leads to problems in getting the id of next atom via j = atom->map(tag[i]+1-atoms_per_mol). But this could be easily circumvented by checking if the tag[i] == natoms, and then getting the next atom local id.

  • The second problem for me was the fact that if j is a ghost atom, then atom->map() doesn’t always give me the closest image. Using domain->closest_image(i,j) solved this problem.

  • The rest is done as how Axel Kohlmeyer proposed it, except that instead of using chunks, I used molecule ids, but that’s about it.

Thanks again and best regards,
Özer

Thank you very much for your help. Following the idea of Axel Kohlmeyer, I
made it work now. As a reference to anyone who might read this in the
future, I want to add two more points on the method proposed by Axel
Kohlmeyer:

- There is one more complication in connecting the first and the last atoms,
and that is with the case when the atom is the last atom in the whole
simulation (when tag[i] = natoms), as it leads to problems in getting the id
of next atom via j = atom->map(tag[i]+1-atoms_per_mol). But this could be
easily circumvented by checking if the tag[i] == natoms, and then getting
the next atom local id.

sorry, but this doesn't make any sense. i suspect, that there is a bug
in your code, that makes it work by chance, but not by construction.
please note, that the order of (local) atoms can change during
neighborlist rebuilds.
if tag[i] == natoms, then nexttag = tag[i]+-atoms_per_mol is
*definitely* a meaningful number and you can call atom->map(nexttag)
and get the first atom of the molecule.
i also don't understand how your solution *can* work correctly. it
doesn't make sense (at least not in the way you describe it).

- The second problem for me was the fact that if j is a ghost atom, then
atom->map() doesn’t always give me the closest image. Using
domain->closest_image(i,j) solved this problem.

i mentioned that explicitly at the very beginning of my response.

axel.

  • quote * loop over the local atoms with the same chunk id. for each atom
    get its local index i and its atom id via atom->tag[i]. then get the
    local id of the next atom via j = atom->map(atom->tag[i]+1) and check
    if it has the same molecule id as atom i. * end quote *

Let’s assume my atom id is i (given that i < nlocal). When I call j = atom->map(atom->tag[i]+1), shouldn’t that create problems because tag[i]+1 simply doesn’t exist?

I invoke j = atom->map(tag[i]+1-atoms_per_mol[mol_index_i]) only when mol_index_i != mol_index_j and when tag[i] == natoms. (natoms in the global sense, not nlocal+nghost, but the total number of atoms in the simulation)

1 Like

* quote * loop over the local atoms with the same chunk id. for each atom
get its local index i and its atom id via atom->tag[i]. then get the
local id of the next atom via j = atom->map(atom->tag[i]+1) and check
if it has the same molecule id as atom i. * end quote *

Let’s assume my atom id is i (given that i < nlocal). When I call j =
atom->map(atom->tag[i]+1), shouldn’t that create problems because tag[i]+1
simply doesn’t exist?

you are making a good point here, i was assuming that atom->map()
should return -1 in this case, which is something you need to check
for in any case, as it may also happen in situations where your
molecule is scattered across two or more subdomains, the communication
cutoff is too short, and thus the j atom isn't present among the ghost
atoms.
however, when checking the implementation of atom->map() it turns out,
that this behavior requires the use of a hash table for atom maps,
something that needs to be explicitly selected. the default is to use
an array, which will result in an off-by-one error and return some
undefined number instead of -1.

axel.

That’s exactly the behavior I observed. It looked like the last atom in the simulation was having random neighbor atoms.