I am currently running non-equilibrium molecular dynamics (NEMD) compression simulations and have encountered an inconsistency with the fix deform erate command that I would like to clarify, especially regarding modern LAMMPS versions. I am using the 2021 LAMMPS version. The reason? I must use an non-updated force field only available in that version.
Simulation Details:
fix 1 all deform 5000 x erate -0.0025 remap x units lattice
Timestep: Δt=0.0002 ps
Purpose: Uniaxial compression (shrinkage) along the x-axis.
The LAMMPS documentation explicitly states that the erate (engineering strain rate) style changes the box dimension linearly according to the Engineering Strain formula based on the original length (L0):
L(t)=L0(1+erate⋅dt)
However, when tracking the box length (L) over a continuous run, I observe that the simulation is calculating the new length incrementally based on the current length (Lt) at each step:
L(t+Δt)=L(t)(1+erate⋅dt)
This iterative approach is mathematically equivalent to the exponential True Strain model over long periods:
L(t)=L0⋅exp(erate⋅t)
This non-linear shrinkage means the command effectively treats the input erate value as a True Strain Rate (trate), causing the box to shrink more slowly than predicted by the linear Engineering Strain formula.
Is this exponential behavior (using Lt for the incremental calculation) the intended and permanent implementation for fix deform erate across modern LAMMPS versions (including 2021 and later)?
Has the contradictory linear formula in the documentation been updated in the latest stable releases to reflect the code’s exponential/True Strain behavior?
Any insight from developers or experienced users would be greatly appreciated. Thank you!
I suggest you create a small and fast to run test input deck with a setup and force field that can be used with both your 2021 version and the latest LAMMPS release and evaluate yourself if the behavior of fix deform is the the same.
The intended behavior is (obviously) what is described in the documentation.
If you believe that you have found an inconsistency, you need to provide proof in the form of a proper minimal reproducer input deck using only features available in the current version of LAMMPS and then explain how exactly you determine the inconsistency.
Please note that for any feature that has been around for a while and is used as much as fix deform, any gross inconsistency as you are claiming it has, is most likely a problem of how you determine or observe the outcome and very unlikely (but not completely impossible) an issue in the implementation. At least that is my observation from many of such reports over the years.
I believe that, as you said, this is an implementation issue. I have noticed that, even although the “fix deform” was effectively configured, for each “run command” set after the fix deform, a new L0 is assigned, and, therefore, a new linear deformation will be execute.
This is documented behavior and consistent with the description of the fix deform and run commands. If you want to not have the fix being re-initialized at every run command, you must use the “start” keyword of the run command. See below for a minimal example:
units metal
region box block 0 10 0 10 0 10
create_box 1 box
mass 1 1.0
pair_style zero 1.0
pair_coeff 1 1
timestep 0.001
variable erate index 10.0
fix 1 all deform 1 x erate ${erate}
variable lx0 equal $(lx)
variable lx equal v_lx0+v_lx0*dt*v_erate*step
thermo_style custom step lx v_lx xlo xhi ly lz
thermo 1
run 100 start 0
run 100 start 0
run 100 start 0
Another option is to use the vel command and calculate what the velocity should be to match the initial erate. Then you don’t have to worry about fix reinitialization.