I think that there might be two bugs in LAMMPS for granular materials in which shear history is stored and restart files are also being used.
The first is that the shear history may not be written correctly to the restart file. The information which is to be written to the restart file is gathered in a buffer; for shear information, the pack_restart function in FixShearHistory is used for this. When called, the pack_restart function adds the data that are stored in the partner, npartner and shearpartner arrays to the buffer for writing. The shearpartner array, which contains the shear displacements, is updated when the pre_exchange function of FixShearHistory is called. When this pre_exchange function is run, the displacements for the current timestep (in the shear array) are copied to shearpartner. The problem is that this pre_exchange function is run only when the neighbour lists are rebuilt, and neighbour lists are not necessarily updated immediately before the restart file is written. Therefore, the shear information that is written could be out of date. I think that this can be fixed by changing the following code in WriteRestart::Command…
if (domain->triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
domain->reset_box();
comm->setup();
comm->exchange();
comm->borders();
if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
write(file);
delete [] file;
…to read as follows, i.e., forcing a neighbour list rebuild:
if (modify->n_pre_exchange) modify->pre_exchange();
if (domain->triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
domain->reset_box();
comm->setup();
if (neighbor->style) neighbor->setup_bins();
comm->exchange();
if (atom->sortfreq > 0) atom->sort();
comm->borders();
if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
if (modify->n_pre_neighbor) modify->pre_neighbor();
neighbor->build();
write(file);
delete [] file;
Perhaps not all of these lines are necessary, and it would be more efficient to invoke these additions only if shear history is stored so an ‘if’ condition to check for this could be included.
A second, and more serious, bug occurs when shear history is being restored from a restart file. The unpack_restart function in FixShearHistory is invoked to set npartner, shearpartner and partner correctly. The init function of the Neighbor class is run and the NeighList constructor is called for each list to be built. The Verlet::setup function is called next, before the neighbour lists are actually built (e.g., by Neighbor::Build calling Neighbor::granular_bin_no_newton). Verlet::setup runs the following line: modify->setup_pre_exchange();, which causes FixShearHistory::pre_exchange to be run. At this stage, list->inum is zero; this is set when the neighbour lists are built. Within the pre_exchange function, all entries in the npartner array are set to zero. The loop is intended to repopulate the npartner array, but because inum is zero, no iterations of the loop are run. Therefore, the npartner array remains filled with zeros. Finally, on neighbour list rebuild (e.g., Neighbor::granular_bin_no_newton), all of the touchptr and shearptr entries are zeroed because of the zeroed npartner entries, and the shear history is lost.
As a temporary repair, I made a small change in FixShearHistory::pre_exchange so that it reads as follows:
if (list->inum > 0)
for (i = 0; i < nlocal; i++) npartner[i] = 0;
This seems to work based on some tests I have done, and ensures that the shear information is preserved. I’m not fully certain of the implications though, e.g., what happens if a restart file is written and read on systems with differing numbers of processors.
I would be interested to hear any thoughts and comments.
Regards,
Kevin Hanley
Postdoctoral researcher
Dept. of Civil and Environmental Engineering
Imperial College London, SW7 2AZ