Athermal quasistatic shear deformation - Fire Algorithm

Divas · January 12, 2024, 5:45pm

Dear all
I am trying to implement Athermal quasistatic shear deformation on a 2d triclinic box with athermal disks. I am using LAMMPS version of july 2021. I am using a 2d periodic box and applying shear in xy plane by using fix deform command with delta keyword, after each strain step I am relaxing the system to the static state using FIRE algorithm. To prepare the data file for the simulation, I have created 5000 particles randomly and then relax the system using FIRE algorithm. During this process, when the system gets strained there is a increase in the shear stress value, but when the relaxation starts after the strain , the value of shear stress is different (getting reduced)at zero step , which should me same as the end of the strain.
Please help me to understand this and please suggest me corrections if I am wrong at any point. I am also attaching my input script and snippet of the log file for the reference.


units		si          #si units of parameters
atom_style	sphere
dimension	2
boundary p p p
newton		off
comm_modify     vel yes

read_data    data.stiffness_10
group 1 type 1
group 2 type 2
set type 1 mass 1
set type 2 mass 1 

neighbor  1.2 bin
neigh_modify	delay 0 every 1 check yes

pair_style granular
pair_coeff * * hooke 10.0 00.0 tangential linear_history 0.0 1.0 0.0  damping mass_velocity


compute peratom all stress/atom NULL pair
compute p all reduce sum c_peratom[1] c_peratom[2] c_peratom[3] c_peratom[4]
#compute mychunk all chunk/atom bin/2d y lower 0.0 y 5.0 units box
variable area equal (102*102)
variable press equal -(c_p[1]+c_p[2])/v_area 
variable sigmaxy equal -c_p[4]/v_area

fix 3 all enforce2d
timestep	0.001
variable i loop 3
label loop

#velocity all ramp vx 0.0 0.0001 y 0.0 80.0 sum yes units box
fix 1 all deform 1 xy delta 0.0001 remap v units box
#fix 100 ave/chunk 1 5 1 mychunk vx 


fix int all nve/sphere disc

thermo 1
thermo_style custom step ke v_press v_sigmaxy cellgamma
dump           50 all custom 1 *.stress_strain id type diameter x y  c_peratom[1] c_peratom[2] c_peratom[3] c_peratom[4]
dump_modify 50 sort id
#restart 1 restart.*.aqs_shear
run 1 pre yes post yes
unfix 1
unfix int
undump 50
thermo 1
min_style  fire
min_modify  tmax 1.0
dump           20 all custom 1 *.stress_relax id type diameter x y  c_peratom[1] c_peratom[2] c_peratom[3] c_peratom[4]
dump           30 all custom 100 *.shearforces id type diameter x y  fx fy tqz vx vy 
dump           40 all atom 100 shear_aqs.dump
dump_modify    20 sort id 
dump_modify    30 sort id 
#restart 1 restart.*.aqs_shear
minimize       0.0  1.0e-6  1 10
undump 20
undump 30
undump 40
next i
jump SELF loop

LAMMPS (30 Jul 2021)

units		si          #si units of parameters
atom_style	sphere
dimension	2
boundary p p p
newton		off
comm_modify     vel yes

read_data    data.stiffness_10
Reading data file ...
  triclinic box = (0.0000000 0.0000000 -0.50000000) to (80.000000 80.000000 0.50000000) with tilt (0.0000000 0.0000000 0.0000000)
  3 by 1 by 1 MPI processor grid
  reading atoms ...
  5000 atoms
  reading velocities ...
  5000 velocities
  read_data CPU = 0.052 seconds
group 1 type 1
2500 atoms in group 1
group 2 type 2
2500 atoms in group 2
set type 1 mass 1
Setting atom values ...
  2500 settings made for mass
set type 2 mass 1
Setting atom values ...
  2500 settings made for mass

neighbor  1.2 bin
neigh_modify	delay 0 every 1 check yes

pair_style granular
pair_coeff * * hooke 10.0 00.0 tangential linear_history 0.0 1.0 0.0  damping mass_velocity


compute peratom all stress/atom NULL pair
compute p all reduce sum c_peratom[1] c_peratom[2] c_peratom[3] c_peratom[4]
#compute mychunk all chunk/atom bin/2d y lower 0.0 y 5.0 units box
variable area equal (102*102)
variable press equal -(c_p[1]+c_p[2])/v_area
variable sigmaxy equal -c_p[4]/v_area

fix 3 all enforce2d
timestep	0.001
variable i loop 3
label loop

#velocity all ramp vx 0.0 0.0001 y 0.0 80.0 sum yes units box
fix 1 all deform 1 xy delta 0.0001 remap v units box
#fix 100 ave/chunk 1 5 1 mychunk vx


fix int all nve/sphere disc

thermo 1
thermo_style custom step ke v_press v_sigmaxy cellgamma
dump           50 all custom 1 *.stress_strain id type diameter x y  c_peratom[1] c_peratom[2] c_peratom[3] c_peratom[4]
dump_modify 50 sort id
#restart 1 restart.*.aqs_shear
run 1 pre yes post yes
Neighbor list info ...
  update every 1 steps, delay 0 steps, check yes
  max neighbors/atom: 2000, page size: 100000
  master list distance cutoff = 2.6
  ghost atom cutoff = 2.6
  binsize = 1.3, bins = 62 62 1
  1 neighbor lists, perpetual/occasional/extra = 1 0 0
  (1) pair granular, perpetual
      attributes: half, newton off, size, history
      pair build: half/size/bin/newtoff
      stencil: full/bin/2d
      bin: standard
