Dear LAMMPS users,
There seems to be a bug in the single functions of pair_gran_hooke_history.cpp and pair_gran_hertz_history.cpp which causes the shear forces to be written out incorrectly. In both functions, the following lines of code are present:
int *touch = list->listgranhistory->firstneigh[i];
double *allshear = list->listgranhistory->firstdouble[i];
for (int jj = 0; jj < jnum; jj++) {
neighprev++;
if (neighprev >= jnum) neighprev = 0;
if (touch[neighprev] == j) break;
}
double shear = &allshear[3neighprev];
The problem is in the loop used to set the value of ‘neighprev’. For particle i, this code cycles through its contacts (of which there are usually jnum). Once particle j has been located, we wish to break out of the loop having found the appropriate value of neighprev. neighprev is then used to find the position in allshear corresponding to the i-j contact.
The condition used to break out of the loop compares j with ‘touch[neighprev]’. j, being a particle id, could have any value (5, 28, 723,…) but entries in ‘touch’ must be either 0 (indicating no contact) or 1. Hence, this break condition is almost never invoked. This is a problem because neighprev almost always has the same value after this loop as it had before the loop: it is incremented jnum times within the loop and reset to zero once. Therefore the same three shear values are always being selected from allshear regardless of j.
I think that this can be fixed very easily by replacing the current break condition with the following:
if (jlist[neighprev] == j) break;
Curiously the jlist pointer is already set in the single function at the moment although it isn’t used anywhere, so perhaps the original intention was to use jlist in this manner.
Regards,
Kevin Hanley
Post-doctoral researcher
Dept. of Civil and Environmental Engineering
Imperial College London, SW7 2AZ