As I recall, my fix was more invasive because the real root issue was
that shear history information had to be resynchronized in the
post-force phase of the integration step. This is a consequence of
the interaction between the neighbor list associated with the
history-dependent potential. The problems I think arise because the
shear neighbor can update *during* the force evaluation - where
changes in particle positions can put particles into (or out of
contact), and mandate an updating of the contact list details in order
to stay correctly synchronized.
Steve - I know that Dr Christoph Kloss already contacted you about the
addition to fix_shear_history.cpp being commented out and hence
I'm not sure that adding the requirement to the WriteRestart::command
for FixShearHistory::pre_exchange to be run is sufficient. I've run
some tests this morning, and it seems to be the case that the shear
history is still not correct unless the neighbour lists are rebuilt
prior to writing the restart file.
I might possibly have one explanation. In the WriteRestart::command
function, Comm::setup, exchange and borders are run to transfer atoms
between processors. If the neighbour lists are not rebuilt, list->inum
and list->ilist remain the same, despite the numbers of atoms on each
processor potentially changing. This will be sorted out when the
neighbour lists are rebuilt by Verlet::setup. However,
FixShearHistory::pre_exchange is run before this, and it makes use of
both inum and ilist. Since these were not updated after the atoms were
transferred (and inum > 0 if the run is continued after the restart
file is written by write_restart), the shear history may therefore be
subtly changed. There may be other effects as well caused by not
rebuilding the neighbour list.
Hi Kevin, Peter -
A few comments.
First, I agree with Kevin that I messed up the patch last
week. A key line was commented out, so I reposted a 18Jun12
patch this AM with the correction.
I also agree with Kevin that there may be an additional issue
with the fix history pre_exchange() being called twice w/out
neighbor lists being updated in-between when doing an explicit
write_restart() call between 2 runs. Need to think about that.
I don't understand Peter's comment that that the per-atom
fix shear/history info needs to be updated every timestep.
The point of the two data structs:
(a) neighbor list shear history info (stored with neighbor list)
(b) fix shear/history info (stored internally in fix shear/history)
is that the shear history info only needs to be current in one
of them. And that it gets transferred from (a) to (b) before
atoms migrate to new procs (or a restart file is written).
I do agree with both of you that there are inconsistencies when
a granular history run is continued (run 100, run 100 instead of
run 200), or a run is restarted from a restart file, meaning
that the code is not producing identical trajectories
in the 2 cases.
But I thought thru this a bit more this AM, and am now semi-convinced
that it may be impossible to continue or restart a run and get
identical trajectories (statistically it should be fine, just not
The reason is basically b/c the granular potentials
are velocity dependent.
Assume the shear history info is correct when you continue/restart (I
agree there may be some lingering bugs in this). If you want to get
the same thermodynamic info on the last step of the previous run (say
step 100) and the first step of the new run (also step 100), then you
need the compute() function in pair gran/hooke/history to give
identical forces at 2 points in the calculation:
(1) the last step 100 of the first run
where force was computed in the middle of the timestep
(2) the setup() call of the 2nd run
But the velocities used by (1) and (2) will not be the same.
That's b/c velocity/omega are updated at the very end of timestep
100 in the first run, which happens between (1) and (2).
Note that this is an issue whether you continue a run
with a 2nd "run" command, or via a restart file.
I don't see any simple way around this. It may be just an unfortunate
fact-of-life for velocity-dependent potentials in a velocity-Verlet
scheme. Nor do I understand why Peter thinks this issue is "fixed" by
his updating of the shear history every timestep.
Comments? Suggestions? Am I missing something here?