Shear Flow - Fix deform - SLLOD - PBC issue

Hello Everyone,

Firstly, I use the latest stable version of LAMMPS (02-Aug-23) on Ubuntu 22.04. I conduct a Kremer-Grest simulation of a polymer melt under shear flow and run on the following issue:

At some random timesteps, Rg jumps from an expected value of ~5 to a value of ~1e+3 sigma. Then it goes back to the expected value.
Upon visualization with wrapped coordinates, I find that at those timesteps the chain “escapes” the simulation box. As shown in the attached picture, the chains can appear in both “broken” and “unbroken” states. The former results in the unreasonable Rg values (when I visualize the unwrapped coordinates, I see again the broken chains).

Here is another example where I used as initial snapshot the LAMMPS benchmark (bench/data.chain). Lx should be much bigger, but I am showing this only for visualization purposes. Initially, the chain is wrapped, and then at the next timestep it is made “unbroken” and escapes the box. Once again, at that timestep, Rg diverges.

Further, this problem occurs when I apply higher shear rates. At low shear rates, close to the Newtonian plateau, no such issue is observed. The higher the shear rate, the problem happens at earlier times. The longer the Lx, the issue arises at later times (for the same shear rate).

All, the simulations run until complete, and I find that the steady-state viscosity and normal stresses are close to what is available in the literature. There is qualitative, and I guess semiquantititave agreement. However, I use Noser-Hoover thermostat while in the paper (Shear thinning behavior of linear polymer melts under shear flow via nonequilibrium molecular dynamics | The Journal of Chemical Physics | AIP Publishing) dpd thermostat is used, which might explain the differences.

Note that in the last 3 shear rates, the chains have become “broken”.

Some more info:

  1. Increased skin distance and communication cutoff (some of these are unnecessary; have tried simpler versions, same problem)
    neighbor 1.0 bin
    neigh_modify every 1 delay 0 check yes
    comm_modify cutoff 3.0 vel yes

  2. Made sure the system was properly equilibrated and Lx is very large. System dimensions are (8L, L, L) with L = 24.5sigma, density 0.85.

  3. To further verify that the initial snapshot is not the cause of the issue, I also used the configuration from LAMMPS benchamarks (bench/in.chain). Same issues were reproduced.

  4. Tried Nose-Hoover with 1 and 3 chains.

  5. Tried SLLOD and p-SLLOD.

Did not expect anything from the last two steps to happen; were more like desperate attempts to solve the issue. Maybe there is something obvious that I am missing in my simulation. I attach the simulation script, where I have isolated the main commands.

Lastly, what might be the issue is that the box size in the x direction is not long enough, given that the chain expands in the x direction and the pbc creates this weird issue. I can increase the size in the x direction even more, but Lx is already ~ 50Rg(eq.). Any other suggestions? Would the thermostat choice (Nose-Hoover) generate this issue, and dpd would not? I do not see why…

Thank you!

in.lammps (2.0 KB)

This is documented behavior and explained in the fix deform documentation. LAMMPS has required constraints for the tilt factors of tilted boxes. If you go over the constraint, the box is flipped and the atoms remapped accordingly including modification of image flags.
It also explains why this is needed, how you can turn it off, and how this can lead to problems, specifically at higher shearing rates and longer trajectories.

Hi Axel,
Thank you for your response.

I have carefully read the fix deform and fix nvt/sllod documentation. Here is part of it:

“If the flip value is set to no, the tilt will continue to change without flipping. Note that if you apply large deformations, this means the box shape can tilt dramatically LAMMPS will run less efficiently, due to the large volume of communication needed to acquire ghost atoms around a processor’s irregular-shaped subdomain. For extreme values of tilt, LAMMPS may also lose atoms and generate an error.”

Flip is by default set to yes; thus, I would not expect the issue to occur.

Perhaps I should use as initial configuration a snapshot from a lower shear rate simulation, where no such issue was found? However, this happens after millions of steps, so the initial velocities not are not causing the problem, and this change would not fix anything… (Here is why I proposed that Problem with sllod energies and averages)

