Law pressure fluctuations.

Dear Lammps users and developers,

I want to make different pressure packages with 1000 spherical particles. I couldn’t get the desired pressure specially with the very low pressures like 0.00008 and 0.00004. After I reached the desired pressure value, again I run the package 1000 000 cycles and during this period again pressure starts to drop from the already set value. I have included my input file and log file herewith and hope somebody help me out to create stable pressure value with the desired value. Following input and output data related to 0.00004 (+/-10%) pressure and I couldn’t hold it when I run final 1000000 time step. As a remedy, I increased the 50 000 time step up to a 500 000 withing the loop - “loopc” and it solved little bit my problem but haven’t completely solve it yet. (Following input file has extra information and I am apologize from you about that.).

INPUT:

Jamming of 2D bidisperse packings: Mu10

variable dmu equal 10
variable xkkt equal 1.0
variable visc equal 1

Different Pressure values in new directories

variable a1 loop 1
label loopa1
variable b1 equal ${a1}+9

if “{b1} == 1" then "variable pdest equal 0.04" & "variable y string p0.04" & "variable cor string ../../0.95" & "variable bscale equal 1.018" if "{b1} == 2” then “variable pdest equal 0.02” &
“variable y string p0.02” &
“variable cor string …/…/p0.04” &
“variable bscale equal 1.018”
if “{b1} == 3" then "variable pdest equal 0.008" & "variable y string p0.008" & "variable cor string ../../p0.02" & "variable bscale equal 1.0115" if "{b1} == 4” then “variable pdest equal 0.004” &
“variable y string p0.004” &
“variable cor string …/…/p0.008” &
“variable bscale equal 1.0055”
if “{b1} == 5" then "variable pdest equal 0.002" & "variable y string p0.002" & "variable cor string ../../p0.004" & "variable bscale equal 1.003" if "{b1} == 6” then “variable pdest equal 0.0008” &
“variable y string p0.0008” &
“variable cor string …/…/p0.002” &
“variable bscale equal 1.001”
if “{b1} == 7" then "variable pdest equal 0.0004" & "variable y string p0.0004" & "variable cor string ../../p0.0008" & "variable bscale equal 1.00055" if "{b1} == 8” then “variable pdest equal 0.0002” &
“variable y string p0.0002” &
“variable cor string …/…/p0.0004” &
“variable bscale equal 1.0003”
if “{b1} == 9" then "variable pdest equal 0.00008" & "variable y string p0.00008" & "variable cor string ../../p0.0002" & "variable bscale equal 1.00016" if "{b1} == 10” then “variable pdest equal 0.00004” &
“variable y string p0.00004” &
“variable cor string …/…/p0.00008” &
“variable bscale equal 1.000041”

variable pstep equal ${pdest}/10.0
#shell mkdir $y
#shell cd $y

Make multiple configurations in new directories

variable a2 loop 10
label loopa2
print “Directory Mu0.001/{y}/Run{a2}”
shell mkdir Run${a2}
shell cd Run${a2}
log log.${a2}

newton off
communicate single vel yes
dimension 2

read_restart {cor}/Run{a2}/restart

neighbor 0.1 bin
neigh_modify check yes

pair_style gran/hooke/history 1.0 {xkkt} 1.0 NULL {dmu} ${visc}
pair_coeff * *

fix 1 all nve/sphere
fix 2 all enforce2d
fix 3 all viscous 0.2

timestep 0.01

compute 1 all erotate/sphere
compute 2 all property/local patom1 patom2
compute 3 all pair/local fx fy fz p1 p2 p3 force p4
thermo_style custom step atoms ke c_1 press vol
compute_modify thermo_temp dynamic yes
thermo 10000

velocity all zero linear

fix 4 all deform 1 x scale {bscale} y scale {bscale} z scale 1.0 &
remap none units box
run 1
unfix 4
reset_timestep 0
run 100000

Loop over minimizations/dynamics

