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