Per MPI rank memory allocation (min/avg/max) = 12.38 | 12.38 | 12.38 Mbytes
Step KinEng v_press v_sigmaxy CellGamma 
       0 5.8814887e-15   0.42476895 0.00053117843           90 
       1 5.8824558e-15   0.42476895 0.00053117843    89.999928 
Loop time of 0.00505827 on 3 procs for 1 steps with 5000 atoms

91.7% CPU use with 3 MPI tasks x no OpenMP threads

MPI task timing breakdown:
Section |  min time  |  avg time  |  max time  |%varavg| %total
---------------------------------------------------------------
Pair    | 0.00048795 | 0.00049006 | 0.00049418 |   0.0 |  9.69
Neigh   | 0          | 0          | 0          |   0.0 |  0.00
Comm    | 4.2293e-05 | 4.3589e-05 | 4.429e-05  |   0.0 |  0.86
Output  | 0.0043866  | 0.0043904  | 0.0043971  |   0.0 | 86.80
Modify  | 8.3509e-05 | 9.2254e-05 | 9.7209e-05 |   0.0 |  1.82
Other   |            | 4.193e-05  |            |       |  0.83

Nlocal:        1666.67 ave        1672 max        1660 min
Histogram: 1 0 0 0 0 0 1 0 0 1
Nghost:        450.333 ave         455 max         441 min
Histogram: 1 0 0 0 0 0 0 0 0 2
Neighs:        11474.3 ave       11537 max       11424 min
Histogram: 1 0 0 1 0 0 0 0 0 1

Total # of neighbors = 34423
Ave neighs/atom = 6.8846000
Neighbor list builds = 0
Dangerous builds = 0
unfix 1
unfix int
undump 50
thermo 1
min_style  fire
min_modify  tmax 1.0
#dump           20 all custom 1 *.stress_relax id type diameter x y  c_peratom[1] c_peratom[2] c_peratom[3] c_peratom[4]
dump           30 all custom 100 *.shearforces id type diameter x y  fx fy tqz vx vy
#dump           40 all atom 100 shear_aqs.dump
#dump_modify    20 sort id
#dump_modify    30 sort id
#restart 1 restart.*.aqs_shear
minimize       0.0  1.0e-6  1 10
  Parameters for fire: 
    dmax delaystep dtgrow dtshrink alpha0 alphashrink tmax tmin    integrator halfstepback 
     0.1        20    1.1      0.5   0.25        0.99    1 0.02 eulerimplicit          yes 
Per MPI rank memory allocation (min/avg/max) = 10.79 | 10.79 | 10.79 Mbytes
Step KinEng v_press v_sigmaxy CellGamma 
       1            0   0.42476919  0.000528954    89.999928 
       2 4.0134557e-11   0.42476919  0.000528954    89.999928 
Loop time of 0.00203243 on 3 procs for 1 steps with 5000 atoms

77.8% CPU use with 3 MPI tasks x no OpenMP threads

Minimization stats:
  Stopping criterion = max iterations
  Energy initial, next-to-last, final = 
                     0                  0                  0
  Force two-norm initial, final = 0.0089593032 0.0089589762
  Force max component initial, final = 0.0023439137 0.0023438106
  Final line search alpha, max atom move = 0.0000000 0.0000000
  Iterations, force evaluations = 1 2

MPI task timing breakdown:
Section |  min time  |  avg time  |  max time  |%varavg| %total
---------------------------------------------------------------
Pair    | 0.0010019  | 0.0012206  | 0.0016093  |   0.8 | 60.05
Neigh   | 0          | 0          | 0          |   0.0 |  0.00
Comm    | 0.00010634 | 0.00010962 | 0.00011139 |   0.0 |  5.39
Output  | 0.00010017 | 0.00010594 | 0.00011717 |   0.0 |  5.21
Modify  | 1.8368e-05 | 2.1484e-05 | 2.4269e-05 |   0.0 |  1.06
Other   |            | 0.0005748  |            |       | 28.28

Nlocal:        1666.67 ave        1672 max        1660 min
Histogram: 1 0 0 0 0 0 1 0 0 1
Nghost:        450.333 ave         455 max         441 min
Histogram: 1 0 0 0 0 0 0 0 0 2
Neighs:        11474.3 ave       11537 max       11424 min
Histogram: 1 0 0 1 0 0 0 0 0 1

Total # of neighbors = 34423
Ave neighs/atom = 6.8846000
Neighbor list builds = 0
Dangerous builds = 0
#undump 20
undump 30
#undump 40
next i
jump SELF loop

srtee · January 13, 2024, 1:59am

Are you trying to replicate the methodology of a paper previously published by another group?

(If not – this is a complex set of manipulations you are trying to do to your system, and I do not know how to gauge your expectations of what the system should do in response. That is, it’s not obvious to me what the “right” or “wrong” answer should be, so it’s also not obvious what you might be doing “wrong” and what you should try changing.)

Divas · January 13, 2024, 4:18am

Dear sir
I am trying to validate the results of the following paper
“Transition rates for slip-avalanches in soft athermal disks under quasi-static simple shear deformations”

I am following the methodology used in the paper mentioned above

srtee · January 13, 2024, 9:45am

Thank you for supplying the paper reference. The changes you need to make are to change remap v in the fix deform instruction:

to remap x (which is also the default, so you can also just remove the remap keyword). You should also remove the fix nve/sphere integrator, as only the FIRE minimiser is responsible for moving the particles.

Please re-read the paper’s Numerical Methods section, especially the second paragraph, and convince yourself that these changes will bring your simulation in line with the described procedure.

Divas · January 13, 2024, 9:54am

Dear Sir
Thank you very much for your reply