variable c loop 400
label loopc
print “LOOP INDEX = $c”
variable p equal press
variable q equal p-{pdest}
variable pdiff equal sqrt($q*q) variable length equal vol^(1.0/2.0) print "Box Length: {length}"
variable bvol equal vol
if “{pdiff} < {pstep}” then “jump …/in.pressure-mu0 break”
if “q > {pstep}” then "variable rscale equal (1.0+0.05/{bvol})^(1.0/2.0)" variable pstep1 equal 1.5*{pstep}
if “q > {pstep1}” then "variable rscale equal (1.0+0.09/{bvol})^(1.0/2.0)" variable pstep1 equal 2.0*{pstep}
if “q > {pstep1}” then "variable rscale equal (1.0+0.12/{bvol})^(1.0/2.0)" variable pstep1 equal 2.5*{pstep}
if “q > {pstep1}” then "variable rscale equal (1.0+0.15/{bvol})^(1.0/2.0)" variable pstep1 equal 3.0*{pstep}
if “q > {pstep1}” then "variable rscale equal (1.0+0.2/{bvol})^(1.0/2.0)" variable pstep1 equal 3.5*{pstep}
if “q > {pstep1}” then "variable rscale equal (1.0+0.25/{bvol})^(1.0/2.0)" variable pstep1 equal 4.0*{pstep}
if “q > {pstep1}” then "variable rscale equal (1.0+0.35/{bvol})^(1.0/2.0)" variable pstep1 equal 4.5*{pstep}
if “q > {pstep1}” then "variable rscale equal (1.0+0.5/{bvol})^(1.0/2.0)" variable pstep1 equal 5.0*{pstep}
if “q > {pstep1}” then "variable rscale equal (1.0+0.8/{bvol})^(1.0/2.0)" variable pstep1 equal 5.5*{pstep}
if “q > {pstep1}” then "variable rscale equal (1.0+1.5/{bvol})^(1.0/2.0)" variable pstep1 equal 6.0*{pstep}
if “q > {pstep1}” then “variable rscale equal (1.0+2.5/${bvol})^(1.0/2.0)”
if “$q == 0.0” then “variable rscale equal 0.999”
if “q < -{pstep}” then "variable rscale equal (1.0-0.05/{bvol})^(1.0/2.0)" variable pstep1 equal 1.5*{pstep}
if “q < -{pstep1}” then "variable rscale equal (1.0-0.09/{bvol})^(1.0/2.0)" variable pstep1 equal 2.0*{pstep}
if “q < -{pstep1}” then "variable rscale equal (1.0-0.12/{bvol})^(1.0/2.0)" variable pstep1 equal 2.5*{pstep}
if “q < -{pstep1}” then "variable rscale equal (1.0-0.15/{bvol})^(1.0/2.0)" variable pstep1 equal 3.0*{pstep}
if “q < -{pstep1}” then "variable rscale equal (1.0-0.2/{bvol})^(1.0/2.0)" variable pstep1 equal 3.5*{pstep}
if “q < -{pstep1}” then "variable rscale equal (1.0-0.25/{bvol})^(1.0/2.0)" variable pstep1 equal 4.0*{pstep}
if “q > {pstep1}” then "variable rscale equal (1.0-0.35/{bvol})^(1.0/2.0)" variable pstep1 equal 4.5*{pstep}
if “q > {pstep1}” then "variable rscale equal (1.0-0.5/{bvol})^(1.0/2.0)" variable pstep1 equal 5.0*{pstep}
if “q > {pstep1}” then "variable rscale equal (1.0-0.8/{bvol})^(1.0/2.0)" variable pstep1 equal 5.5*{pstep}
if “q > {pstep1}” then "variable rscale equal (1.0-1.5/{bvol})^(1.0/2.0)" variable pstep1 equal 6.0*{pstep}
if “q > {pstep1}” then “variable rscale equal (1.0-2.5/${bvol})^(1.0/2.0)”
if “q == -{pdest}” then “variable rscale equal 0.999”
print "{rscale}" fix 4 all deform 1 x scale {rscale} y scale ${rscale} z scale 1.0 &
remap none units box
run 1
unfix 4
run 50000

next c
jump …/in.pressure-mu0 loopc
label break

End Loops

print “BREAK”
variable c delete
variable p equal press
variable q equal p-{pdest}
if “$q > 0.0” then “variable rscale equal 1.000001”
if “$q < 0.0” then “variable rscale equal 0.999999”
if "q == 0.0" then "variable rscale equal 0.9999999" fix 4 all deform 1 x scale {rscale} y scale ${rscale} z scale 1.0 &
remap none units box
run 1
unfix 4