Can you please explain to me what I am missing? How can I fix it?

Thank you!

I don’t see any use of the LAMMPS inbuilt computes to calculate the radius of gyration. So you are using your own post-processing script?

In which case, you need to validate that your post-processing analysis is in fact physically and mathematically correct. It seems you are getting incorrect values whenever a polymer’s boundary crossings change. You will have to work out how to correct the analysis.

This is a known and old aspect of molecular dynamics – how to “unwrap” chains correctly for mathematical analysis of chain properties. It applies even to standard, fixed-volume simulations, as long as they are periodic. Consider a 1nm-long polymer freely diffusing inside a 100nm-wide box, and suppose the radius of gyration is naively calculated using “within-box” coordinates. Then the Rg will be calculated as about 1nm when the polymer is far from any box edges but about 100nm (off the top of my head) when the polymer is crossing any box edges.

My logic would conclude the opposite. Have you made an experiment to see if there is any difference when you set flip to no?

Thank you for your responses!

I am using both LAMMPS built-in methods and my own code. Both reproduce the same issue.
Here is part of the script to calculate Rg

compute        clmol poly chunk/atom molecule compress yes
compute        clrg  poly gyration/chunk clmol tensor
variable       lrg1 equal ave(c_clrg[1])
variable       lrg2 equal ave(c_clrg[2])
variable       lrg3 equal ave(c_clrg[3])
variable       lrg4 equal ave(c_clrg[4])
variable       lrg5 equal ave(c_clrg[5])
variable       lrg6 equal ave(c_clrg[6])
variable       lrg equal v_lrg1+v_lrg2+v_lrg3
fix            flrg all ave/time ${filefreq} 1 ${filefreq} v_lrg v_lrg1 v_lrg2 v_lrg3 v_lrg4 v_lrg5 v_lrg6 file rg_poly.out mode scalar

I think I do not get unphysical values whenever a polymer is crossing the boundary. I am getting those values after millions of steps (when the polymer is crossing the boundary). So, what I show in the 2nd picture is at long times. The longer the box in the x-dimension the later the issue arises.

Unwrapping or wrapping coordinates is not the issue. I am exporting the trajectory in unwrapped (and wrapped), so the molecules should be whole, and no weird Rg values should appear.
The snapshot I show is also in wrapped coordinates, using VMD to wrap them.

Yes, I have made that experiment. As stated in the documentation of fix deform, there are two issues :

  1. it is computationally more expensive (orders of magnitudes more).
  2. More importantly, it crashes.
    ERROR on proc 7: Bond atoms 94221 94222 missing on proc 7 at step 211555 (…/ntopo_bond_all.cpp:59)

Probably, this option is more feasible for start-up shear or something like that.

Furthermore, people have been using the flip set to yes, and having no such issues, so there has to be something else. Maybe there are more restrictions on the tilt factors? Again the simulation box dimensions are quite large; ~ 70Rg,8Rg,8Rg, where Rg is Rg at equilibrium,

Okay, you may be running into LAMMPS’s compile-dependent image flag size limitations.

On the default setting, image flags may silently under- or overflow when they leave the range [-512,511]. This would correspond with your simulation running fine for millions of steps – and still running fine afterwards, but image-flag-specific analyses going wrong. This would also correspond with higher shear rates leading to errors earlier (since that means roughly the same distance is traversed before the error occurs).

You can try to avoid this by either manually “re-imaging” between runs, or recompiling LAMMPS with the CMake flag -D LAMMPS_SIZES=bigbig, and see if that solves the issue.


Thank you! The image flags were creating the issue. I recompiled LAMMPS with the flag you suggested, and everything runs smoothly.
I would suggest that this is added to the fix deform documentation. It will save a couple of weeks from somebody that stumbles upon the same issue.

Also, the density I had was not 0.85, but 0.8. This is the reason the viscosity values were a bit off. Choosing Nose-Hoover or DPD cannot result in such differences. Just mentioning it so that it is recorded somewhere I guess.

Thank you again!