Simple shear in the athermal, quasistatic limit

Hello, Lammps junta,
I am trying to implement athermal simple shear in a 2D triclinic box. I am attaching the minimal script that I am using to do so. However, I am skeptical about the correctness of the implementation.
The idea is such that I created a 2D system composed of N atoms which are equilibrated at T=2 and then cooled at some glass rate to T=0. Subsequently, the system is sheared by tilting xy using fix deform and the energy minimization was performed to get the new glass. But, this last energy minimization was not successful as the stopping criterion is throwing inesearch alpha is zero. The script is:

dimension       2
units           lj
atom_style      atomic
boundary        p p p
pair_style      lj/cut 2.5
pair_modify 	shift yes  #/truncated and shifted only potential
pair_coeff 	1 1 1.0 1.0 2.5
pair_coeff 	2 2 0.5 0.88 2.2
pair_coeff 	1 2 1.5 0.8 2.0
mass 		* 1.0
velocity        all create 2.0 32451 dist gaussian 	
neighbor        0.3 bin
neigh_modify    delay 0 every 1 check yes
comm_modify 	vel yes
fix             2 all enforce2d
min_style       sd
minimize        0.0e-12 1.0e3 1000000 1000000
timestep        0.01

fix             1 all nvt temp 2.0 2.0 $(100*dt)
run             10000
unfix 1

fix             1 all nvt temp 2.0 0.00001 $(100*dt)
run		30000
unfix 1
unfix 2
undump 5

change_box 	all triclinic
fix	 	1 all deform 1 xy delta 2.0  units box
fix      	2 all nvt temp 0.00001 0.00001 $(100*dt)
run 		1

#/after cooling, run minimization using CG algo
min_style       cg
minimize        1.0e-12 0.0 1000000 1000000

Any suggestions will really help me with this implementation. So looking forward to your comments and suggestions.

Hi @syed,

I don’t quite understand what you mean when you say

What do you mean by not successful? If alpha is zero, then it means that you’re atoms do not move and the energy can’t be reduced more by the algorithm.


  • without any information on your system it is quite hard to be more helpful
  • Your script as you give it will not run, the undump 5 commands does not have a dump previously defined.

Hey @Germain thanks for the reply. And sorry for that undump command. I deleted some commands from my original script to make this minimal.
In one of my projects, I am studying mechanical yield in amorphous materials like glasses. I use athermal quasistatic shearing (AQS protocol) on a glassy state obtained at T=0, wherein the system is subject to small shear steps and then allowed, after each step, to find a new mechanically stable minimum of its potential energy.
To achieve this in LAMMPS, I cooled down my system at a rate to form a glass. Then subjected the glass to simple shear by tilting the xy of the triclinic box and then running the CG algorithm to reach the lowest PE minima. But, the energy minimization is leading me to that lowest possible energy state. So, my worry is whether this AQS protocol is possible and whether am I implementing it correctly. How is the PBC implemented and why I am not able to reach that lowest energy state of the sheared system?

Again I don’t quite understand what you want to achieve. When you say

The CG algorithm leads you to the nearest local potential energy minimum. Is it not what you expect? How can you tell that you’re not reaching the lowest energy state? This questions are very close to the one I already asked.

Apart from that without any reference from the literature, it is hard to tell you if you’re doing anything wrong with regard to the physics or the implementation. If you have any paper in mind or previous results, the best way to tell that you are doing things correctly is to reproduce the results.

On one hand you are concerned about accuracy, but then you use cutoffs that produce quite inaccurate results. Please note that shifting the energy to zero at the cutoff does not remove the discontinuity in the force, which is a particular problem for advanced minimization algorithms as its extrapolation schemes cannot anticipate those.

Why this command? I see nothing in your input that requires it.

This is overly complex. Why not tilt the box directly with the change_box command?

It is a very bad habit to post scripts in a forum when asking for help without testing them.

The best way to check this is to try an reproduce published results.

What does the PBC implementation have to do with this?

From which do you infer this?

I was shearing the system using fix deform but now I understand that I can use change_box instead.

By PBC I mean how LAMMPS will work out the PBC during minimization or any dynamics for this sheared simulation box or there will be no movement across the boundaries at all. @akohlmey @Germain

You are not answering my question. Your implication that LAMMPS won’t handle PBC correctly because it is a tilted box is rather offensive, as it implies that the people programming LAMMPS are careless and don’t properly test what they are doing (yes, there are bugs occasionally, but you are implying that the basic functionality is not working), and it also implies that many LAMMPS users have not noticed it for many years (which is equally offensive).

I cannot explain the details of how LAMMPS applies PBC here, but it is explained in detail in the recently published LAMMPS paper and also in the manual at: 4.4.1. Partitioning — LAMMPS documentation and 4.4.2. Communication — LAMMPS documentation