reset_timestep 0
run 1

Dump out IMAGE

dump 1 all image 1000001 image.*.pnm diameter diameter &
axes yes 0.8 0.02 zoom 1.5
dump_modify 1 backcolor white pad 9

Dump out COORDINATE information

dump 2 all custom 1000001 dump.cor id type xs ys zs radius

Dump out FORCE information

dump 3 all local 1000001 dump.lst index c_2[1] c_2[2] &
c_3[1] c_3[2] c_3[3] &
c_3[4] c_3[5] c_3[6] &
c_3[7] c_3[8]

Dump out VELOCITY information

dump 4 all custom 1000001 dump.vel vx vy vz omegax omegay omegaz

run 1000000

write_data data.restart
write_restart restart

Loop over configurations

clear
shell cd …/
next a2
jump ./in.pressure-mu0 loopa2

Loop over pressure values

clear
#shell cd …/
next a1
jump ./in.pressure-mu0 loopa1

OUTPUT (LOG FILE)

run 1000000
Memory usage per processor = 8.78267 Mbytes
Step Atoms KinEng 1 Press Volume
1 1000 2.5870705e-12 8.2492659e-10 4.2148245e-05 1499.3169
10000 1000 1.943907e-12 8.24178e-10 4.151136e-05 1499.3169
20000 1000 1.4407515e-12 8.2420492e-10 4.0990862e-05 1499.3169
30000 1000 1.3703846e-12 8.2394905e-10 4.0525445e-05 1499.3169
40000 1000 2.0721774e-12 8.2446355e-10 3.9841855e-05 1499.3169
50000 1000 2.8454535e-12 8.2477524e-10 3.8954572e-05 1499.3169
60000 1000 1.9510726e-12 8.2433761e-10 3.8334498e-05 1499.3169
70000 1000 2.0330681e-12 8.2424004e-10 3.7580088e-05 1499.3169
80000 1000 1.8533629e-12 8.2408555e-10 3.6870917e-05 1499.3169
90000 1000 1.6590961e-12 8.2397834e-10 3.6124815e-05 1499.3169
100000 1000 1.6471704e-12 8.178316e-10 3.5492283e-05 1499.3169
110000 1000 1.9964254e-12 8.1828545e-10 3.4864475e-05 1499.3169
120000 1000 1.606612e-12 8.1787739e-10 3.4310546e-05 1499.3169
130000 1000 2.215521e-12 8.1848671e-10 3.3656182e-05 1499.3169
140000 1000 1.6897598e-12 8.1812737e-10 3.3091016e-05 1499.3169
150000 1000 1.0957467e-12 8.175119e-10 3.2664569e-05 1499.3169
160000 1000 1.0438994e-12 8.1751155e-10 3.2265455e-05 1499.3169
170000 1000 9.3292552e-13 8.1746584e-10 3.1877325e-05 1499.3169
180000 1000 9.529973e-13 8.175393e-10 3.1445112e-05 1499.3169
190000 1000 9.1685337e-13 8.1754338e-10 3.1041145e-05 1499.3169
200000 1000 5.2950313e-13 8.1722393e-10 3.0755626e-05 1499.3169
210000 1000 3.5597569e-13 8.1710583e-10 3.0622155e-05 1499.3169
220000 1000 2.9366042e-13 8.1709091e-10 3.0519242e-05 1499.3169
230000 1000 2.7760723e-13 8.1709259e-10 3.0410945e-05 1499.3169
240000 1000 2.8039894e-13 8.170996e-10 3.0296637e-05 1499.3169
250000 1000 1.6060704e-13 8.1702902e-10 3.0190368e-05 1499.3169
260000 1000 1.2488739e-13 8.1701435e-10 3.0111185e-05 1499.3169
270000 1000 1.5217197e-13 8.1702254e-10 3.0028791e-05 1499.3169
280000 1000 2.3585822e-13 8.1710818e-10 2.9929015e-05 1499.3169
290000 1000 1.6755e-13 8.1704468e-10 2.9851738e-05 1499.3169
300000 1000 1.7525368e-13 8.1701498e-10 2.9764878e-05 1499.3169
310000 1000 3.4640046e-13 8.1705725e-10 2.9642935e-05 1499.3169
320000 1000 2.8339923e-13 8.1704262e-10 2.9505876e-05 1499.3169
330000 1000 2.6375733e-13 8.1703543e-10 2.9383283e-05 1499.3169
340000 1000 2.8930907e-13 8.170523e-10 2.9270564e-05 1499.3169
350000 1000 2.5753287e-13 8.1704618e-10 2.9164519e-05 1499.3169
360000 1000 2.623974e-13 8.1705348e-10 2.9023209e-05 1499.3169
370000 1000 2.7142926e-13 8.1706517e-10 2.8878188e-05 1499.3169
380000 1000 2.7749747e-13 8.1707095e-10 2.8736328e-05 1499.3169
390000 1000 2.9115074e-13 8.1708441e-10 2.8587026e-05 1499.3169
400000 1000 4.5125493e-13 8.1735548e-10 2.8408583e-05 1499.3169
410000 1000 3.3622808e-13 8.1710791e-10 2.8196337e-05 1499.3169
420000 1000 4.5066075e-13 8.173388e-10 2.8015033e-05 1499.3169
430000 1000 6.0096294e-13 8.1729216e-10 2.7788513e-05 1499.3169
440000 1000 5.5894108e-13 8.1720966e-10 2.7488051e-05 1499.3169
450000 1000 2.367608e-13 8.1700655e-10 2.7287757e-05 1499.3169
460000 1000 3.5621116e-13 8.173905e-10 2.7128147e-05 1499.3169
470000 1000 1.5378261e-13 8.1696379e-10 2.6970624e-05 1499.3169
480000 1000 2.8492802e-13 8.1715323e-10 2.6829504e-05 1499.3169
490000 1000 2.0407529e-13 8.1707788e-10 2.6711335e-05 1499.3169
500000 1000 1.8882901e-13 8.1704556e-10 2.6575743e-05 1499.3169
510000 1000 1.9177393e-13 8.1704837e-10 2.6465415e-05 1499.3169
520000 1000 1.2265363e-13 8.1695114e-10 2.6370778e-05 1499.3169
530000 1000 2.4097861e-13 8.1711079e-10 2.6269131e-05 1499.3169
540000 1000 1.4348468e-13 8.169709e-10 2.6201273e-05 1499.3169
550000 1000 1.3191044e-13 8.1695576e-10 2.6126451e-05 1499.3169
560000 1000 2.0220744e-13 8.1696476e-10 2.6038674e-05 1499.3169
570000 1000 1.1215104e-13 8.1682794e-10 2.59816e-05 1499.3169
580000 1000 1.3824537e-13 8.1686436e-10 2.5913435e-05 1499.3169
590000 1000 8.3237114e-14 8.1681711e-10 2.5884399e-05 1499.3169
600000 1000 1.020082e-13 8.1683463e-10 2.5831161e-05 1499.3169
610000 1000 7.9375452e-14 8.168188e-10 2.5800201e-05 1499.3169
620000 1000 8.4808637e-14 8.1682344e-10 2.5747739e-05 1499.3169
630000 1000 9.7133398e-14 8.1684277e-10 2.5709901e-05 1499.3169
640000 1000 8.3460101e-14 8.1682248e-10 2.5661242e-05 1499.3169
650000 1000 1.3439869e-13 8.168757e-10 2.5617054e-05 1499.3169
660000 1000 8.4301545e-14 8.1681921e-10 2.5572085e-05 1499.3169
670000 1000 1.3398231e-13 8.1686451e-10 2.5522182e-05 1499.3169
680000 1000 8.2556047e-14 8.168159e-10 2.5482253e-05 1499.3169
690000 1000 1.3813486e-13 8.1690069e-10 2.5429986e-05 1499.3169
700000 1000 9.4409437e-14 8.1684022e-10 2.5393928e-05 1499.3169
710000 1000 8.4339368e-14 8.1682013e-10 2.5336107e-05 1499.3169
720000 1000 9.4711249e-14 8.1682786e-10 2.5279599e-05 1499.3169
730000 1000 1.0553971e-13 8.1685752e-10 2.5236075e-05 1499.3169
740000 1000 8.2245446e-14 8.1682003e-10 2.5185847e-05 1499.3169
750000 1000 9.5993806e-14 8.1682453e-10 2.513292e-05 1499.3169
760000 1000 9.2589773e-14 8.1361214e-10 2.5089763e-05 1499.3169
770000 1000 7.6098881e-14 8.1359339e-10 2.5039967e-05 1499.3169
780000 1000 8.4653283e-14 8.1359705e-10 2.4995269e-05 1499.3169
790000 1000 7.3128455e-14 8.1359105e-10 2.4968449e-05 1499.3169
800000 1000 6.6911868e-14 8.1358965e-10 2.4932158e-05 1499.3169
810000 1000 7.0202138e-14 8.1359441e-10 2.4911515e-05 1499.3169
820000 1000 5.5732443e-13 8.1394618e-10 2.4671967e-05 1499.3169
830000 1000 4.0068819e-13 8.1380758e-10 2.4463569e-05 1499.3169
840000 1000 3.8322241e-13 8.1381167e-10 2.4295049e-05 1499.3169
850000 1000 3.0260032e-13 8.1352682e-10 2.4157239e-05 1499.3169
860000 1000 4.3241904e-13 8.1359164e-10 2.3958588e-05 1499.3169
870000 1000 5.7446141e-13 8.1366938e-10 2.3698974e-05 1499.3169
880000 1000 4.2497331e-13 8.136032e-10 2.3424663e-05 1499.3169
890000 1000 3.507537e-13 8.1354143e-10 2.3208992e-05 1499.3169
900000 1000 2.9542445e-13 8.135048e-10 2.3012035e-05 1499.3169
910000 1000 2.2962245e-13 8.134803e-10 2.2793612e-05 1499.3169
920000 1000 4.303875e-13 8.1376202e-10 2.2485836e-05 1499.3169
930000 1000 3.3424117e-13 8.1361682e-10 2.2289041e-05 1499.3169
940000 1000 5.8942582e-13 8.1414245e-10 2.2114959e-05 1499.3169
950000 1000 2.3436585e-13 8.1347405e-10 2.1927854e-05 1499.3169
960000 1000 1.5694237e-13 8.1345691e-10 2.1762032e-05 1499.3169
970000 1000 1.9195992e-13 8.1342509e-10 2.1653073e-05 1499.3169
980000 1000 2.3163231e-13 8.1353393e-10 2.1551734e-05 1499.3169
990000 1000 7.8291667e-14 8.1239501e-10 2.1421499e-05 1499.3169
1000000 1000 7.4942584e-14 8.1240906e-10 2.1275339e-05 1499.3169
1000001 1000 7.5009833e-14 8.1240918e-10 2.1275321e-05 1499.3169
Loop time of 238.641 on 1 procs for 1000000 steps with 1000 atoms

Pair time () = 164.794 (69.0554) Neigh time () = 0 (0)
Comm time () = 6.03731 (2.52987) Outpt time () = 0.0813289 (0.03408)
Other time (%) = 67.7278 (28.3806)

Nlocal: 1000 ave 1000 max 1000 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 157 ave 157 max 157 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 2071 ave 2071 max 2071 min
Histogram: 1 0 0 0 0 0 0 0 0 0

Total # of neighbors = 2071
Ave neighs/atom = 2.071
Neighbor list builds = 0
Dangerous builds = 0

write_data data.restart
write_restart restart

Loop over configurations

clear
shell cd …/
next a2
jump …/in.pressure-mu1 loopa2
print “Directory Mu10/{y}/Run{a2}”
Directory Mu10/p0.00004/Run3
shell mkdir Run${a2}
shell mkdir Run3
shell cd Run${a2}
shell cd Run3
log log.${a2}
log log.3

No one is going to decipher a complicated
input script such as that, except you.

Explain what you are trying to do in words
and what you think isn’t working like you expect.

Note that graunular particles are typically very
hard, and they have dissipative forces in their
force field. So running dynamics and trying
to equilibrate to some target pressure is an
odd thing to do.

Steve