What's the purpose of shearupdate in a granular pair potential?

Hi all,

For example in hertz/history, I notice:

        if (shearupdate) {
          shear[0] += vtr1*dt;
          shear[1] += vtr2*dt;
          shear[2] += vtr3*dt;

Are all history-dependent pairwise quantities stored here? Also can these quantities be stored in shear vector irrespective of shearupdate?


The shearupdate flag is set to zero during the “setup()” phase of a run. This is necessary since you do not update positions but just recover/generate missing data like forces and torques from the initial configuration before you start the actual MD and propagation of particles. If you do a “run 0” command you only do the setup() part of a run without any propagation.

1 Like

Hi Axel,

I think I might have run into an issue pertaining to shearupdate. Specifically related to my code, I boosted the size of the shearvector from its default 3 to 12 (I have a lot of history dependent quantities), and commanded for the shearvector to be updated only if shearupdate==1.

What I noticed is that if my input script has 2 run commands, the shear vector at the end of the 1st run command is different to the shear vector of the 2nd run command as verified by dumping out the shear vector.

For example, if I have a run command for 20 timesteps, followed by another run command for 5 timesteps, I get the following shearvector quantities at the end of the 20th timestep

 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 
0.000213412582525 0.000202741953399 10052.153373580449625 1.079405422295767 
0.000213412582525 10052.153373580449625 1.117713565356370 0.000000000000000

and I get this shearvector at the beginning of the ‘setup phase’ of the second run command:

 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 
0.000213412582525 0.000202741953399 -10052.153373580449625 -1.079405422295767 
0.000213412582525 -10052.153373580449625 -1.117713565356370 0.000000000000000

I notice that some quantities are negative although the magnitude is the same, and I did not perform any operations in my pair_gran source file in between outputting both of these quantities. Would you happen to know why this sign change happens?

I am confident that the shear vector is being stored correctly for the first run command as it was verified with an independent script with hand calculations but it is just this transition between two run commands that is giving this issue.

Many thanks as always!

No. Without looking at the actual code and doing runs with instrumented binaries, valgrind and more (in other words: throwing the full debugging broadside at it), I can only speculate.

Hello Axel,

I understand. I was wondering if you would happen to know what this bit of code does as found in PairGranHookeHistory:: init_style():

  if (history) {
    irequest = neighbor->request(this,instance_me);
    neighbor->requests[irequest]->id = 1;
    neighbor->requests[irequest]->half = 0;
    neighbor->requests[irequest]->granhistory = 1;
    neighbor->requests[irequest]->dnum = 3;

My potential is a child of pair_gran_hooke_history, and I was wondering if the last line of dnum=3 has any bearing on the shear vector?

Thank you!

This requests a neighbor list. But this is outdated and not compatible with current LAMMPS.
Please see:

1 Like