Shear Flow - Fix deform - SLLOD - PBC issue

[alex3]

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)

[akohlmey]

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.

[alex3]

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!

[srtee]

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.

[akohlmey]

@alex3

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

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

[alex3]

Thank you for your responses!

Shern,
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

# RADIUS OF GURATION TENSOR FOR TYPE POLY
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.

Axel,
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,

[srtee]

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.

[alex3]

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!

[zx108s221]

Hi Alex,

I am learning how to do the viscosity-shear rate simulation recently.
My problem is that I dont know how to get the right viscosity under a certain shear rate.
I following the LAMMPS example of NEMD via fix deform and fix nvt/sllod (VISCOSITY folder).
Here is the key part of my script. The simulaion is done with 900 water molecules.

# deformation
timestep        ${dtr}
change_box      all triclinic
kspace_style    pppm/tip4p 1.0e-6

variable        srate equal 1e2*${fs2s}  # strain rate from s^-1 to fs^-1
variable        vs equal -pxy*${atm2Pa}/${srate}*${fs2s}  # viscosity in Pas

fix             1 all nvt/sllod temp ${T} ${T} ${Tdampn_r}
fix             2 all deform 1 xy erate ${srate} units box remap v

thermo_style    custom step temp press ke pe etotal v_vs

run             50000

I tried to get the viscosity-shear rate curve of liquid water first, but the output is like

   Step          Temp          Press          KinEng         PotEng         TotEng          v_vs     
    200000   303.71308      222.5452       1628.6545     -10198.22      -8569.5658      199766.31    
    210000   299.58122     -438.96819      1606.4975     -10303.893     -8697.3959      193449.56    
    220000   296.2246      -632.80357      1588.4977     -10276.765     -8688.2669      66770.096    
    230000   304.69467     -510.86307      1633.9183     -10199.677     -8565.7592      175896.56    
    240000   298.49909      237.44175      1600.6946     -10202.469     -8601.7741      730028.86    
    250000   311.07246      288.6201       1668.119      -10300.535     -8632.4159      70681.944  

As you can see, the variable vs is no near to 0.001Pas at all.
Could you please give me any advice on how to apply the shear rate and calculate viscosity properly?
Thanks in advance.

[srtee]

From your pasted snippet, we cannot make any precise inferences (we can’t even say if fs2s and atm2pa have the correct numerical values).

[zx108s221]

Thanks for answering. Here is the full script.

units           real
dimension       3
boundary        p p p 
atom_style      full

#----------Convert LAMMPS real units to SI units----------

variable        kB equal 1.3806504e-23
variable        NA equal 6.02214e23
variable        kcal2J equal 4184
variable        A02m equal 1.0e-10
variable        fs2s equal 1.0e-15
variable        atm2Pa equal 101325

#----------Initial Data settings----------

read_data       900H2O.data
set             type 1 charge 0.5564
set             type 2 charge -1.1128

#----------Force field parameters----------

pair_style      lj/cut/tip4p/long 2 1 1 1 0.1546 12.0 11.0
pair_modify     mix arithmetic tail yes
bond_style      harmonic
angle_style     harmonic
dihedral_style  charmm
kspace_style    pppm/tip4p 1.0e-6
neighbor        2.0 bin
neigh_modify    every 1 delay 10 check yes

# TIP4P-2005 H2O H=1,O=2
# https://doi.org/10.1063/1.2121687
pair_coeff      1 1 0.0 0.0             # H-H
pair_coeff      2 2 0.1852 3.1589       # O-O
bond_coeff      1 0 0.9572              # H-O
angle_coeff     1 0 104.52              # H-O-H

#----------Group----------

group           1 type 2                # Oh2o

#----------Calculations----------

variable        i equal 0
variable        ransed equal 19930206+${i}
variable        T equal 300                             # temp in Kelvin
variable        P equal 1                               # pressure in atm
variable        dte equal 1.0                           # equilibrium
variable        dtr equal 1.0                           # running data
variable        si equal 10                             # sample interval
variable        ci equal 1000                           # correlation interval
variable        ti equal ${si}*${ci}                    # total interval
variable        Nevery equal 10                         # Nevery
variable        Nrepeat equal 1000                      # Nrepeat
variable        Nfreq equal ${Nevery}*${Nrepeat}        # Nfreq
variable        th equal 10000                          # thermo interval
variable        Tdampb equal ${dte}*10
variable        Tdampc equal ${dtr}*1000
variable        Tdampn_e equal ${dte}*200
variable        Tdampn_r equal ${dtr}*200
variable        Pdampn_e equal ${dte}*2000
variable        Pdampn_r equal ${dtr}*2000

#----------Log----------

#log             log${i}-${T}K_${P}atm.lammps

#----------Run----------

timestep        ${dte}
thermo          ${th}
thermo_style    custom step temp press ke pe etotal density

velocity        all create ${T} ${ransed} rot yes dist gaussian

fix             shakeh2o all shake 1.0e-6 100 0 b 1 a 1

fix             b all temp/berendsen ${T} ${T} ${Tdampb} 
fix             e all nve

run             100000

unfix           b
unfix           e

fix             1 all npt temp ${T} ${T} ${Tdampn_e} iso ${P} ${P} ${Pdampn_e}

run             100000

unfix           1

# deformation
timestep        ${dtr}
change_box      all triclinic
kspace_style    pppm/tip4p 1.0e-6

variable        srate equal 1e2*${fs2s}  # strain rate from s^-1 to fs^-1
variable        vs equal -pxy*${atm2Pa}/${srate}*${fs2s}  # viscosity in Pas

thermo_style    custom step temp press ke pe etotal v_vs

fix             1 all nvt/sllod temp ${T} ${T} ${Tdampn_r}
fix             2 all deform 1 xy erate ${srate} units box remap v

fix             3 all ave/time ${si} ${ci} ${ti} v_vs ave running file vs-${i}.log

run             100000

My concern lies in two ways.
First, is the deformation process and code correct?
Second, pxy for such a small liquid simulation system dont have a stable value around 1atm. So is it possible to calculate shear viscosity using the formula ‘-pxy/srate’? If isnot, is it a good choice to just calculate shear viscosity as the usual ways(e.g. Green-Kubo method for pxy)?

Thanks!

[Germain]

Hi @zx108s221,

I am passing by to say that this question is off topic and that I would advise you to generally avoid necroposting (see the forum guidelines for more info).

You are more likely to get relevant answer by 1. posting in the relevant category (Science Talk seems to be the most relevant, it might fit the use of the LAMMPS beginner tag in LAMMPS but this is not really a LAMMPS issue per-se), 2. make a new post detailing your issue in its title.

But this question is a methodology question and it would be wiser to dig the relevant scientific literature. The Todd and Daivis book referenced in the LAMMPS sllod manual page is an excellent starting point.