[lammps-users] pressure is not consistent when continuing a simulation

Hi All,

I found the pressure from the last step in a run is not the same as the first step in the continuing run, however, all other thermodynamics quantities match. Since I have some rigid bodies in my system, I have tried to use the Lammps example (examples/rigid) to check what’s wrong. I found the same thing happens with that example file.

I have modified the input file as the following way:

  1. insert the following two line before “neigh_modify exclude group clump1 clump1”
    neighbor 1 bin
    neigh_modify every 1 delay 0 check no
    The lammps would complain about non-zero dangerous builds without this modification.

  2. At the end, I change the run steps to :

run 100
run 100

I expect I would get the same value for all thermodynamics quantities at the 100th time step from both run commands. How ever the pressure is not consistent. The following is the result from the log file:

Step Temp E_pair E_mol TotEng Press
0 115.29439 5235.9179 0 5272.2142 -2.7403788
50 14910.685 571.71558 0 5265.82 32.006171
100 16298.442 136.66184 0 5267.653 16.444229
Loop time of 0.0300425 on 8 procs for 100 steps with 81 atoms

Pair time () = 0.000577599 (1.9226) Neigh time () = 0.00348008 (11.5838)
Comm time () = 0.0232758 (77.4761) Outpt time () = 0.000445336 (1.48235)
Other time (%) = 0.00226372 (7.53507)

Nlocal: 10.125 ave 44 max 0 min
Histogram: 6 0 0 0 0 0 1 0 0 1
Nghost: 41.375 ave 51 max 36 min
Histogram: 2 0 2 0 2 1 0 0 0 1
Neighs: 108.75 ave 410 max 0 min
Histogram: 6 0 0 0 0 0 0 0 0 2

Total # of neighbors = 870
Ave neighs/atom = 10.7407
Neighbor list builds = 100
Dangerous builds = 0
run 100
Memory usage per processor = 1.5774 Mbytes
Step Temp E_pair E_mol TotEng Press
100 16298.442 136.66184 0 5267.653 20.042934
150 16682.606 17.490511 0 5269.4219 14.900344
200 16733.929 1.372872 0 5269.4617 14.569267

Does anyone have this experience?

All the best,

Dongsheng

I think this is simply because the contribution that fix rigid
makes to the virial (and hence to the pressure) is estimated
on the first step of a run. This is b/c in velocity Verlet,
the forces at the beginning of a run are not used for
an entire timestep but for 1/2 step. So fix rigid and fix shake
(for example) have to estimate the total contribution to
the virial, as 2x the 1/2 step contribution in this case. This
is not the same as how it would be done in a continuing
run, so the answers are a bit different.

Steve